Embedding stochastic differential equations into neural networks via
  dual processes
We propose a new approach to constructing a neural network for predicting
expectations of stochastic differential equations. The proposed method does not
need data sets of inputs and outputs; instead, the information obtained from
the time-evolution equations, i.e., the corresponding dual process, is directly
compared with the weights in the neural network. As a demonstration, we
construct neural networks for the Ornstein-Uhlenbeck process and the noisy van
der Pol system. The remarkable feature of learned networks with the proposed
method is the accuracy of inputs near the origin. Hence, it would be possible
to avoid the overfitting problem because the learned network does not depend on
training data sets.
Naoki Sugishita, Jun Ohkubo
# Embedding stochastic differential equations into neural networks via dual processes

###### Abstract

We propose a new approach to constructing a neural network for predicting expectations of stochastic differential equations. The proposed method does not need data sets of inputs and outputs; instead, the information obtained from the time-evolution equations, i.e., the corresponding dual process, is directly compared with the weights in the neural network. As a demonstration, we construct neural networks for the Ornstein-Uhlenbeck process and the noisy van der Pol system. The remarkable feature of learned networks with the proposed method is the accuracy of inputs near the origin. Hence, it would be possible to avoid the overfitting problem because the learned network does not depend on training data sets.

## 1 Introduction

Many researchers have studied the evolution of dynamical or stochastic systems in physics and other research fields. Recently, the amount of data has exponentially increased, and there are many studies on dynamical systems based on these large data sets. One of the aims of these studies is to transform observed data into predictive models of the physical world, and neural networks are a hopeful candidate for this aim. Of course, a simple application of conventional neural networks would not be enough because the physical world has many characteristics due to various constraints. For example, it would be beneficial to equip features of the time-evolution into the learning steps to make the prediction more accurate. Some ideas have appeared recently; the numerical scheme for time-evolution with multi-step time-stepping schemes is employed [1]. Some works focused on network architectures. Reference [2] discussed a network architecture for the time-evolution partial differential equation. In [3, 4], recurrent neural networks were applied for model predictive control. There is a study to deal with fluid flow simulation with long short-term memory (LSTM) [5]. Other types of discussions based on the universal approximation theorem were given in [6], in which DeepONet was also proposed.
The physics-informed machine learning has the same research direction. In [7, 8], frameworks for learning partial differential equations were discussed; see a review [9] for this topic. There are studies in which conservation laws are combined with learning; the inclusion of the Hamiltonian structures makes learning more stable [10]. A recent paper gives a good review of this topic in the introduction [11]. However, compared with the studies on deterministic dynamical systems, there has been little discussion about systems with noise, i.e., stochastic systems. Such stochastic systems are the topic of the present paper.

One of the problems of machine learning for physical subjects is the data sets. For example, when one constructs a neural network for prediction from an input coordinate, it is necessary to prepare data with various initial conditions. Note that stochastic systems require large data sets compared with deterministic systems; in stochastic systems, we must consider statistics of predictions, and the calculation of expected values takes high computational cost. Hence, we need more computational effort in the data preparation and the learning steps. However, if we have the information of equations governing the system as the prior knowledge, it could be possible to reduce the size of the data sets and costs for learning. Of course, it will be possible to apply additional online learning steps to the trained network to obtain a more accurate one. Hence, the key question of the present paper is as follows: How should we cooperate with the information on the time-evolution equation for stochastic systems? If we can embed the prior knowledge directly without data sets, the method would complement the conventional ones.

In the present paper, the aim is to predict statistics after time evolution in stochastic systems, especially stochastic differential equations. For this aim, we propose a new

