1 Introduction↩︎

We consider the following stochastic differential equation \[\label{main95equation} \left\{ \begin{array}{ll} \displaystyle{ \,{\mathrm d}X(t) = a(t,X(t-))\,{\mathrm d}t + \sigma \,{\mathrm d}L(t), \;t\in [0,T]},\\ X(0)=x_0, \end{array} \right.\tag{1}\] where \(x_0\in\mathbb{R}^d\), \(\sigma\in\mathbb{R}^{d\times m}\), the drift coefficient \(a:[0,T]\times\mathbb{R}^d\to\mathbb{R}^d\) is assumed to be at least continuous, \((L(t))_{t\geq 0}\) denotes an \(m\)-dimensional Lévy process, and \(X(t-)=\lim\limits_{s\to t, s<t}X(s)\), see [1]. We investigate a non-standard approach to approximating solutions of such equations. In particular, we propose a novel deep learning-based methodology that enables artificial neural networks to learn and replicate the trajectories of the underlying SDE. Hence, the main question we pose in this paper is as follows:

Can artificial neural networks learn the stochastic dynamics given by the SDE 1 ?

We know the answer is positive for deterministic ODEs and PDEs, see, for example, [2], [3]. However, learning stochastic dynamics appears to be a much more involved task. The main problem is that artificial neural networks are deterministic functions. Therefore, we must be precise about exactly what ANNs are approximating when learning solutions to SDEs. The direct approach, known from the ODE case, is not possible for SDE. In our approach, we will use the fact that the solution \(X\) lives on the path space \(D([0,T],\mathbb{R}^m)\) - the Skorokhod space - and can be expressed as a deterministic functional \(\Psi\) of the trajectories of the driving process \(L\), see [4]. This fact allows us to define a suitable loss function for artificial neural network of a feedforward type. After minimizing the loss (i.e., training), we obtain the artificial neural network that approximates the function \(\Psi\). As a byproduct, we get a new approximation algorithm for the SDEs 1 driven by the additive Lévy noise.

One of the first attempts in that direction (solving SDEs via artificial neural networks) was [5]. The thesis [5] introduces a neural network-based discretization of the Wong–Zakai approximation for solving SDEs. In this approach, the Brownian motion is replaced by a smooth interpolation, transforming the SDE into a deterministic ODE system. A deep neural network is trained to approximate the ODE solutions by minimizing a loss function reflecting the SDE’s drift–diffusion dynamics. Once trained, the model can generate approximate trajectories of the underlying stochastic process, offering a data-driven alternative to classical numerical schemes. Numerical experiments show that while the approach successfully reproduces pathwise dynamics, it currently underperforms classical schemes such as Euler–Maruyama and Milstein in terms of convergence rate and computational efficiency, indicating the need for further methodological refinement. In [6], the proposed method solves stochastic partial differential equations (SPDEs) by combining the Wiener chaos expansion with neural network approximation. The Wiener chaos expansion of the solution is first truncated and deterministic or random neural networks replace the propagators (its deterministic coefficient functions). Deterministic networks are fully trained, while random networks use fixed random hidden-layer parameters, requiring only the linear output layer to be optimized. By leveraging the universal approximation property of neural networks, this framework provides theoretical convergence guarantees and practical algorithms for approximating a wide class of SPDEs. Finally, the article [7] introduces SPINODE, a stochastic physics-informed neural ordinary differential equation framework designed to learn unknown hidden physics terms inside stochastic differential equations from trajectory data. It represents the hidden physics via neural networks and propagates uncertainty through the known SDE structure to obtain deterministic ODEs governing the time evolution of statistical moments (e.g., means and covariances). Using ODE solvers within the neural ODE paradigm and training via moment-matching with mini-batch gradient descent and adjoint sensitivity, the method fits the neural networks to data-estimated moments. They conclude that SPINODE is flexible and scalable for multivariate stochastic systems with multiplicative noise, but note that learning requires many repeated trajectories.

In this paper, we propose a different approach to that known from [5], [6], [7]. Namely, instead of replacing Brownian motion \(W\) (or, in our case the Lévy process) with some approximation, we transform our semimartingale equation 1 into a suitable random ODE, with a noise process embedded in the right-hand side function. This methodology allows us to cover the SDEs driven with additive Lévy noise and is more direct than the method proposed in the papers mentioned above. Compared to the existing literature, we extend the PINN approach (see, for example, [8], [3]) to define its stochastic counterpart, called stochastic physics-informed neural networks (StPINNs). We provide a suitable theoretical framework and we present results of numerical experiments.

2 Preliminaries↩︎

For \(x\in\mathbb{R}^d\) by \(\|x\|\) we mean the Euclidean norm, while for a matrix \(A\in\mathbb{R}^{d\times m}\) by \(\|A\|\) we mean the Frobenius norm. Let \((\Omega,\mathcal{F},\mathbb{P})\) be a complete probabilistic space, and let \((\mathcal{F}_t)_{t\geq 0}\) be a filtration on that space that satisfies the usual conditions, see [9]. We set \(\mathcal{F}_{\infty}=\sigma\Bigl(\bigcup_{t\geq 0}\mathcal{F}_t\Bigr)\). For two sub-\(\sigma\)-fields \(\mathcal{A}_1,\mathcal{A}_2\) of \(\mathcal{F}\) we define \(\mathcal{A}_1\vee\mathcal{A}_2=\sigma(\mathcal{A}_1\cup\mathcal{A}_2)\). We assume that \((L(t))_{t\in [0,T]}\) is an \(m\)-dimensional Lévy process on \((\Omega,\mathcal{F},\mathbb{P})\) and with respect to the filtration \((\mathcal{F}_t)_{t\geq 0}\), see Section 4, Chapter I in [9]. By \(\mathcal{D}^1\left([0,T]\times\Omega;\mathbb{R}^d\right)\) we denote the space of all \(\mathbb{R}^d\)-valued stochastic processes \((y(t))_{t\in [0,T]}\) that are adapted to \((\mathcal{F}_t)_{t\geq 0}\), have a.s. continuous trajectories with \((\mathcal{F}_t)_{t\geq 0}\)-adapted and càdlàg derivative processes \((y'(t))_{t\in [0,T]}\), and have the finite norm \[\|y\|_1=\Bigl(\mathbb{E}\Bigl[\sup\limits_{0\leq t\leq T}\|y(t)\|^2\Bigr]\Bigr)^{1/2}+\Bigl(\mathbb{E}\int\limits_0^T\|y'(t)\|^2dt\Bigr)^{1/2}<+\infty.\] (See also Lemma 2 in Appendix.)
We impose the classical assumptions for the drift coefficient \(a:[0,T]\times\mathbb{R}^d\to\mathbb{R}^d\):

• \(a\in C\left([0,T]\times\mathbb{R}^d;\mathbb{R}^d\right)\)

• there exists \(L_a\in [0,+\infty)\) such that for all \(t\in [0,T],y_1,y_2\in\mathbb{R}^d\) \[\|a(t,y_1)-a(t,y_2)\|\leq L_a\|y_1-y_2\|.\]

Since the Lévy process \((L(t))_{t\in [0,T]}\) is a semimartingale (see Corollary at page 55 in [9]), by Theorem 6, page 255 in [9], under the assumptions (A1), (A2) the equation 1 has a unique strong solution \(X\). Moreover, since \(X\) has càdlàg paths, we get for all \(t\in [0,T]\) \[\int\limits_0^t a(s,X(s-))ds=\int\limits_0^t a(s,X(s))ds \quad a.s.\] Hence, instead of 1 in the rest of the paper we will be considering \[\label{main95equation2} \left\{ \begin{array}{ll} \displaystyle{ \,{\mathrm d}X(t) = a(t,X(t))\,{\mathrm d}t + \sigma \,{\mathrm d}L(t), \;t\in [0,T]},\\ X(0)=x_0, \end{array} \right.\tag{2}\]

3 Transforming SDE to RODE↩︎

We introduce the process \[\label{transform951} Y(t)=X(t)-\sigma L(t), \;t\in [0,T],\tag{3}\] and by 2 with probability one for all \(t\in [0,T]\) \[Y(t)=x_0+\int\limits_0^t a(s,Y(s)+\sigma L(s))ds.\] Hence, the stochastic process \((Y(t))_{t\in [0,T]}\) solves the following random ordinary differential equation (RODE) \[\label{main95equation23} \displaystyle{ \mathbb{P}\Bigl(\forall_{t\in [0,T)}Y'(t) = f(t,Y(t),L(t)), Y(0)=x_0\Bigr)=1},\tag{4}\] with \[f(t,y,w)=a(t,y+\sigma w), \quad (t,y,w)\in [0,T]\times\mathbb{R}^d\times\mathbb{R}^m,\] and where for \(t\in[0,T)\) by \(Y'(t)\) we mean the right-hand side derivative. The value of \(Y'(T)\) can be defined arbitrarily. Here we take \(Y'(T):=Y'(T-)\) a.s., since \(Y'(T-)=f(T,Y(T),L(T-))\) is a.s. finite. The process \((Y(t))_{t\in [0,T]}\) is adapted, has continuous trajectories, the derivative process \((Y'(t))_{t\in [0,T]}\) has càdlàg paths, and is also adapted. Note that the trajectories of \(Y\) are, in general, of higher smoothness than the trajectories of \(X\). For example, if \(L\) is the Wiener process, then the trajectories of \(X\) are nowhere differentiable, while the trajectories of \(Y\) are even continuously differentiable on \([0,T]\) a.s.

The proof of the following fact is straightforward and therefore omitted.

From Theorem 3.1 in [4] and Fact 1, we get that there exists a Skorokhod measurable function \(\Psi:D([0,T],\mathbb{R}^m)\to D([0,T],\mathbb{R}^d)\) such that \[\label{univ95rep951} Y=\Psi(L),\tag{5}\] where \(Y\) is the unique solution to 4 . This give us a pathwise representation of the solution \(Y\) on the path space, i.e.: \[Y(\omega)=\Psi(L(\omega)),\] where, for a fixed \(\omega\in\Omega\), \(L(\omega)\in D([0,T],\mathbb{R}^m)\) is a fixed trajectory of the Lévy process \(L\). Moreover, from 3 we have that \[X(\omega)=\bar\Psi(L(\omega)):=\Psi(L(\omega))+\sigma L(\omega),\] where \(\bar\Psi:D([0,T],\mathbb{R}^m)\to D([0,T],\mathbb{R}^d)\) is also Skorokhod measurable. Therefore,

the main idea is to use a deep neural network \(\mathcal{N}(\mathrm{w},\cdot)\), which input will be a finite approximation of trajectories of the Lévy process \(L\), to approximate \(\Psi\).

In this paper, we focus on the feedforward architecture for \(\mathcal{N}(\mathrm{w},\cdot)\). However, other approaches are also possible. (This will be a topic of our future work.)

4 Reformulation of solving SDE as a stochastic optimization problem↩︎

4.1 Theoretical loss function↩︎

Under assumption (A1), we define the nonlinear functional \[\mathcal{\bar L}:\mathcal{D}^1\left([0,T]\times\Omega;\mathbb{R}^d\right)\to [0,+\infty],\] by \[\mathcal{\bar L}(y)=\mathbb{E}\Bigl[H(y,\tau,L)\Bigr], \;y\in \mathcal{D}^1\left([0,T]\times\Omega;\mathbb{R}^d\right), \label{eq:loss95function}\tag{6}\] where \[H(y,\tau,L)=\|x_0-y(0)\|^2+T\cdot\|y'(\tau)-f(\tau,y(\tau),L(\tau))\|^2,\] and the random variable \(\tau\sim U(0,T)\), defined on \((\Omega,\mathcal{F},\mathbb{P})\), is independent \(\mathcal{F}_{\infty}\). We refer to \(\mathcal{\bar L}\) as the theoretical loss function.

Below, we state sufficient conditions for the finiteness of \(\mathcal{\bar L}.\)

Below, we prove further properties of the theoretical loss function \(\mathcal{\bar L}\).

From the coercivity we have that \(\mathcal{\bar L}(y)+\infty\) as \(\|y\|_1\to+\infty\). By the consistency property, we have that the sequence of approximations \((y_n)_{n\in\mathbb{N}}\) converges to the solution \(Y\) iff \(\mathcal{\bar L}(y_n)\to 0\) as \(n\to+\infty\).

4.2 Optimization problem statement↩︎

By Proposition 1, Lemma 1, and Corollary 2, we arrive at the following result for the theoretical loss function \(\mathcal{\bar L}\).

By the results above, we can express the solution to 4 in terms of finding a global minimizer of \(\mathcal{\bar L}\). Note, however, that the corresponding minimization problem \[\label{inf95dim95rodes} \inf\limits_{u\in \mathcal{D}^1\left([0,T]\times\Omega;\mathbb{R}^d\right)}\mathcal{\bar L}(u),\tag{10}\] is infinite-dimensional. Hence, for the practical implementation, instead of considering 10 , we aim at the following finite-dimensional minimization problem \[\inf\limits_{\mathrm{w}\in\mathbb{R}^N}\mathcal{L}_n(\mathrm{w}),\] where, for a given \(n\in\mathbb{N}\), \[\mathcal{L}_n(\mathrm{w})=\mathbb{E}\Bigl[H_n(\mathcal{N}(\mathrm{w},\cdot),\tau,L)\Bigr],\] \[\begin{align} \label{eq:finte95loss95function} &H_n(\mathcal{N}(\mathrm{w},\cdot),\tau,L)=\|\mathcal{N}(\mathrm{w},0,\bar L_{\Delta_n})-x_0\|^2\notag\\ &+T\cdot\Bigl\|\frac{\partial}{\partial t}\mathcal{N}(\mathrm{w},t,\bar L_{\Delta_n})\Bigl|_{t=\tau}-f(\tau,\mathcal{N}(\mathrm{w},\tau,\bar L_{\Delta_n}),\tilde{L}_n(\tau))\Bigl\|^2, \end{align}\tag{11}\] the function \(\mathcal{N}:\mathbb{R}^N\times [0,T]\times\mathbb{R}^{m\times n}\to\mathbb{R}^d\) is the feedforward artificial neural network (see, for example, [3]) of sufficient regularity (see Theorem 1), \[\Delta_n=\{0=t_0<t_1<\ldots<t_n=T\},\] \[\bar L_{\Delta_n}=\{L(t_1),L(t_2),\ldots, L(t_n)\},\] and \((\tilde{L}_n(t))_{t\in [0,T]}\) is:

• either the Lévy bridge approximation of \(L\), i.e.: \(\tilde{L}_n(t)\) is a sample drawn from the conditional law of \(L(t) \;| \;\sigma(\bar L_{\Delta_n})\) see [14],

• or the step process based on \(\bar L_{\Delta_n}\), i.e., \[\tilde{L}_n(t)=\sum\limits_{j=0}^{n-1}L(t_j)\mathbf{1}_{[t_j,t_{j+1})}(t)+L(T)\mathbf{1}_{\{t_n\}}(t).\]

Note that, for general Lévy processes, conditional distributions and simulation of Lévy bridges can be highly nontrivial, especially in the cases of infinite jumps activity (infinite variation). The exact Lévy bridge simulation is possible, for example, in the case of the Wiener process, and the (compound) Poisson process. Despite of the mentioned drawbacks, the exact simulation \(\tilde{L}_n\) from the Lévy bridge has the following advantage \[\label{popul95loss95prop}\mathcal{L}_n(\mathrm{w})=\mathcal{\bar L}(y_{n,\mathrm{w}})\tag{12}\] where \(y_{n,\mathrm{w}}(t,\omega)=\mathcal{N}(\mathrm{w},t,\bar L_{\Delta_n}(\omega))\), see the following fact.

The adaptedness of \(y_{n,\mathrm{w}}\) can be provided by enforcing causality in the neural network architecture by, for example, plugins the inputs \(L(t_j)\) up to the current time \(t\in [t_j,t_{j+1})\) and filling the rest of the inputs (for \(t_j> t\)) with zeros.

We refer to \(\mathcal{L}_n\) as to approximated loss function.

4.3 SGD algorithm↩︎

For fixed \(n\in\mathbb{N}\), by applying the Robbins-Monro stochastic approximation algorithm [15], we get \[\begin{align} \label{SGD95H95n} &\mathrm{w}^0\in\mathbb{R}^N,\notag\\ &\mathrm{w}^{k+1}=\mathrm{w}^k-\eta_k\cdot\nabla_{\mathrm{w}} H_n(\mathcal{N}(\mathrm{w},\cdot),\tau_k,L_k)\Bigl|_{\mathrm{w}=\mathrm{w}^k} \end{align}\tag{13}\] for \(k=0,1,2,\ldots\), where \((\tau_k)_{k\geq 0}\) is an iid sequence from \(U(0,T)\), \((L_k)_{k\geq 0}\) is a sequence of independent Lévy processes (that is also independent of \((\tau_k)_{k\geq 0}\)), and we impose the following assumption on the sequence of learning rates \(\{\eta_k\}_{k\geq 0}\) \[\sum\limits_{k\geq 0}\eta_k=+\infty, \;\sum\limits_{k\geq 0}\eta^2_k<+\infty.\] After training the neural network and obtaining \(\textrm{w}^*\), we get the neural network such that \[\mathcal{N}(\mathrm{w}^*,t,\bar L_{\Delta_n}(\omega))\approx (\Psi(L(\omega)))(t),\] and we define the approximation to \(X\) as \[\bar X_n(t,\omega)=\mathcal{N}(\mathrm{w}^*,t,\bar L_{\Delta_n}(\omega))+\sigma \tilde{L}_n(t,\omega), \quad t\in [0,T], \omega\in\Omega.\]

The resulting scheme, together with the above derivations, introduces the so-called Stochastic Physics-Informed Neural Networks (StPINNs). We have the following Céa-type quasi-optimality result.

Note that when \(\sigma=0\), then we recover the classical PINNs, see [3]. In this special case the neural network that we train does not depend on \(\bar L_{\Delta_n}\), so \(\mathcal{N}:\mathbb{R}^N\times [0,T]\to\mathbb{R}^d\). Moreover, \(X=Y\in C^1([0,T];\mathbb{R}^d)\), \(\|X\|_1=\sup\limits_{0\leq t\leq T}\|X(t)\|+\Bigl(\int\limits_0^T\|X'(t)\|^2dt\Bigr)^{1/2}\) and the approximate loss \(\mathcal{L}(\mathrm{w})=\mathcal{\bar L}(\mathcal{N}(\mathrm{w},\cdot))\) does not depend on \(n\). Hence, we have the following corollary.

The main differences between the derivation of StPINNs and PINNs are the use of the universal representation theorem that gives us the existence of the function \(\Psi\) in 5 , and the construction of the loss function for the artificial neural network \(\mathcal{N}(\mathrm{w},\cdot,\cdot)\) that, after minimization, allows us to approximate the function \(\Psi\). It was not the case in the classical PINNs.

4.4 Some insight into the multiplicative noise case - Doss-Sussman transformation↩︎

Our diffusion removal transformation is only suitable for the SDEs with the additive noise as in 1 and cannot be directly applied to the SDEs with the multiplicative noise. In the scalar case and \(L=W\) being a scalar Wiener process let us consider the SDE \[\label{scalar95m95equation4} \left\{ \begin{array}{ll} \displaystyle{ \,{\mathrm d}X(t) = a(t,X(t))\,{\mathrm d}t + \sigma X(t)dW(t), \;t\in [0,T]},\\ X(0)=x_0>0, \end{array} \right.\tag{14}\] To transform the above equations into the RODE of the form 4 we proceed as follows. First, consider the SDE \[\label{scalar95m95equation5} \left\{ \begin{array}{ll} \displaystyle{ \,{\mathrm d}Z(t) = (e^{-Z(t)}a(t,e^{Z(t)})-\sigma^2/2)\,{\mathrm d}t + \sigma dW(t), \;t\in [0,T]},\\ Z(0)=\ln(x_0), \end{array} \right.\tag{15}\] and then use the transformation 3 , i.e., \(Y(t)=Z(t)-\sigma W(t)\), to transform the SDE 15 to the RODE \[\label{main95equation231} \displaystyle{ \mathbb{P}\Bigl(\forall_{t\in [0,T)}Y'(t) = f(t,Y(t),W(t)), Y(0)=\ln(x_0)\Bigr)=1},\tag{16}\] where \[f(t,y,w)=e^{-(y+\sigma w)}a(t,e^{y+\sigma w})-\sigma^2/2, \quad (t,y,w)\in [0,T]\times\mathbb{R}\times\mathbb{R}.\] Then, we can apply our StPINN method to approximate \(Y\). Since \(X(t)=e^{Z(t)}\), what can be verified by the Itô formula, finally, we get the approximation of \(X\).

Note that both transformations are the particular case of the well-known Doss-Sussman transformation, see [16], [17]. In our future work we investigate the possibility of applying the D-S transformation to the multidimensional SDEs with multiplicative noise to extend the StPINNs approach derived in this paper.

5 Numerical experiments↩︎

In this section, we focus on numerical experiments. All the code is available at
https://github.com/MarcinBaranek/Research/tree/master/SDE-ANN-solver
We use a discretization of the time interval \([0,T]\) as a uniform mesh with \(n\) elements. For numerical purpose we take \(L\) as the Wiener process. As the exact solution \(X\) of 1 we treat the solution from the Euler algorithm calculated with precision \(float64\) and the number of discretization points equal to \(2^{17}\). By \(\bar{X}_n\) we mean the solution calculated by the neural network with precision \(float32\) and information about \(L\) at \(n\) points. To see how perform the network we define errors \[err_{2,n}([0,T]) = \sqrt{\frac{1}{M}\sum_{i=1}^M \frac{T}{n}\sum_{j=1}^n\vert \bar{X}(t_j, L_i) - X(t_j,L_i)\vert^2}\approx\Bigl(\mathbb{E}\int\limits_0^T|\bar X_n(s)-X(s)|^2ds\Bigr)^{1/2},\] and the error at the end-point \(T\) is approximated by \[err_{2,n}(T) = \sqrt{\frac{1}{M}\sum_{i=1}^M \vert \bar{X}(T, L_i) - X(T,L_i)\vert^2}\approx (\mathbb{E}|\bar X_n(T)-X(T)|^2)^{1/2}.\] Those errors are calculated with precision \(float64\). We take \(M=10000\), unless specified otherwise.

We consider three exemplary equations, where for the first two \(L\) is a scalar Wiener process, and for the last example \(L\) is a four dimensional Wiener process.

• \[\label{eq:example1} \left\{ \begin{array}{ll} \displaystyle{ \,{\mathrm d}X(t) = 5(0.4-X(t-))\,{\mathrm d}t + 0.61 \,{\mathrm d}L(t), \;t\in [0,1]},\\ X(0)=-0.3. \end{array} \right.\tag{17}\]

• \[\label{eq:example2} \left\{ \begin{array}{ll} \displaystyle{ \,{\mathrm d}X(t) = 5(0.4-\sin(X(t-)))\,{\mathrm d}t + 0.61 \,{\mathrm d}L(t), \;t\in [0,1]},\\ X(0)=-0.3. \end{array} \right.\tag{18}\]

• \[\left\{ \begin{array}{ll} \displaystyle{ dX_t = 5\left(\begin{bmatrix}0.4\\0.7\\-1.2\end{bmatrix} - X_t\right) \, dt + \begin{bmatrix} 0.61&0.3&0&0.1\\ -0.3&-0.72&0.45&0\\ 0&0&0.34&-0.81 \end{bmatrix} \, dL_t,\quad t\in[0,1]},\\ X(0)= \begin{bmatrix} -0.3\\-0.3\\-0.3 \end{bmatrix}. \end{array} \right.\]

5.1 Implementation details↩︎

We employ a fully connected deep neural network with three hidden layers, using \(\tanh\) as the activation function. A sample implementation of this network is shown in Listing [lst:network]. We assume that the network input has shape \((\mathrm{Batch\;size}, 1+n)\), where \(n\) is the number of points in a uniform mesh over the interval \([0, T]\). The output of the network has shape \((\mathrm{Batch\;size}, 1)\).

import keras

class Network(keras.Model):
    def __init__(
            self, n_points: int, initial_value: float,
            derivative_initial_value: float
    ):
        super().__init__()
        hidden_layers = [
            keras.layers.Dense(n_points, activation="hard_silu"),
            keras.layers.Dense(512, activation="hard_silu"),
            keras.layers.Dense(256, activation="hard_silu"),
            keras.layers.Dense(128, activation="hard_silu"),
        ]
        self.initial_value = initial_value
        self.derivative_initial_value = derivative_initial_value
        self.hidden_layers = hidden_layers
        self.out = keras.layers.Dense(1)

    def call(self, x):
        t = x[..., :1]
        h = x
        for layer in self.hidden_layers:
            h = layer(h)
        f_theta = self.out(h)
        return self.initial_value\;
            + self.derivative_initial_value * t\;
            + f_theta * t ** 2 / 2.

The neural network \(\mathcal{N}\) is trained to minimize the loss function defined in 11 . To accelerate the training process and enforce the initial condition, we constrain the network to satisfy \[\mathcal{N}(\mathrm{w},0,\bar{L}_{\Delta_n})=x_0.\] This constraint is imposed through the transformation \[\bar{\mathcal{N}}(\mathrm{w},t,\bar{L}_{\Delta_n})=x_0 + t\mathcal{N}(\mathrm{w},t,\bar{L}_{\Delta_n}).\] which guarantees that the initial condition is satisfied identically. As a consequence, the first term in 11 is identically zero for all network parameters. Furthermore, prior to training, the expected value of the time derivative of the solution at \(t=0\) is known and equals \(f(0,x_0,0)\). By combining this information with the imposed initial condition, the neural network architecture can be expressed as \[\bar{\mathcal{N}}(\mathrm{w},t,\bar{L}_{\Delta_n})=x_0 +tf(0,x_0,0) +\frac{t^2}{2}\mathcal{N}(\mathrm{w},t,\bar{L}_{\Delta_n}).\] This behavior is implemented in lines 31–33 of Listing [lst:network].

In Figure 1 we present an example of network architecture. The number of neurons in each dense layer is our guess, and we believe the similar result are obtainable with different architecture.

5.2 training details.↩︎

We train the network using randomly sampled trajectories of \(L\) at \(n\) points, i.e., \(\bar L_{\Delta_n}\), together with independently generated time values \(\tau\). The loss function is defined in Equation 6 . In our setup, one epoch corresponds to training the network on a batch of 64 sampled trajectories of \(L\).

In Figure 2, we observe the evolution of the training loss and the approximated errors. The approximation errors are plotted on a logarithmic scale. We hypothesize that the neural network is unlikely to overfit, provided that we use randomly sampled trajectories of \(L\), therefor the network is unable to memorize the training set. The hypothesis is supported empirically by Figure 2. The training loss stabilizes over time, and the approximation error exhibits a similar behavior. Consequently, no increase in the approximation error is observed during training.

The network \(\mathcal{N}\) appears to be overly influenced by the trajectories \(\bar{L}_{\Delta_n}\) and struggles to solve the ODE (case when \((\bar{L}_{\Delta_n} \equiv 0)\)). The input contains a lot of numbers representing \(\bar{L}_{\Delta_n}\) and only one representing the time, therefore the network may doesn’t learn how important the time is. The example solution is in Figure 3 (a). To make the network paying a more attention to the time, we train it with randomly generated trajectories \(\bar{L}_{\Delta_n}\), perform weights optimization, train with trivial trajectories \(\bar{L}_{\Delta_n}\equiv0\), perform weights optimization. The result in the ODE case are presenting in Figure 3 (b).

We observed improved performance; however, there remains room for further optimization. Notably, the network struggles to accurately approximate the solution toward the end of the interval. To address this issue, we introduced a weighting factor by multiplying the loss function by \(\tau\), thereby increasing the penalty as \(\tau\) grows. The effects of this modification are illustrated in Figure 4.

5.3 Results↩︎

This subsection focuses on visualizing the model’s performance. We examine the trajectories generated based on the model output and contrast them with the corresponding exact trajectories. This comparison offers insight into how well the learned dynamics replicate the system’s true behavior.

In Figure 5, we present the performance of the model trained to solve equation in example 1. The green trajectory comes from the exact solution, and the blue is calculated by the neural network. The graphs with trajectories are entitled with values of \(err_{2,512}([0,T])\) and \(err_{2,512}(T)\), but in this case, the errors are calculated for \(M=1\).

The next chart (Fig. 6) presents the performance of the model trained to solve example 2. The results indicate that the model’s performance is highly similar to that obtained for example 1, despite the presence of a nonlinear drift.

Additionally, we examine the behavior of the neural network for Example 1, using as the process \(L\) both the Poisson process and a composite process consisting of the sum of a Poisson process and a Wiener process. In Figure 7 we present sample trajectories for such selection of process \(L\). The charts on the left side shows trajectories driven by a Poisson process and charts on the right side presents trajectories driven by a sum of Poisson and Wiener processes.

We examine the behavior of the neural network for Example 2, noting that this case corresponds to a nonlinear setting. In this experiment, the driving process \(L\) is again considered in two variants: a Poisson process and a composite process given by the sum of a Poisson process and a Wiener process. In Figure 8 we present sample trajectories obtained for this configuration. The charts on the left-hand side show trajectories generated when the system is driven by a Poisson process, while the charts on the right-hand side illustrate trajectories corresponding to the composite process formed by the sum of Poisson and Wiener processes.

We examine the multidimensional Example 3 in two cases. In the first case, \(L\) is a four-dimensional Wiener process, while in the second case, it is a four-dimensional process given by the sum of a Wiener and a Poisson process. The results are presented in Figure 9.

6 Appendix↩︎

We present here auxiliary properties of the processes \(Y\).

Directly from Proposition 4 we get the following corollary.

We believe that the lemma below is well-known. However, we were unable to find a direct statement or proof, so, for the reader’s convenience, we provide its justification.

7 Conclusion and future work↩︎

In this paper, we have introduced the StPINNs, the extension of the classical PINNs. Our approach allows us to solve SDEs 1 , driven by an additive Lévy process, with artificial neural networks. The derivation performed within this paper shows that there is a rich mathematical structure behind StPINNs, and throughout mathematical investigations should be done in order to state convergence theorems, etc., rigorously. In our future work, we plan to address some of these problems and also consider the SDEs with multiplicative noise case.

References↩︎