Figure 1: Conventional and proposed methods for the learning of stochastic systems. A most naive approach is to generate data with Monte Carlo methods; target statistics are evaluated from the sampled data sets. Employing dual processes enables us to skip the sampling step; it is possible to evaluate the statistics directly by solving the backward Kolmogorov equation. These two approaches give a data set with the inputs and outputs, and then neural networks are learned with conventional backpropagation. By contrast, the proposed method employs direct comparison with the information of the dual processes.
method to learn neural networks without generating data sets for statistics. The key of the proposal is the usage of dual processes. The duality of stochastic processes has long been studied in statistical physics and mathematical physics [12], and recent studies clarified that the conventional backward Kolmogorov equation gives dual processes straightforwardly. Figure 1 summarizes the proposal. In a most naive approach, we generate data with Monte Carlo samplings and calculate target statistics for various initial conditions. As shown later, one can evaluate the target statistics directly from the dual process. In both methods, after obtaining a data set with pairs of initial conditions and the target statistics, the backpropagation gives leaned neural networks. By contrast, the proposed method employs direct comparison with the information of the dual processes. Hence, we utilize an optimization procedure instead of backpropagation. Since there is no need to evaluate statistics from the samplings, it would be possible to write that the proposed method embeds stochastic differential equations directly into neural networks. We demonstrate the proposed method with one-dimensional and two-dimensional noisy systems, which will clarify the learned features different from the conventional approach.

The remaining part of the present paper proceeds as follows. In Section 2, we review the method to evaluate statistics of the stochastic differential equations without any sampling. Section 3 gives the main proposal to construct neural networks directly from the stochastic differential equations. Two numerical demonstrations are given in section 4. Finally, section 5 concludes this paper.

## 2 Numerical method to evaluate statistics without sampling

In the present paper, we only focus on the stochastic differential equations and the statistics after the time-evolution. As for the basics of the stochastic differential equations, see [13] for example. Here, Let \(\mathbf{x}(t)\) be a \(D\)-dimensional random vector which obeys the following stochastic differential equation:

\[d\mathbf{x}=\mathbf{a}(\mathbf{x})dt+B(\mathbf{x})d\mathbf{W}(t), \tag{1}\]

where \(\mathbf{a}(\mathbf{x})\) is a vector function called the drift coefficient and \(B(\mathbf{x})\) is a matrix-valued function called the diffusion coefficient. \(\mathbf{W}(t)\) is a vector of Wiener processes.

Since one cannot generally solve stochastic differential equations analytically, the Monte Carlo method is employed to sample trajectories from stochastic differential equations. The famous method is the Euler-Maruyama method; for example, see [14]. In the Euler-Maruyama method, a time-discretization is necessary; the sampling steps are time-consuming when we want to evaluate statistics with high accuracy for various initial conditions.

As stated in the introduction, the usage of duality in stochastic processes has long been studied in physics [12]. While the usage was mainly restricted to exactly solvable cases, recent studies give a simple derivation of dual processes from stochastic differential equations [15]; there is also an algorithm based on combinatorics to evaluate
the statistics of stochastic differential equations [16]. Here, there is no need to know the details of the duality because a discussion based on the backward Kolmogorov equation is enough for the aim of the present paper. Hence, we briefly review the key points which are enough to understand the proposals in the next section.

First, we consider the following Fokker-Planck equation instead of the stochastic differential equation [13]:

\[\frac{\partial}{\partial t}p(\mathbf{x},t)=\mathcal{L}p(\mathbf{x},t), \tag{2}\]

where \(p(\mathbf{x},t)\) is the probability density function and \(\mathcal{L}\) is the time-evolution operator defined as

\[\mathcal{L}=-\sum_{i}\frac{\partial}{\partial x_{i}}a_{i}(\mathbf{x})+\frac{1}{2} \sum_{i,j}\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}\left[B(\mathbf{x},t)B( \mathbf{x})^{\mathrm{T}}\right]_{i,j}. \tag{3}\]

The aim here is to evaluate the \(m\)-th order moment of the \(i\)-th element of \(x(t)\); i.e.,

\[M_{i}^{(m)}(\mathbf{x}_{0},t)=\mathbb{E}_{p(\mathbf{x},t|\mathbf{x}_{0},t_{0})}[x_{i}^{m}], \tag{4}\]

where \(p(\mathbf{x},t|\mathbf{x}_{0},t_{0})\) is the probability density function of \(\mathbf{x}(t)\) with the initial condition \(\mathbf{x}(t_{0})=\mathbf{x}_{0}\). The initial condition is written as

\[p(\mathbf{x},t_{0})=\delta(\mathbf{x}-\mathbf{x}_{0}), \tag{5}\]

where \(\delta(\cdot)\) is the Dirac delta function.

Second, we derive the backward Kolmogorov equation instead of the Fokker-Planck equation; as for the backward Kolmogorov equation, see [17]. Since the time-evolution of the probability density function is formally denoted as

\[p(\mathbf{x},t|\mathbf{x}_{0},t_{0})=e^{\mathcal{L}(t-t_{0})}p(\mathbf{x},t_{0}), \tag{6}\]

the derivation of the time-evolution operator for the backward Kolmogorov equation is easily understood as follows:

\[M_{i}^{(m)}(\mathbf{x}_{0},t) =\int x_{i}^{m}p(\mathbf{x},t|\mathbf{x}_{0},t_{0})d\mathbf{x} \tag{7}\] \[=\int x_{i}^{m}\left\{e^{\mathcal{L}(t-t_{0})}p(\mathbf{x},t_{0}) \right\}d\mathbf{x}\] \[=\int x_{i}^{m}\left\{e^{\mathcal{L}(t-t_{0})}\delta(\mathbf{x}-\mathbf{x }_{0})\right\}d\mathbf{x}\] \[=\int\left\{e^{\mathcal{L}^{\dagger}(t-t_{0})}x_{i}^{m}\right\} \delta(\mathbf{x}-\mathbf{x}_{0})d\mathbf{x}\] \[=\int\varphi(\mathbf{x},t)\delta(\mathbf{x}-\mathbf{x}_{0})d\mathbf{x}\] \[=\varphi(\mathbf{x}_{0},t),\]

where \(\mathcal{L}^{\dagger}\) is the adjoint operator of \(\mathcal{L}\) and \(\varphi(\mathbf{x},t)\) is a solution of a time evolution equation

\[\frac{\partial}{\partial t}\varphi(\mathbf{x},t)=\mathcal{L}^{\dagger}\varphi(\bm {x},t). \tag{8}\]
The adjoint operator \(\mathcal{L}^{\dagger}\) is written as

\[\mathcal{L}^{\dagger}=\sum_{i}a_{i}(\mathbf{x})\frac{\partial}{\partial x_{i}}+\frac{ 1}{2}\sum_{i,j}\left[B(\mathbf{x})B(\mathbf{x})^{\mathrm{T}}\right]_{i,j}\frac{\partial ^{2}}{\partial x_{i}\partial x_{j}}, \tag{9}\]

which corresponds to the time-evolution operator for the backward Kolmogorov equation. In (7), we employed integration by parts and the fact that the probability density functions will vanish at \(x_{i}\to\pm\infty\). Note that the initial condition for (8) is

\[\varphi(\mathbf{x},t_{0})=x_{i}^{m}, \tag{10}\]

which corresponds to the target statistic. Here, we focus on the value of the solution of (8) at a certain coordinate \(\mathbf{x}_{0}\), i.e., \(\varphi(\mathbf{x}_{0},t)\); the value immediately gives the expectation of the target statistic, \(M_{i}^{(m)}(\mathbf{x}_{0},t)\).

Third, we expand \(\varphi(\mathbf{x},t)\) with monomial basis functions, \(\mathbf{x}^{\mathbf{n}}=\prod_{i}x_{i}^{n_{i}}\), so that

\[\varphi(\mathbf{x},t)=\sum_{\mathbf{n}\in\mathbb{N}^{D}}P_{i}^{(m)}(\mathbf{n},t)\mathbf{x}^{ \mathbf{n}}, \tag{11}\]

where \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) are the expansion coefficients. By substituting this equation into (7), we have

\[M_{i}^{(m)}(\mathbf{x}_{0},t)=\sum_{\mathbf{n}\in\mathbb{N}^{D}}P_{i}^{(m)}(\mathbf{n},t) \mathbf{x}_{0}^{\mathbf{n}}. \tag{12}\]

Hence, it is enough to obtain the coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) to evaluate the expectation \(M_{i}^{(m)}(\mathbf{x}_{0},t)\). Employing the basis expansion in (11) with the time-evolution equation in (8), we have

\[\sum_{\mathbf{n}\in\mathbb{N}^{D}}\left\{\frac{\partial}{\partial t}P_{i}^{(m)}( \mathbf{n},t)\right\}\mathbf{x}^{\mathbf{n}}=\mathcal{L}^{\dagger}\left\{\sum_{\mathbf{n}\in \mathbb{N}^{D}}P_{i}^{(m)}(\mathbf{n},t)\mathbf{x}^{\mathbf{n}}\right\}, \tag{13}\]

which gives the simultaneous ordinary differential equations for \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) by comparing the coefficients of the basis expansion. Note that the initial condition should be

\[P_{i}^{(m)}(\mathbf{n},t_{0})=\delta_{n_{i},m}, \tag{14}\]

where \(\delta_{\cdot,\cdot}\) is the Kronecker delta function. The initial condition stems from \(\varphi(\mathbf{x},t_{0})\) in (10). We will denote examples of the simultaneous ordinary differential equations for \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) later, which will help the reader understand the above discussion.

Here are some comments on the above discussion. The first comment is on duality in stochastic processes: the derived equations for \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) do not satisfy the law of conservation of probability, so the coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) are not probabilities. As shown in [18], it is possible to recover the probabilistic characteristics to extend the discussions, which leads to the connection with the duality relation in stochastic processes. However, as written later, the coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) are enough to learn neural networks; there is no need to recover the probabilistic characters here. The second comment is related to the target statistics. The above discussion is limited to the target statistics in (4), i.e., a simple moment. While the restriction simplifies the
discussion, we should mention that other statistics, such as correlations, can also be evaluated.

In summary, the moment \(M_{i}^{(m)}(\mathbf{x}_{0},t)\) is evaluated by solving the simultaneous ordinary differential equations derived from (13), without samplings of stochastic processes. Our goal is to construct a neural network to predict the target moment. The learning procedure needs expectations for various initial coordinates. Furthermore, stochastic cases need many samples for a single initial coordinate to evaluate the expected value, which requires high computational costs. Then, the above approach based on the simultaneous ordinary differential equations enables us to avoid the samplings. However, we still need to prepare a data set, pairs with an input coordinate and the target moment, to use conventional learning frameworks in previous studies. Is there a more efficient method suitable for stochastic differential equations? Next, we propose a simple way to answer this question and directly embed the information in the equations to neural networks via \(\{P_{i}^{(m)}(\mathbf{n},t)\}\).

## 3 Proposed method

Let us consider the following neural network with a single hidden layer: The input is the initial coordinate of stochastic process \(\mathbf{x}_{0}\), the output is an estimate of the moment \(M_{i}^{(m)}(\mathbf{x}_{0},t)\), the number of nodes of the hidden layer is \(n\), and the activation function is a sigmoid function \(\sigma(x)=1/(1+e^{-x})\). While this neural network is not deep, the universal approximation theorem guarantees the approximation ability when we use many hidden nodes [19, 20]. As demonstrated later, the simple structure works well for examples with nonlinear coefficients.

The key of the proposed method is the direct comparison of the coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) in (12) with the corresponding components in the neural network. Although the direct comparison is a simple idea, it yields preferable learning results, as discussed later.

First, the output of the neural network, \(y\in\mathbb{R}\), given the input \(\mathbf{x}\in\mathbb{R}^{D}\) can be represented as

\[y =\mathbf{q}^{\rm T}\mathbf{\sigma}(R\mathbf{x}+\mathbf{s}) \tag{15}\] \[=\sum_{i=1}^{n}q_{i}\sigma\left(\sum_{j=1}^{D}R_{ij}x_{j}+s_{i} \right),\]

where \(\mathbf{q},\mathbf{s}\in\mathbb{R}^{n},R\in\mathbb{R}^{n\times D}\) are the weights of the neural network. Using the Taylor expansion of \(\sigma(x)\) up to the \(N\)-th order, we can approximate the output with the power series of the input \(\mathbf{x}\) as follows:

\[y\approx\sum_{\{\mathbf{l}\in\mathbb{N}^{D}|\sum_{j}l_{j}\leq N\}} \sum_{k=\sum_{j}l_{j}}^{N}\sum_{i=1}^{n}\left(\begin{array}{c}k\\ l_{1},\ldots,l_{D},k-\sum_{j}l_{j}\end{array}\right)\frac{\sigma^{(k)}(0)}{k! }q_{i}\mathbf{r}_{i}^{l}s_{i}^{k-\sum_{j}l_{j}}\mathbf{x}^{\mathbf{l}}, \tag{16}\]

where \(\sigma^{(k)}(\cdot)\) denotes the \(k\)-th derivative of \(\sigma(\cdot)\) and \(\mathbf{r}_{i}\) is the \(i\)-th row vector of the
weight matrix \(R\). Note that a notation for multinomial coefficients

\[\left(\begin{array}{c}k\\ k_{1},\ldots,k_{r}\end{array}\right)=\frac{k!}{k_{1}!\ldots k_{r}!} \tag{17}\]

is used in (16).

At this stage, the correspondence between (12) and (16) is clear; both have the basis-expansion form for \(\mathbf{x^{n}}\) or \(\mathbf{x^{l}}\). Hence, it is possible to employ a direct comparison between them. For example, one can use the sum of the square errors of coefficients in (12) and (16) as a cost function:

\[C_{i,t}^{(m)}(\mathbf{q},R,\mathbf{s})=\sum_{\{\mathbf{l}\in\mathbb{N}^{D}|\sum_{j}l_{j} \leq N\}}\left(P_{i}^{(m)}(\mathbf{l},t)-P_{\mathrm{NN}}(\mathbf{l},\mathbf{q},R,\mathbf{s}) \right)^{2}, \tag{18}\]

where

\[P_{\mathrm{NN}}(\mathbf{l},\mathbf{q},R,\mathbf{s})=\sum_{k=\sum_{j}l_{j}}^{N}\sum_{i=1}^ {n}\left(\begin{array}{c}k\\ l_{1},\ldots,l_{D},k-\sum_{j}l_{j}\end{array}\right)\frac{\sigma^{(k)}(0)}{k! }q_{i}\mathbf{r^{l}}_{i}^{k-\sum_{j}l_{j}}. \tag{19}\]

By minimizing this cost function \(C_{i,t}^{(m)}\), we obtain the neural network that estimates the moment \(M_{i}^{(m)}(\mathbf{x}_{0},t)\).

Note that the number of simultaneous ordinary differential equations derived from (13) is infinite. Of course, we cannot evaluate the infinite number of coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\). However, with the comparison in (18), it is enough to consider a finite set of \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) for \(\left\{\mathbf{n}\in\mathbb{N}^{D}|\max_{i}n_{i}\leq N\right\}\); the order of Taylor approximation, \(N\), relates to the upper bound. From this fact, it is possible to approximate the simultaneous ordinary differential equations for \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) with the limited range of indexes \(\{\mathbf{n}\in\mathbb{N}^{D}|\max_{i}n_{i}\leq N\}\); we employ this approximation in the next section. Note that it would be natural to determine the order of Taylor approximation from the number of coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) used in the cost function in (18). The coefficients stem from the expansion of the function \(\varphi(\mathbf{x},t)\) in (11) for the backward Kolmogorov equation. Hence, the number of coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) is related to the approximation accuracy of \(\varphi(\mathbf{x},t)\). In other words, it is enough to use cutoffs that yield enough coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) of the function \(\varphi(\mathbf{x},t)\) with the accuracy one wants to approximate.

What kind of features can we expect from the proposed method? Here, we focus on the Taylor-type basis expansion in the proposed method. Hence, one could expect that the approximation performance is better when \(\mathbf{x}_{0}\) is closer to the origin. Next, we confirm this conjecture with numerical demonstrations.

## 4 Numerical examples

In this section, we demonstrate the proposed method for two stochastic differential equations. The first example is a famous one-dimensional model, i.e., the Ornstein-Uhlenbeck process [13]. The analytical solutions for moments are known for the Ornstein-Uhlenbeck process, and it is easy to compare the performance. The second
example is the noisy van der Pol system, which is a two-dimensional model. The noisy version of the original van der Pol system [21] has already been used in filtering [22, 23] and recent studies related to data analysis for dynamical systems [24].

In both cases, the simultaneous ordinary differential equations for \(P_{i}^{(m)}(\mathbf{n},t)\) are solved numerically via the functions of the Scipy package, scipy.integrate.solve_ivp. The minimization procedure for the cost function (18) is performed with scipy.optimize.least_squares, in which we randomly generate the initial values for \(\mathbf{q},R,\mathbf{s}\) from a uniform distribution with the range \([-1,1)\). The initial time \(t_{0}\) is \(0\), and the time interval for each case is denoted in table 1. Table 1 also shows the network size, \(n\), and the order of Taylor approximation, \(N\).

### Ornstein-Uhlenbeck process

The stochastic differential equation for the Ornstein-Uhlenbeck process is given as follows [13]:

\[dx=-\gamma xdt+\sigma dW(t). \tag{20}\]

Hence, the time-evolution operator \(\mathcal{L}\) for the Fokker-Planck equation in (2) is given as

\[\mathcal{L}=\frac{\partial}{\partial x}\gamma x+\frac{\partial^{2}}{\partial x ^{2}}\frac{\sigma^{2}}{2}, \tag{21}\]

which leads to the adjoint operator in (9) as

\[\mathcal{L}^{\dagger}=-\gamma x\frac{\partial}{\partial x}+\frac{\sigma^{2}}{ 2}\frac{\partial^{2}}{\partial x^{2}}. \tag{22}\]

Since there is only one variable \(x\), the coefficients are \(\{P_{i}^{(m)}(n_{1},t)\}\) for \(n_{1}=0,1,2,\cdots\). The derived simultaneous ordinary differential equations are as follows:

\[\frac{d}{dt}P_{i}^{(m)}(n_{1},t)=-\gamma n_{1}P_{i}^{(m)}(n_{1},t)+\frac{ \sigma^{2}}{2}(n_{1}+2)(n_{1}+1)P_{i}^{(m)}(n_{1}+2,t) \tag{23}\]

for \(n_{1}=0,1,2,\cdots,N\). Note that we employ the finite cutoff with \(N\), and set \(P_{i}^{(m)}(n_{1},t)=0\) for \(n_{1}>N\). Here, as denoted in Table 1, we set \(N=12\).

As for the numerical demonstrations, we construct two neural networks to estimate the first and second-order moments. As for the first-order moment case, the initial condition for (23) is set as \(P_{i}^{(m)}(n_{1},0)=\delta_{n_{1},1}\); \(P_{i}^{(m)}(n_{1},0)=\delta_{n_{1},2}\) is used for the second-order moment case.

\begin{table}
\begin{tabular}{l c c c} \hline \hline  & Time interval \(t\) & \begin{tabular}{c} Number of hidden \\ layer nodes \(n\) \\ \end{tabular} & 
\begin{tabular}{c} Order of Taylor \\ approximation \(N\) \\ \end{tabular} \\ \hline Ornstein-Uhlenbeck & 1 & 4 & 12 \\ noisy van der Pol & 0.1 & 8 & 17 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Parameter settings used in the two examples.
The analytical solutions for moments are as follows [13]:

\[M^{(1)}(x_{0},t) =x_{0}e^{-\gamma t}, \tag{24}\] \[M^{(2)}(x_{0},t) =\left(x_{0}e^{-\gamma t}\right)^{2}+\frac{\sigma^{2}}{2\gamma} \left(1-e^{-2\gamma t}\right). \tag{25}\]

Here, we set \(\gamma=\sigma=1\). Figure 2 shows the numerical results obtained from the proposed method. We also depict the analytical results. As for the proposed method, after the learning procedures in section 3, various inputs \(x_{0}=x(t=0)\) are used to obtain the corresponding outputs, and we drew the curves in figure 2.

The numerical results indicate that the first-order moment is estimated well even in the small neural network. As for the second-order moment, the accuracy decreases as the input value is further away from the origin. This behavior is just the expected one from the characteristics of the Taylor approximations, as discussed in section 3.

We here note that it is possible to change the origin. The expansion of \(\varphi(\mathbf{x},t)\) in (11) is around zero. As discussed in [25], it is sometimes beneficial to shift the origin, which yields the different values for the coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\). Then, the learned networks show good agreements around the shifted origins; we have numerically checked these characteristics. From these results, it is possible to say the idea employed in the backward Kolmogorov equation and the optimization work well.

The Ornstein-Uhlenbeck process is linear. Next, we use a two-dimensional example with nonlinear coefficients.

### Noisy van der Pol system

The noisy van del Pol system obeys the following stochastic differential equations:

\[d\left(\begin{array}{c}x_{1}\\ x_{2}\end{array}\right)=\left(\begin{array}{c}x_{2}\\ \epsilon x_{2}(1-x_{1}^{2})-x_{1}\end{array}\right)dt+\mathrm{diag}(\nu_{11}, \nu_{22})d\mathbf{W}(t), \tag{26}\]

where \(\epsilon\) and \(\nu_{11},\nu_{22}\) are parameters, and \(\mathrm{diag}(\cdot)\) means the diagonal matrix with the corresponding diagonal elements.

Figure 2: Estimated moments of the Ornstein-Uhlenbeck process by the neural networks constructed with the proposed method. (A) First order moment. (B) Second order moment. The dashed lines correspond to The analytical solutions.
The corresponding time-evolution operator \(\mathcal{L}\) and its adjoint one \(\mathcal{L}^{\dagger}\) are as follows:

\[\mathcal{L} =-\frac{\partial}{\partial x_{1}}x_{2}-\frac{\partial}{\partial x_{2 }}\left(\epsilon x_{2}\left(1-x_{1}^{2}\right)-x_{1}\right)+\frac{\partial^{2} }{\partial x_{1}^{2}}\frac{\nu_{11}^{2}}{2}+\frac{\partial^{2}}{\partial x_{2 }^{2}}\frac{\nu_{22}^{2}}{2}, \tag{27}\] \[\mathcal{L}^{\dagger} =x_{2}\frac{\partial}{\partial x_{1}}+\left(\epsilon x_{2}\left(1 -x_{1}^{2}\right)-x_{1}\right)\frac{\partial}{\partial x_{2}}+\frac{\nu_{11}^{ 2}}{2}\frac{\partial^{2}}{\partial x_{1}^{2}}+\frac{\nu_{22}^{2}}{2}\frac{ \partial^{2}}{\partial x_{2}^{2}}. \tag{28}\]

Then, the simultaneous ordinary differential equations for \(\{P_{i}^{(m)}(n_{1},n_{2},t)\}\) become

\[\frac{d}{dt}P_{i}^{(m)}(n_{1},n_{2},t) =(n_{1}+1)P_{i}^{(m)}(n_{1}+1,n_{2}-1,t)+\epsilon n_{2}P_{i}^{(m)} (n_{1},n_{2},t)\] \[\quad-\epsilon n_{2}P_{i}^{(m)}(n_{1}-2,n_{2},t)-(n_{2}+1)P_{i}^{ (m)}(n_{1}-1,n_{2}+1,t)\] \[\quad+\frac{\nu_{11}}{2}(n_{1}+2)(n_{1}+1)P_{i}^{(m)}(n_{1}+2,n_{ 2},t)\] \[\quad+\frac{\nu_{22}}{2}(n_{2}+2)(n_{2}+1)P_{i}^{(m)}(n_{1},n_{2}+ 2,t). \tag{29}\]

As in the Ornstein-Uhlenbeck case, a finite cutoff is employed; \(P_{i}^{(m)}(n_{1},n_{2},t)=0\) for \(n_{1}>N\) or \(n_{2}>N\).

In the demonstration, we set \(\epsilon=1\) and \(\nu_{11}=\nu_{22}=1\). For the noisy van der Pol systems, there is no analytical result. Hence, we use the estimated values obtained from the method in section 2 as the approximate true value.

The proposed method is the approach via the right side depicted in figure 1. For comparison, we perform the conventional machine learning approach with backpropagation. As explained in section 1, a naive approach, the left one in figure 1, requires samplings with the Monte Carlo method to estimate the expected values. However, it takes high computational costs. Hence, we here employ the approach via the middle one in figure 1; the expected values for various initial conditions are evaluated with the method in section 2. We make a data set with the pairs of an input coordinate and the target moment. Then, the neural network is learned from the data set. Of course, the Monte Carlo method with a considerably large data set gives the same learned network. The comparison approach requires various initial coordinates for the data set; we generated them from the uniform distribution with the range \([-4,4]^{2}\). The data size is \(250,000\); the size is chosen to give similar total estimation errors, as explained later. The learning procedure is performed with PyTorch with AdaDelta.

Again, note that the proposed method does not need any preparation for the data set. Instead, the coefficients \(\{P_{i}^{(m)}(\mathbf{n},t)\}\) are numerically evaluated and used.

The results are shown in figure 3. The upper ones correspond to the approximate true results obtained by the method in section 2. The middle ones are the results of the proposed method. The bottom ones correspond to those by the conventional learning approach based on backpropagation. As for the first-order moments of \(x_{1}\) and \(x_{2}\), it is difficult to see the differences among the three cases; the learning results are good enough. As for the second-order moment for \(x_{2}\), we see the differences from the true one. The proposed method gives a similar shape to the true one, while the color is a little thinning away from the origin. On the other hand, the result of backpropagation is worse even near the origin.
To examine this behavior in more detail, we evaluate the mean squared error between the learned results and approximate true ones for the distance from the origin. We calculate the outputs of the learned neural network values for various inputs with the \(100\times 100\) mesh grids for polar coordinates with the range \([0,4]\times[0,2\pi)\). The results are shown in figure 4. Note that the proposed method and the backpropagation method give similar total errors with distances from \(0\) to \(4\); the data size was chosen so that they would be. Note that we do not intend exactly to yield the same total errors because we here want to focus on the behavior of the errors with respect to the distance. In figure 4, we see a clear difference between the results; the proposed method gives more accurate estimations near the origin. This characteristic reflects the feature of the Taylor expansion in the proposed method.

## 5 Conclusion

We proposed a new learning method for neural networks, which target is the statistics of the stochastic differential equations. The proposed method directly compares coefficients obtained from the dual process and the weights in the neural networks.

Figure 3: Estimated moments for the noisy van der Pol system. (Upper) Approximately true results. (Middle) Results by the proposed method. (Bottom) Results by the conventional learning approach based on backpropagation.
There is no need for sampling procedures for the stochastic processes, and the proposed method gives a different framework from previous ones based on backpropagation. Although the idea is naive, the numerical demonstration shows remarkable features of the proposed approach; we have more accurate estimations near the origin, and the errors increase with distance. This feature stems from the Taylor-type basis expansion in the derivation of the method. Of course, it is easy to shift the origin of the state space using Ito's lemma, as commented in section 4.1. Hence, we can select another coordinate as the origin to estimate the moments around the selected coordinate more accurately. Problems such as overfitting are less likely to occur compared with the conventional learning methods from data sets because the proposed method does not require sampling of stochastic processes. These features will be hopeful in some applications in which one wants to embed the information of system equations into neural networks.

The present work is the first attempt utilizing direct comparisons, and the proposed method and the numerical examples would be enough as the demonstration. Of course, there are remaining works, as follows.

First, the estimations of the proposed method are less accurate except for inputs near the origin. Different cost functions with some focused coordinates could improve it.

Second, we should seek suitable numerical optimization methods to make the algorithm faster. It is interesting to investigate how the performance changes for neural

Figure 4: The estimated errors. The solid lines correspond to the results for the proposed method, and the dashed lines correspond to those for the backpropagation.
networks with other structures because we used a simple neural network with one hidden layer in the present work. Although the cases of multiple hidden layers are possible in principle, they require solving rather difficult nonlinear optimization problems. Hence, further development of numerical solvers would be necessary.

Third, practical applications would be investigated. For example, the following situation seems to be the case quite often: One knows the form of the time-evolution equations, but there are some ambiguities about the values of parameters. In such a situation, we can use the proposed method to construct an initial neural network based on the information of the equations; then, we apply additional learning steps to the initialized neural network. These procedures will reduce the size of data sets.

Fourth, we need further work to make the method practical in high-dimensional systems. We confirmed that the algorithm based on combinatorics [16, 25] can deal with at least four or five-dimensional systems rapidly. Furthermore, it is possible to evaluate the information of moments for systems with several dozen variables using the tensor-train format; see, for example, [26] for the tensor-train format. Of course, the optimization procedure will take longer computational time, and we need to develop further studies in the research community of optimizations. Although the practical applications for higher-dimensional cases are beyond the scope of the present paper, we hope that the present work will motivate further collaboration between statistical physics and machine learning.

The idea proposed in the present paper is the first step for future work. We hope that the present work will open up a new way to use the information of the time-evolution equation directly to learn neural networks.

This work was supported by JSPS KAKENHI Grant Number JP21K12045.
