# Weak Collocation Regression method: fast reveal hidden stochastic dynamics
from high-dimensional aggregate data

Liwei Lu  Zhijun Zeng  Yan Jiang  Yi Zhu  Pipi Hu Corresponding author,
<EMAIL_ADDRESS>
###### Abstract

Revealing hidden dynamics from the stochastic data is a challenging problem as
randomness takes part in the evolution of the data. The problem becomes
exceedingly complex when the trajectories of the stochastic data are absent in
many scenarios. Here we present an approach to effectively modeling the
dynamics of the stochastic data without trajectories based on the weak form of
the Fokker-Planck (FP) equation, which governs the evolution of the density
function in the Brownian process. Taking the collocations of Gaussian
functions as the test functions in the weak form of the FP equation, we
transfer the derivatives to the Gaussian functions and thus approximate the
weak form by the expectational sum of the data. With a dictionary
representation of the unknown terms, a linear system is built and then solved
by the regression, revealing the unknown dynamics of the data. Hence, we name
the method with the Weak Collocation Regression (WCR) method for its three key
components: weak form, collocation of Gaussian kernels, and regression. The
numerical experiments show that our method is flexible and fast, which reveals
the dynamics within seconds in multi-dimensional problems and can be easily
extended to high-dimensional data such as 20 dimensions. WCR can also
correctly identify the hidden dynamics of the complex tasks with variable-
dependent diffusion and coupled drift, and the performance is robust,
achieving high accuracy in the case with noise added.

Keywords— stochastic dynamics, the weak form, collocation of kernels, Fokker-
Planck equation, aggregate data, revealing hidden dynamics

##  1 Introduction

Nowadays, a large amount of data has been collected in different realms, and
revealing the hidden dynamics buried in the data is an essential topic in
science and engineering. Such dynamics can be described mathematically through
differential equations. In this regard, differential equations derived from
the so-called first principle successfully describe nature’s phenomena during
the past several centuries. The examples include the Navier-Stokes equation in
hydrodynamics for fluid dynamics [1], Schrödinger equation in quantum
mechanics for probability current [2], Black-Scholes in computational finance
for option pricing [3]. Often, the data we collect is a tremendous amount of
time series dataset, which accurately describes stochastic differential
equations. For instance, Brownian motion for molecular dynamics and geometric
Brownian motion for the financial market [4]. On the other hand, leveraging
explosive growth of available data, machine learning, especially deep
learning, has recently attained tremendous success in computer vision, natural
language processing, and many other topics in computer science [5].

Benefiting from these two, plenty of work spring up, and within them, one most
important focus is pouring new structures onto modeling the hidden dynamics
from data. The physical informed neural network (PINN) adds physical
constraints to the data by adding the residual of the differential equations
to the loss and making learning coefficients of the unknown terms of the
governing equations reliable [6]. Kutz et al. proposed a framework named
“SINDy” by combining regression and sparse identification to reveal nonlinear
dynamical systems [7]. The time-series data always contain many missing points
and flaws with the high noise, making analysis hard and tricky. Hu et al.
proposed using symbolic ODE (ordinary differential equations) to reveal hidden
dynamics from time series data [8] with the integration form, making the
learning process more stable. More related works about inferring differential
equations, see [9, 10, 11, 12, 13, 14, 15]. However, in a more general
scenario, the data come from the stochastic process where the determined-
dynamic system is not feasible because of the randomness. Revealing hidden
dynamics from the stochastic data calls for new structures suited to this
circumstance.

In the realm of dynamic systems, aggregate data refers to a data format in
which the trajectory of each individual modeled by the evolution of state is
not available, but rather a sample from the state distribution at a particular
time point is available [16]. For example, the data collected for single-cell
DNA sequence analysis, bird migration, and social gathering are aggregate
data, as the individual trajectories for a long time are impossible to follow
with only the collection of the states of different individuals obtained. In
contrast, trajectory data includes all the individual information along the
time. Thus, in other literature, trajectory and aggregate data are also called
“paired” and “unpaired” data, respectively [17].

For trajectory data, there exist many methods developed such as Hamiltonian
neural networks [18], Hidden Markov Model (HMM) [19], Kalman Filter (KF) [20],
Particle Filter (PF) [21] and related works [16, 22]. Revealing determined
hidden dynamics from data above can also be regarded as one special case of
trajectory data when stochastic term vanishes. Thus, plenty of works
investigate the trajectory data by revealing buried determined or under-
determined dynamics from real data.

However, few works are investigated for aggregate data because of the absence
of individual trajectories. By leveraging the Fokker-Planck equation, the
governing function of the probability density function of the Winner process
variable, researchers investigate the revealing hidden stochastic dynamics
from aggregate data. Zhou et al. propose a novel method using the weak form of
the Fokker Planck Equation (FPE), a partial differential equation, to describe
the density evolution of data in a sampled form, which is then combined with
the Wasserstein generative adversarial network (WGAN) in the training process
[23]. Chen et al. proposed a method that leverages the physical-informed
structure of Fokker-Planck equations and the approximation of the probability
density function and reveals the hidden dynamics from samples in several time
snapshots [24]. Liu et al. [17] integrated the distance measure such as
Wasserstein distance used in WGAN [25] and the forward numerical solution of
the parameterized stochastic differential equations as the Physics-Informed
Deep Generative Models to reveal the hidden dynamics of the data.

In this work, leveraging the weak form of the Fokker-Planck equation, we
proposed a Weak Collocation Regression (WCR) method to reveal the hidden
stochastic dynamics from nonequally-spaced temporal aggregate data. The
Fokker-Planck equation describes the time evolution of the probability density
function of the random variable in the Brownian motion. By the weak form of
the Fokker-Planck equation, one can reduce the aggregate data to the low
dimensional temporal sequence where the drift and diffusion term information
is buried. A Monte-Carlo sum is introduced to calculate the weak form with a
gaussian kernel as the test function. Temporal derivatives are reduced by
Linear Multistep Methods (LMMs), and then the linear system constructed by
collocation of the kernels gives an exact approximation of the stochastic
dynamics. The benefits are three folds.

  1. 1.

Comparable performance. Our method has a low computational cost, comparable
accuracy, and restrained error, and the dimensional curse can be lessened.
Numerical experiments show that the 1-dimensional problem can be easily
revealed within $0.02$ second on MacBook, while the computational time of the
three or 4-dimensional problem can be limited to seconds. We can significantly
reduce the error by changing the time interval and sample number. Directly
numerically solving the inverse problem of stochastic differential equations
may encounter the curse of dimensionality [26], where the tendency of
numerical techniques requires a high computational cost growing exponentially
with respect to the dimension of the variables. Benefits from the Monte-Carlo
summation and the random collocations methods, the curse of dimensionality is
lessened. Further, our method has natural potential for parallel computation,
which also reduces the cost.

  2. 2.

Robustness in different quality Data. Our method can handle a small amount of
non-equally noise data without trajectory information and even partially
missing. No trajectories are needed; the only thing needed is the summation of
all the points at each given time. Our method hence has natural permutation
invariance with respect to the observations of different individuals at the
same time snapshot. The measurement error can be reduced to a great extent
since the summation of the aggregate data is resistant to white noise.
Partially missing data is allowed because we only require the summation over
each time snapshot. Numerical experiments illustrate the stability of our
methods against random noise. Without trajectories, the method can also be
applied to different numbers of individuals at different time snapshots.

  3. 3.

Complex tasks. Our method can overcome the coupled drift term derived from
Sombrero potential, variable-dependent diffusion term, and high-contrast
problem and can be extended to a broader range. Numerical experiments show
that it obtained the same result if we expanded each drift and diffusion term
in a high-order polynomial. It shows that in dealing with the complex task
above, our method has a broader representation and shows superior performance
in revealing the hidden dynamics with a mild amount of data and high accuracy.

##  2 Methodology

In this paper, we consider the scenarios of the so-called aggregate data or
unpaired data in some literature, where the trajectory information of each
individual modeled by the evolution of the state is unavailable. Only the
collection of the $N_{t_{j}}$ samples
$\mathbb{X}_{t_{j}}=\\{\boldsymbol{x}^{i}_{t_{j}}\\}_{i=1}^{N_{t_{j}}}$ from
the distribution of state $X_{t}$ at certain time point $t_{j},j=1,2,\cdots,L$
is available, i.e., the snapshots.

For simplicity, we denote the set of the aggregate data as

$\mathbb{X}\overset{\Delta}{=}\\{\mathbb{X}_{t_{j}}\\}_{j=1}^{L}=\big{\\{}\\{\boldsymbol{x}^{i}_{t_{j}}\\}_{i=1}^{N_{t_{j}}}\big{\\}}_{j=1}^{L},$
(1)

where $t_{j}$ is the time index of the $j$-th snapshot in the total $L$
snapshots, $\boldsymbol{x}_{t_{j}}^{i}$ is the $i$-th sample of the variable
$X_{t}$ at time $t_{j}$. To model the stochastic dynamics of the data, the
general stochastic differential equations (SDEs) for Brownian motion are
considered

$dX_{t}=\boldsymbol{\mu}(X_{t},t)dt+\boldsymbol{\sigma}(X_{t},t)dW_{t},$ (2)

with drift term $\boldsymbol{\mu}(X_{t},t)\in\mathbb{R}^{d}$ and diffusion
term $\sigma(X_{t},t)\in\mathbb{R}^{d\times w}$ . Here, $X_{t}$ is the
$d$-dimensional random variable of the data, and $W_{t}$ is the
$w$-dimensional standard Brownian motion.

The density function $p(x,t)$ of the above random variable $X_{t}$ can be
described by the Fokker-Planck equation (FP), see [27], and we restate the
result in Lemma 1.

######  Lemma 1.

Suppose $\\{X_{t}\\}$ solves the SDEs (2), then the possibility density
function $p(x,t)$ of the random variable $X_{t}$ satisfies the following
d-dimensional Fokker-Planck equation by the Itô integral

$\partial_{t}p=-\nabla\cdot(\boldsymbol{\mu}p)+\sum_{i,j}^{d}\partial_{ij}(D_{ij}p),$
(3)

where $\boldsymbol{x}\in\Omega\subset\mathbb{R}^{d}$,
$t\in[0,T]\subset\mathbb{R}$, $p=p(\boldsymbol{x},t)\in\mathbb{R}$ is the
probability density function with $\int_{\Omega}p(\boldsymbol{x},t)dx=1$,
$\boldsymbol{\mu}=[\mu_{1}(\boldsymbol{x},t),\mu_{2}(\boldsymbol{x},t),\cdots,\mu_{d}(\boldsymbol{x},t)]^{T}$,
and the diffusion matrix $[D_{ij}]=[D_{ij}(\boldsymbol{x},t)]$ is given by

$D=\frac{1}{2}\boldsymbol{\sigma}\boldsymbol{\sigma}^{T}.$ (4)

The Fokker-Planck equation bridges the gap between stochastic dynamics and the
distributions of the data samples by sharing drift term and diffusion relation
in SDE (2) and FP equation (3). To reveal the hidden stochastic dynamics (2),
one can reveal the unknown terms in FP equation (3) by the samples following
the corresponding probability distribution.

However, there is a huge gap between the density function and the data
samples. Directly modeling the data by the Fokker-Planck equation usually
needs the temporal and spatial derivatives of the density function $p(t,x)$
with respect to $t$ and $x$, requiring a large number of samples for
constructing the smooth density function. It becomes even impossible to
construct the smooth density function when the dimension increases. To
overcome this difficulty, in this work, we propose a framework using the weak
form with collocation integral kernels instead of reconstructing density
function to fast reveal hidden stochastic dynamics with a mild amount of data.
We briefly introduce the methods here. With the weak form, the partial
derivatives moved to the explicit kernel function by integration by parts,
making the computation of the spatial derivatives much easier. Followed by the
Monte-Carlo summation, the integration over space can be given by the
summation over the samples. Moreover, the terms with temporal derivatives can
be well reduced by the linear multistep method (LMMs). Finally, the linear
system is built with a basis expansion of the unknown drift and diffusion, and
the sparse regression approximates the unknown terms. Figure 1 exhibits the
picture of the methodology.

The weak form, LMMs, and the regression model will be discussed in detail in
the following several subsections.

Figure 1: The diagram of the weak collocation regression method. The aggregate
data set $\mathbb{X}$ on panel (b) is the collection of $L$ snapshots of
samples at time $t_{1},t_{2},\ldots,t_{L}$ from one unknown stochastic
process. We model this process by the stochastic differential equations in
panel (a) with unknown drift $\boldsymbol{\mu}(X_{t},t)$ and diffusion
$\sigma(X_{t},t)$ terms. By sampling Gaussian kernels in panel (c), the weak
form in panel (d) gives the algebraic relation of the unknown terms and the
data set for each kernel. By the LMMs and the expansion of the unknown terms,
a linear system is built and further combined together to form a large system
over all of the collocation kernels. Finally, the STRidge method gives the
sparse regression of the drift and diffusion terms; hence, the hidden dynamics
are revealed.

###  2.1 Leverage the weak form of the Fokker-Planck equation

In PDE theories, weak solutions from the weak form of PDEs are vital in
analysis; see Evans’ book [28] for detail. The weak form is introduced to
solving PDEs with neural networks by the so-called WAN method [29] and then
extended to reveal the unknown parameters in PDEs [30]. Borrowing the ideas
from the Galerkin methods, we alter the FP equations’ weak form by using the
kernel function collocations as test functions. The generalities of test
functions are approximately fulfilled through the kernel collocation, unlike
the max update of the test function in WAN. Without maximum optimization for
the test function, minimum optimization for the unknown terms in the dynamics
achieves a more stable performance with high robustness.

The kernel function $\psi(x):\mathbb{R}^{d}\to\mathbb{R}$ is given as the test
function. For example, a typical choice can be the Gaussian function with the
form

$\phi(\boldsymbol{x},\boldsymbol{\rho},\boldsymbol{\gamma})=\Pi_{i=1}^{d}\frac{1}{\gamma_{i}\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{x_{i}-\rho_{i}}{\gamma_{i}}\right)^{2}},$
(5)

where $\boldsymbol{x}=(x_{1},\ldots,x_{d})^{T}\in\Omega$,
$\boldsymbol{\rho}=(\rho_{1},\ldots,\rho_{d})^{T}$ and
$\boldsymbol{\gamma}=(\gamma_{1},\ldots,\gamma_{d})^{T}$ are the expectation
and standard deviation of the Gaussian function.

Multiply $\psi(\boldsymbol{x})$ on both sides of equation (3), integrate by
parts, and one obtains

$\begin{split}\frac{d}{dt}\int_{\Omega}p(\boldsymbol{x},t)\psi(\boldsymbol{x})d\boldsymbol{x}=\int_{\Omega}p(\boldsymbol{x},t)\boldsymbol{\mu}(\boldsymbol{x},t)\cdot\nabla\psi(\boldsymbol{x})d\boldsymbol{x}+\int_{\Omega}p(\boldsymbol{x},t)\sum_{rs}^{d}D_{sr}\partial_{rs}\psi(\boldsymbol{x})d\boldsymbol{x},\end{split}$
(6)

where $\partial_{rs}=\frac{\partial^{2}}{\partial x_{r}\partial x_{s}}$ and
$r,s$ are indices of the dimension. Here we have used the fact that
$p(x,t)|_{\partial\Omega}\to 0$ holds for the large enough domain $\Omega$.

As $p(\boldsymbol{x},t)$ is the density of the variable $X_{t}$ at time $t$
with $\int p(\boldsymbol{x},t)d\boldsymbol{x}=1$, we can rewrite equation (6)
as the expectation form

$\frac{d}{dt}\mathbb{E}_{\boldsymbol{x}\sim
p(\boldsymbol{x},t)}[\psi(\boldsymbol{x})]=\mathbb{E}_{\boldsymbol{x}\sim
p(\boldsymbol{x},t)}[\boldsymbol{\mu}(\boldsymbol{x},t)\cdot\nabla\psi(\boldsymbol{x})]+\mathbb{E}_{\boldsymbol{x}\sim
p(\boldsymbol{x},t)}[\sum_{rs}^{d}D_{sr}\partial_{rs}\psi(\boldsymbol{x})],$
(7)

where $\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x},t)}$ is the expectation
over the probability density function $p(\boldsymbol{x},t)$.

###  2.2 Approximate the weak form using data

In the scenarios of real applications, the probability distribution
$p(\boldsymbol{x},t)$ is hard to obtain, but the data of many samplings from
the distribution is always feasible, thanks to modern technology and
instruments. By the law of large numbers, the expectation over the
distribution can be approximated by the summation over the samplings of the
variable related to the distribution, i.e.,

$\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\psi(\boldsymbol{x}_{t}^{i})\sim\mathbb{E}_{\boldsymbol{x}\sim
p(\boldsymbol{x},t)}[\psi(\boldsymbol{x})],$ (8)

where $N_{t}$ is the number of the samplings of the variable $X_{t}$ at time
$t$, and $\boldsymbol{x}^{i}_{t}$ is the $i$-th sampling.

For simplicity, we denote the data set over time as
$\mathbb{X}_{t}\overset{\Delta}{=}\\{\boldsymbol{x}^{i}_{t}\\}_{i=1}^{N_{t}}$.
Thus, with the data set of the samplings $\mathbb{X}_{t}$, the weak form of
the FP equation (3) yields

$\frac{d}{dt}\biggl{(}\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\psi(\boldsymbol{x}^{i}_{t})\biggr{)}=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\boldsymbol{\mu}(\boldsymbol{x}_{t}^{i},t)\cdot\nabla\psi(\boldsymbol{x}^{i}_{t})+\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}(\sum_{rs}^{d}D_{sr}\partial_{rs}\psi(\boldsymbol{x}^{i}_{t}))+\epsilon,$
(9)

where we have used the approximations of the expectations by the summations of
the samples and the discussion of the error term $\epsilon$ can be found in
the appendix. The drift and diffusion information is related to the data set
$\mathbb{X}_{t}$ in equation (9). The dimension of the data has been reduced
to one by the summation over the samplings with only the dependence on the
time $t$ left. This form reveals that our method is naturally suited for the
high-dimensional problem.

For simplicity, we denote

$\frac{d}{dt}y(\mathbb{X}_{t})=f(\mathbb{X}_{t},\boldsymbol{\mu},D),$ (10)

where

$y(\mathbb{X}_{t})=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\psi(\boldsymbol{x}^{i}_{t})$

and

$f(\mathbb{X}_{t},\boldsymbol{\mu},D)=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\boldsymbol{\mu}(\boldsymbol{x}_{t}^{i},t)\cdot\nabla\psi(\boldsymbol{x}^{i}_{t})+\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}(\sum_{rs}^{d}D_{sr}\partial_{rs}\psi(\boldsymbol{x}^{i}_{t}))$

are scalar functions varying over time $t$.

###  2.3 Approximate the temporal derivatives

We show a method like Tikhonov differentiation to directly calculate the
derivatives in the appendix because its results have a certain gap with the
method used in the text. Other than considering the approximation of the
derivatives, we note that the equation (10) is a 1-dimensional parameterized
ordinary differential equation where the solutions on the discrete times were
given with $y(\mathbb{X})$. The inverse problem of finding unknown terms in
the FP equation now reduces to the inverse problem of the 1-dimensional
ordinary differential equation (10). Many methods recently have been developed
to reveal the dynamics such as SINDy [31], PINN [6] and ODENet [8].
Considering the simplicity, in this work, we use the implicit form of Linear
Multistep Methods (LMMs) here to construct an explicit form of equation (10).

Linear Multistep Methods (LMMs) have been developed as popular numerical
schemes for the integration of the ordinary differential equations for the
known dynamic systems [32] with well-established mathematical theories [33].
Raissi et al. constructed multistep neural networks for the data-driven
discovery of nonlinear dynamical systems leveraging LMMs [34]. Recently Du et
al. applied LMMs in learning hidden dynamics from the data of a given state
with theoretical analysis of the stabilities and the convergence for the
inverse problem [35, 23]. In their work, an equally-spaced version of LMMs is
used. However, real data in the record usually contain many missing points and
flaws, making the equally-spaced version of LMM hard to apply. In this work,
with the help of variable step size Adams methods [36], we can deal with the
non-equally spaced temporal data for the discovery of the hidden dynamics.

For equally-spaced temporal data, Implicit Adams methods of trapezoidal rule,
Milne method, 2-step backward differentiation method (BDF2), and Adams-Moulton
methods are used in this work, some of them are listed below, and others can
be found in the appendix.

$\text{Trapezoidal rule:}\quad y_{n+1}-y_{n}=\frac{h}{2}(f_{n+1}+f_{n}),$ (11)
$\text{Milne method:}\quad
y_{n+1}-y_{n-1}=\frac{h}{3}(f_{n+1}+4f_{n}+f_{n-1}),$ (12)

where $y_{n}=y(t_{n})$, $f_{n}=f(t_{n},\theta)$ and $h\equiv
h_{n}=t_{n+1}-t_{n}$ for $n=1,2,\cdots,L-1$.

For more general temporal data without equally spaced time, the 2-step BDF-
formula of variable step size methods is as follows

$y_{n+1}-\frac{\left(1+\omega_{n}\right)^{2}}{1+2\omega_{n}}y_{n}+\frac{\omega_{n}^{2}}{1+2\omega_{n}}y_{n-1}=h_{n}\frac{1+\omega_{n}}{1+2\omega_{n}}f_{n+1},$
(13)

where $\omega_{n}=h_{n}/h_{n-1}$.

Recall that the Implicit Adams methods of the trapezoidal rule require only
two adjacent time points. Thus the variable step-size version reads

$y_{n+1}-y_{n}=\frac{h_{n}}{2}(f_{n+1}+f_{n}).$ (14)

Our numerical experiments show that the one-step rule of the LMMs performs
well when the samples on only three-time snapshots are available. In contrast,
the LMMs methods requiring two steps fail in such a scenario even if they have
a higher order of accuracy. Other numerical schemes, e.g., the schemes listed
in the appendix, have also been tested but with a worse performance than the
LMMs.

Either by Tikhonov differentiation or the LMMs, we derive the discrete form of
equation (10) as the following form

$\boldsymbol{\hat{y}}(\mathbb{X})=\boldsymbol{\hat{f}}(\mathbb{X},\boldsymbol{\mu},D),$
(15)

where $\boldsymbol{\hat{y}}(\mathbb{X})$ is a vector constructed by
$\boldsymbol{y}(\mathbb{X})=\left\\{y(\mathbb{X}_{t_{j}})\right\\}_{j=1}^{L}$,
and $\boldsymbol{\hat{f}}(\mathbb{X},\boldsymbol{\mu},D)$ is a vector with the
same size given by
$\boldsymbol{f}(\mathbb{X},\boldsymbol{\mu},D)=\left\\{f(\mathbb{X}_{t_{j}},\boldsymbol{\mu},D)\right\\}_{j=1}^{L}$.
Now the equation (15) gives us an algebraic equation, from which, by an apt
ansatz form for the drift and diffusion terms, we can reveal
$\boldsymbol{\mu}$ and $D$ with the data set $\mathbb{X}$. Although we list
Tikhonov differentiation in the appendix, it is much worse than LMMs in our
experiments. Thus, we omit the results leveraging the Tikhonov
differentiation.

###  2.4 Build the regression model

In this subsection, we investigate the sparse regression and collocations of
the kernels to solve $\boldsymbol{\mu}$ and $D$ from (15) with the data set
$\mathbb{X}$. The polynomial expressions are adopted for approximating the
drift and diffusion terms, with the basis coefficients to be determined. Let

$\Lambda=\\{1,x_{1},x_{2},\cdots,x_{d},x_{1}x_{1},x_{1}x_{2},\cdots,x_{d}^{p}\\}^{T}$
(16)

denote the $p$-th order complete polynomials with respect to the variable
$\boldsymbol{x}=(x_{1},\cdots,x_{d})^{T}$. The number of the terms of
$\Lambda$ is $b\overset{\Delta}{=}|\Lambda|=\begin{pmatrix}p+d\\\
p\end{pmatrix}$. Here, for simplicity, we suppose that $\boldsymbol{\mu}$ and
$D$ are independent of $t$. And the components $\mu_{i}$ of
$\boldsymbol{\mu}=[\mu_{i}]$ and $D_{ij}$ of $D=[D_{ij}]$ expand as

$\mu_{i}=\sum_{j=1}^{b}\zeta_{ij}^{\mu}\Lambda_{j},\quad\text{and
}D_{ij}=\sum_{k=1}^{b}\zeta_{ijk}^{D}\Lambda_{k},$ (17)

where $\Lambda_{j}$ is the $j$-th component in the basis set $\Lambda$.

Now, equation (10) has the following matrix form

$\frac{d}{dt}\boldsymbol{y}(\mathbb{X})=B(\mathbb{X})\boldsymbol{\zeta}$ (18)

where

$\boldsymbol{\zeta}=\\{\underbrace{\zeta_{11}^{\mu},\cdots,\zeta_{1b}^{\mu},\cdots,\zeta_{d1}^{\mu},\cdots,\zeta_{db}^{\mu}}_{\text{flatten
of
}\boldsymbol{\zeta}^{\mu}},\underbrace{\zeta_{111}^{D},\cdots,\zeta_{11b}^{D},\cdots,\zeta_{dd1}^{D},\cdots,\zeta_{ddb}^{D}}_{\text{flatten
of }\boldsymbol{\zeta}^{D}}\\}^{T},$ (19)

and the flatten scheme of $\boldsymbol{\zeta}^{\mu}$ and
$\boldsymbol{\zeta}^{D}$ is shown as

$\boldsymbol{\zeta}^{\mu}:\begin{pmatrix}\lx@xy@svg{\hbox{\raise
0.0pt\hbox{\kern
8.04486pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern
0.0pt\offinterlineskip\halign{\entry@#!@&&\entry<EMAIL_ADDRESS>0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{11}^{\mu}}$}}}}}}}{\hbox{\kern 16.7133pt\raise
0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{12}^{\mu}}$}}}}}}}{\hbox{\kern 41.35675pt\raise
0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\cdots}$}}}}}}}{\hbox{\kern 63.52519pt\raise
0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{1b}^{\mu}}$}}}}}}}{\hbox{\kern-7.9875pt\raise-24.76108pt\hbox{\hbox{\kern
0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{21}^{\mu}}$}}}}}}}{\hbox{\kern
16.7133pt\raise-24.76108pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern
3.0pt\raise 0.0pt\hbox{$\textstyle{\zeta_{22}^{\mu}}$}}}}}}}{\hbox{\kern
41.35675pt\raise-24.76108pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\cdots}$}}}}}}}{\hbox{\kern
63.52519pt\raise-24.76108pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{2b}^{\mu}}$}}}}}}}{\hbox{\kern-6.75pt\raise-48.94716pt\hbox{\hbox{\kern
0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\cdots}$}}}}}}}{\hbox{\kern
17.9508pt\raise-48.94716pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern
3.0pt\raise 0.0pt\hbox{$\textstyle{\cdots}$}}}}}}}{\hbox{\kern
41.35675pt\raise-48.94716pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\cdots}$}}}}}}}{\hbox{\kern
64.56436pt\raise-48.94716pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\cdots}$}}}}}}}{\hbox{\kern-8.04486pt\raise-73.13324pt\hbox{\hbox{\kern
0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{d1}^{\mu}}$}}}}}}}{\hbox{\kern
16.65594pt\raise-73.13324pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{d2}^{\mu}}$}}}}}}}{\hbox{\kern
41.35675pt\raise-73.13324pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\cdots}$}}}}}}}{\hbox{\kern
63.46783pt\raise-73.13324pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{db}^{\mu}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
63.52519pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
63.52519pt\raise-24.76108pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
63.46783pt\raise-73.13324pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
7.98752pt\raise-21.98418pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
6.75pt\raise-46.6598pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
8.04488pt\raise-70.4071pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces}}}}\ignorespaces\\\
\end{pmatrix},\quad\boldsymbol{\zeta}^{D}:\begin{pmatrix}\lx@xy@svg{\hbox{\raise
0.0pt\hbox{\kern
9.44487pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern
0.0pt\offinterlineskip\halign{\entry@#!@&&\entry<EMAIL_ADDRESS>0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern 22.34805pt\raise
0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern 41.30861pt\raise
0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{11b}^{D}}$}}}}}}}{\hbox{\kern 72.70485pt\raise
0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern 98.16762pt\raise
0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern 117.18552pt\raise
0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{d1b}^{D}}$}}}}}}}{\hbox{\kern-3.0pt\raise-23.26053pt\hbox{\hbox{\kern
0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
15.96056pt\raise-23.26053pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{112}^{D}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern
25.34805pt\raise-46.52106pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{$\scriptstyle{\hbox{\kern 6.0pt\vrule
height=3.0pt,depth=3.0pt,width=0.0pt}}$}}}}}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern
47.49777pt\raise-23.26053pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
72.70485pt\raise-23.26053pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
91.72275pt\raise-23.26053pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{d12}^{D}}$}}}}}}}{\hbox{\kern
123.43205pt\raise-23.26053pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern-9.3875pt\raise-46.52106pt\hbox{\hbox{\kern
0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{111}^{D}}$}}}}}}}{\hbox{\kern
22.34805pt\raise-46.52106pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
47.49777pt\raise-46.52106pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
66.25998pt\raise-46.52106pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{d11}^{D}}$}}}}}}}{\hbox{\kern
98.16762pt\raise-46.52106pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
123.43205pt\raise-46.52106pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern-3.0pt\raise-69.78159pt\hbox{\hbox{\kern
0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
22.34805pt\raise-69.78159pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
41.25124pt\raise-69.78159pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{1db}^{D}}$}}}}}}}{\hbox{\kern
72.70485pt\raise-69.78159pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
98.16762pt\raise-69.78159pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
117.12816pt\raise-69.78159pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{ddb}^{D}}$}}}}}}}{\hbox{\kern-3.0pt\raise-93.04211pt\hbox{\hbox{\kern
0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
15.90318pt\raise-93.04211pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{1d2}^{D}}$}}}}}}}{\hbox{\kern
47.49777pt\raise-93.04211pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
72.70485pt\raise-93.04211pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
91.66539pt\raise-93.04211pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{dd2}^{D}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern
75.70485pt\raise-93.04211pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{$\scriptstyle{\hbox{\kern 6.0pt\vrule
height=3.0pt,depth=3.0pt,width=0.0pt}}$}}}}}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern
123.43205pt\raise-93.04211pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern-9.44487pt\raise-116.30264pt\hbox{\hbox{\kern
0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{1d1}^{D}}$}}}}}}}{\hbox{\kern
22.34805pt\raise-116.30264pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
47.49777pt\raise-116.30264pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
66.20262pt\raise-116.30264pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\zeta_{dd1}^{D}}$}}}}}}}{\hbox{\kern
98.16762pt\raise-116.30264pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern
123.43205pt\raise-116.30264pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
45.12842pt\raise-4.94444pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
45.12842pt\raise-74.72603pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
121.03856pt\raise-4.94444pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
121.03856pt\raise-74.72603pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
2.1709pt\raise-66.78159pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
8.58067pt\raise-104.44487pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
77.88454pt\raise-66.78159pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
84.32025pt\raise-104.44487pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
47.49777pt\raise-26.02322pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces\ignorespaces\ignorespaces{}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern
123.43205pt\raise-26.0115pt\hbox{\hbox{\kern 0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces}}}}\ignorespaces\end{pmatrix}$

and $B(\mathbb{X}):=\begin{pmatrix}\boldsymbol{b}_{1}^{T}(\mathbb{X})\\\
\vdots\\\ \boldsymbol{b}_{L}^{T}(\mathbb{X})\end{pmatrix}$ is a known
coefficient matrix with size $L\times(db+d^{2}b)$,
$\boldsymbol{y}(\mathbb{X}):=\begin{pmatrix}y_{1}(\mathbb{X})\\\ \vdots\\\
y_{L}(\mathbb{X})\end{pmatrix}$ is a column vector with size $L$.

By applying linear multistep methods (LMMs) on the temporal derivative in
(18), one can obtain a linear system about the coefficient vector
$\boldsymbol{\zeta}$ as

$A(\mathbb{X})\boldsymbol{\zeta}=\boldsymbol{\hat{y}}(\mathbb{X}),$ (20)

where $A(\mathbb{X})$ and $\boldsymbol{\hat{y}}(\mathbb{X})$ are constructed
by $B(\mathbb{X})$ and $\boldsymbol{y}(\mathbb{X})$ in (18) because of the
linearity of the system, and the coefficient vector $\boldsymbol{\zeta}$ are
collected and vectorized from all the coefficients
$\left\\{\left\\{\zeta_{ij}^{\mu}\right\\}_{j=1}^{b}\right\\}_{i=1}^{d}$ and
$\left\\{\left\\{\zeta_{ijk}^{D}\right\\}_{k=1}^{b}\right\\}_{i,j=1}^{d}$ in
the expansions of $\boldsymbol{\mu}$ and $D$ shown in (19). Now, equation (20)
gives the relation between the data set $\mathbb{X}$ and the unknown
parameters in the polynomial expansions of the hidden dynamics.

######  Example 1.

To clarify, we give an example of how the linear system is obtained by
applying LMMs such as the trapezoidal rule (11) on the equation (18) with
temporal derivatives. By the trapezoidal rule, the $i$-th row vector
$\boldsymbol{a}^{T}_{i}(\mathbb{X})$ of the matrix $A(\mathbb{X})$ and the
$i$-th scalar $\hat{y}_{i}(\mathbb{X})$ of the vector
$\hat{\boldsymbol{y}}(\mathbb{X})$ are calculated as follows

$\boldsymbol{a}^{T}_{i}(\mathbb{X})=\frac{h}{2}(\boldsymbol{b}^{T}_{i}(\mathbb{X})+\boldsymbol{b}^{T}_{i+1}(\mathbb{X})),\quad\hat{y}_{i}(\mathbb{X})=y_{i+1}(\mathbb{X})-y_{i}(\mathbb{X}),\quad
i=1,2,\cdots,L-1$ (21)

Hence, the matrix $A(\mathbb{X})\in\mathbb{R}^{(L-1)\times(db+d^{2}b)}$ and
vector $\boldsymbol{\hat{y}}(\mathbb{X})\in\mathbb{R}^{L-1}$ in the linear
system can be assembled by

$A(\mathbb{X}):=\begin{pmatrix}\boldsymbol{a}^{T}_{1}(\mathbb{X})\\\ \vdots\\\
\boldsymbol{a}^{T}_{L-1}(\mathbb{X})\end{pmatrix}\quad\text{and
}\boldsymbol{\hat{y}}(\mathbb{X}):=\begin{pmatrix}\hat{y}_{1}(\mathbb{X})\\\
\vdots\\\ \hat{y}_{L-1}(\mathbb{X})\end{pmatrix}.$ (22)

Here, the rows of the matrix $A(\mathbb{X})$ or the elements of the vector
$\boldsymbol{\hat{y}}(\mathbb{X})$ are constructed by the rows of
$B(\mathbb{X})$ or the elements of $\boldsymbol{y}(\mathbb{X})$ where the
number of rows is reduced by one because of the trapezoidal rule.

A sparse regression applied to the linear system then reveals the data’s
hidden dynamics. However, before that, we would discuss the collocation
strategies to improve the robustness and accuracy. Recall that in this work,
the test function $\psi$ is considered as the Gaussian function
$\psi=\phi(\boldsymbol{x},\boldsymbol{\rho},\boldsymbol{\gamma})\overset{\Delta}{=}\Pi_{i=1}^{d}\frac{1}{\gamma_{i}\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{x_{i}-\rho_{i}}{\gamma_{i}}\right)^{2}}$
given by (5), one can easily obtain the the specific form of the linear system
(20) as

$A(\mathbb{X},\boldsymbol{\rho},\boldsymbol{\gamma})\boldsymbol{\zeta}=\boldsymbol{\hat{y}}(\mathbb{X},\boldsymbol{\rho},\boldsymbol{\gamma})$
(23)

by the replacement of
$\psi(x)=\phi(\boldsymbol{x},\boldsymbol{\rho},\boldsymbol{\gamma})$.

The collections of the test functions is taken as
$\mathbb{C}^{d}_{\boldsymbol{\rho},\boldsymbol{\gamma}}=\left\\{\phi(\cdot,\boldsymbol{\rho}_{m},\boldsymbol{\gamma}_{m})\right\\}_{m=1}^{M}$
where $\boldsymbol{\rho}_{m}\in\mathbb{R}^{d}$ and
$\boldsymbol{\gamma}_{m}\in\mathbb{R}^{d}$. Taking the $m$-th Gaussian
function $\phi(\cdot,\boldsymbol{\rho}_{m},\boldsymbol{\gamma}_{m})$ as the
test function, the linear system with respect to $\boldsymbol{\zeta}$ yields

$A_{m}\boldsymbol{\zeta}=\boldsymbol{\hat{y}}_{m},$ (24)

where $A_{m}=A(\mathbb{X},\boldsymbol{\rho}_{m},\boldsymbol{\gamma}_{m})$ and
$\boldsymbol{\hat{y}}_{m}=\boldsymbol{\hat{y}}(\mathbb{X},\boldsymbol{\rho}_{m},\boldsymbol{\gamma}_{m})$.
The values of the unknown coefficients $\boldsymbol{\zeta}$ can be revealed by
solving this linear system. However, integral with one test function
$\phi(\boldsymbol{x},\boldsymbol{\rho}_{i},\boldsymbol{\gamma}_{i})$ reveals
only part of the information from the given data. Hence, we build the linear
systems over the whole test function collection
$\mathbb{C}^{d}_{\boldsymbol{\rho},\boldsymbol{\gamma}}$. Namely, we solve the
linear equation

$\tilde{A}\boldsymbol{\zeta}=\boldsymbol{\tilde{y}},$ (25)

where

$\tilde{A}=\begin{pmatrix}A_{1}\\\ A_{2}\\\ \vdots\\\
A_{M}\end{pmatrix},\;\text{and}\;\boldsymbol{\tilde{y}}=\begin{pmatrix}\boldsymbol{\hat{y}}_{1}\\\
\boldsymbol{\hat{y}}_{2}\\\ \vdots\\\ \boldsymbol{\hat{y}}_{M}\end{pmatrix}$
(26)

are constructed by the matrix $\\{A_{m}\\}_{m=1}^{M}$ and the vector
$\\{\hat{\boldsymbol{y}}_{m}\\}_{m=1}^{M}$ over all of the test functions in
the collection $\mathbb{C}^{d}_{\boldsymbol{\rho},\boldsymbol{\gamma}}$.

To better approximate the entire test functional space, the collection
$\mathbb{C}^{d}_{\boldsymbol{\rho},\boldsymbol{\gamma}}$ needs to be large
enough, which means a huge amount of computational cost. To overcome this
difficulty, we borrow the ideas of the functional basis from the Galerkin
methods and the collocation points from the collocation methods. It is key to
generate the collection of the test functions
$\mathbb{C}^{d}_{\boldsymbol{\rho},\boldsymbol{\gamma}}$ taking into account
the efficiency and the accuracy. For simplicity, we randomly generate $M$
parameters of $\boldsymbol{\rho}_{m}$ in the range of the values of the data
$\mathbb{X}$. And the uniform standard derivation parameter
$\boldsymbol{\gamma}_{m}=\boldsymbol{1}$ is taken. Thus the test function
collection $\mathbb{C}^{d}_{\boldsymbol{\rho},\boldsymbol{\gamma}}$ is
obtained and the large linear system (25) is then built.

Followed by the linear regression and sparse identification algorithm
(STRidge) [7], the explicit form of the solved dynamics fits the data.
Specifically, to give the sparse results, the linear regression is applied but
with the hard threshold $\eta$ taken, i.e., the minor component of
$\boldsymbol{\zeta}$ less than $\eta$ would be set zero, and then the linear
regression continues until converges. Distributing the components of
$\boldsymbol{\zeta}$ to the drift and diffusion terms gives the explicit form
of the hidden stochastic dynamics.

The pseudo-code of the algorithm is exhibited in Algorithm 1. Figure 1 shows
the whole procedure of the methodology.

######  Remark 1.

Another choice of the ansatz for the drift and diffusion terms can be the
neural network because of its so-called universal approximation properties.
However, it would require nonlinear optimization, which enhances the
complexity of the algorithm, and left for further investigations beyond the
scope of this work. We focus on introducing the framework of the weak form to
reveal the hidden dynamics in this work.

Result: The explicit form of the governing stochastic equation of the data.

Input: Aggregate data set $\mathbb{X}$.

1 Sample $M$ expectations $\\{\boldsymbol{\rho}_{m}\\}_{m=1}^{M}$ in the
region $\Omega$ which contains all of the samples in the data set
$\mathbb{X}$;

2 Set the expectations $\\{\boldsymbol{\gamma}_{m}\\}_{m=1}^{M}$ as hyper-
parameters with default 1 in each dimension;

3 Generate the Gaussian function collections
$\mathbb{C}^{d}_{\boldsymbol{\rho},\boldsymbol{\gamma}}=\left\\{\phi(\cdot,\boldsymbol{\rho}_{m},\boldsymbol{\gamma}_{m})\right\\}_{m=1}^{M}$
using above expectations and derivations;

4

5for  _$m <=M$_ do

6 Compute and assemble the vector
$\boldsymbol{y}(\mathbb{X})=\begin{pmatrix}\vdots\\\
\frac{1}{N_{t_{l}}}\sum_{i=1}^{N_{t_{l}}}\phi(\boldsymbol{x}^{i}_{t_{l}},\boldsymbol{\rho}_{m},\boldsymbol{\gamma}_{m})\\\
\vdots\\\ \end{pmatrix}$ by the left term of the weak form (9) where the test
function
$\psi=\phi(\boldsymbol{x}^{i}_{t_{l}},\boldsymbol{\rho}_{m},\boldsymbol{\gamma}_{m})$
and $l$ is the index of the $l$-th time snapshots;

7

8 Compute and assemble the matrix $B(\mathbb{X})$ in (18) by the right term of
the weak form (9) over all of the time snapshots with the help of the
polynomial basis expansions of each entry of drift vector
$\boldsymbol{\mu}=[\mu_{i}]$ and diffusion matrix $D=[D_{ij}]$
$\mu_{i}=\sum_{j=1}^{b}\zeta_{ij}^{\mu}\Lambda_{j},\quad\text{and
}D_{ij}=\sum_{k=1}^{b}\zeta_{ijk}^{D}\Lambda_{k};$

9

10 Assemble the vector $\boldsymbol{\zeta}$ by
$\boldsymbol{\zeta}^{\mu}=[\zeta_{ij}^{\mu}]$ and
$\boldsymbol{\zeta}^{D}=[\zeta_{ijk}^{D}]$ using the flatten scheme in (19);

11

12 Compute the matrix $A_{m}$ and $\hat{\boldsymbol{y}}_{m}$ by the matrix
$B(\mathbb{X})$ and $\boldsymbol{y}(\mathbb{X})$ using the LMMs, e.g., the
trapezoidal rule (11);

13

14 end for

15Construct the large linear system
$\tilde{A}\boldsymbol{\zeta}=\boldsymbol{\tilde{y}}$ with
$\tilde{A}=(A_{1}^{T},\ldots,A_{M}^{T})^{T}$ and
$\boldsymbol{\tilde{y}}=(\boldsymbol{\hat{y}}_{1}^{T},\ldots,\boldsymbol{\hat{y}}_{M}^{T})^{T}$;

16 Compute $\boldsymbol{\zeta}$ by a sparse linear regression of
$\tilde{A}\boldsymbol{\zeta}=\boldsymbol{\tilde{y}}$;

Reconstruct the drift $\boldsymbol{\mu}$ and diffusion $D$ with the solved
coefficient $\boldsymbol{\zeta}$.

missingfnum@algorithm1Algorithm 1 Weak Collocation Regression method (WCR).

##  3 Numerical experiments

Data acquisition.  All the raw data used in this work for the experiments of
revealing the hidden dynamics are obtained by integrating the given SDEs

$d{X}_{t}=\boldsymbol{\mu}_{t}dt+\boldsymbol{\sigma}_{t}dW_{t},\quad t\geq 0$

from $t=0$ to $t=T$ with Euler-Maruyama scheme

$\tilde{{X}}_{(i+1)\delta t}=\tilde{{X}}_{i\delta
t}+\boldsymbol{\mu}_{t}\delta t+\boldsymbol{\sigma}_{t}\sqrt{\delta
t}\boldsymbol{\mathcal{N}}_{i}.$ (27)

Here $\delta t$ is the time step of the numerical scheme,
$\boldsymbol{\mathcal{N}}_{i}$ are i.i.d standard Gaussian random variables,
and the initial values are sampled from a given distribution such as a
Gaussian function. The same random seed is adopted among different
experiments. The experimental data is then sampled from these trajectories,
i.e., $N_{i}$ points are sampled at each time snapshot to remove the
trajectory information, and only $L$ time snapshots are chosen as the
experimental data, denoted as
$\mathbb{X}=\\{\mathbb{X}_{i}\\}_{i=1}^{L}=\\{\\{\boldsymbol{x}_{i}^{j}\\}_{j=1}^{N_{i}}\\}_{i=1}^{L}$.
The time snapshots can be non-equally spaced where the time interval $\Delta
t_{i}=t_{i+1}-t_{i}$ varies, and if the time interval is equal, we denote the
interval as $\Delta t$. The random noise is added to the raw data with the
noise level $\delta$ as
$\hat{\boldsymbol{x}}_{i}^{j}=\boldsymbol{x}_{i}^{j}+\delta\mathcal{U}_{i}^{j}\boldsymbol{x}_{i}^{j}$
where where $\mathcal{U}_{i}^{j}$ is a random variable.
Experimental setups.  In this work, we use the collocations of the Gaussian
functions as the test functions for the weak form. We sample these Gaussian
kernels by randomly sampling the expectations $\boldsymbol{\rho}_{m}$,
$m=1,2,\ldots,M$ using Latin Hypercube Sampling method [37] in the region of
the data, i.e., the hypercube containing all of the data. Moreover, the
standard derivation $\boldsymbol{\gamma}=\gamma I_{d}$, where $I_{d}$ is the
identity matrix, is chosen as a hyperparameter with default $\gamma=1$ in each
sampled Gaussian function. In all cases, our experiments show that all of the
nonzero coefficients are correctly identified, and all of the zero terms are
eliminated by the sparse regression. Hence in this work, we define the Maximum
Relative Error (MRE) of nonzero terms

$\text{MRE}=\max_{\theta_{i}\neq
0}\frac{|\hat{\theta}_{i}-\theta_{i}|}{|\theta_{i}|}$

as the criterion for the evaluation of the experimental results, where
$\theta_{i}$ represents the $i$-th parameter of the drift and diffusion terms
and $\hat{\theta}_{i}$ is the learned parameter from data. All the experiments
were done on the MacBook Pro 2021 with an M1 chip. We summarize the notations
used throughout the experiments in Table 1.

Variable | definition | Variable | definition
---|---|---|---
$L$ | Number of time snapshots | $\Delta t$ | Time interval of snapshots
$\gamma$ | Gaussian variance | $\boldsymbol{\rho}_{m}$ | Expectation of the $m$-th Gaussian kernel
$M$ | Gaussian sample number | $N_{i}$ | Sample number at $i$-th time snapshot
MRE | Maximum Relative Error | $\delta$ | Noise level
Table 1: Notations used throughout the experiments

###  3.1 Typical 1-dimensional problem

One-dimensional stochastic problem widely exists in the scientific and
engineering fields, such as population growth, asset price, and investments
[38]. To better illustrate the abilities of our WCR method for dealing with
complex tasks, in this section, we focus on revealing the hidden dynamics from
1d aggregate data avoiding the difficulties brought by the dimension of the
data. We consider the following 1d model

$dX_{t}=\mu(X_{t})dt+\sigma(X_{t})dW_{t},$ (28)

where drift $\mu(X_{t})$ and diffusion $\sigma(X_{t})$ terms reduce to scalar
functions. Here we mainly focus on the three cases: (I) Cubic polynomial
problem with only three snapshots; (II) Variable-dependent diffusion problem;
(III) Quintic polynomial problem with high contrast. Furthermore, these cases
give a direct illustration of the good performance of WCR on complex tasks.

(I) Cubic polynomial problem with three snapshots. In this case, the raw data
of the experiment are generated by the 1d cubic polynomial form

$dX_{t}=(X_{t}-X_{t}^{3})dt+dW_{t}$ (29)

in from $t=0$ to $t=1$ with initial values at $t=0$ sampled from a Gaussian
distribution $\mathcal{N}(0,0.1)$ using Euler-Maruyama scheme (27). We
obtained $10,000$ samples at each snapshot, and only three snapshots at (a)
$t=0.1,0.3,0.5$ and (b) $t=0.2,0.5,1$ are adopted as the experimental data. In
this case, we expand the drift term as the third-order basis expansion form

$\mu(x)=\lambda_{0}+\lambda_{1}x+\lambda_{2}x^{2}+\lambda_{3}x^{3},$

and the diffusion term is treated as a tunable parameter
$D(x)=\frac{1}{2}\sigma^{2}=D_{0}$. We have used the same experimental setup
for the same problem of Chen’s work [24].

Only three snapshots are available in the data set (a) and (b). Further, the
time intervals of the snapshots in (b) are not equal. Thus the approximation
of the temporal derivatives requires a reliable scheme to overcome the
difficulties. Here we use the variable step-size version of Implicit Adams
methods of trapezoidal rule (14), which can be applied to a small number of
time snapshots and non-equally spaced time series data. For the collocation
kernels, 20 Gaussian kernels are sampled by sampling the expectations of the
Gaussian function using the LHS method with the standard variance set as
$\gamma=0.85$.

The results are shown in Figure 2 and Table 2. We have compared the results
with Chen’s results [24] as state of the art (SOTA) to show that for the same
problem, the WCR method can achieve a comparable accuracy but with a much less
computational cost. All of the nonzero terms are correctly identified, and all
of the zero terms are eliminated by the sparse regression. The MRE of (a) and
(b) are less than $4.6\%$ and $0.42\%$, respectively. Most notably, the
experiments of applying WCR on the data set (a) and (b) are all completed
within 0.02s on a MacBook Pro 2021 with an M1 chip.

For a mild amount of the data, the WCR method still achieves remarkable
performances. We sampled 1000 points for each snapshot in case (b) above
without changing other setups to test our framework. In this case, the WCR
method still achieves a good result with the MRE less than $4\%$ within only
0.006s on the MacBook pro.

Figure 2: The results of 1d cubic polynomial problem compared with Chen’s work (SOTA). Reveal the unknown drift and diffusion terms with 10000 samples of $X_{t}$ at different time snapshots: (a) Observations at $t=0.1,0.3,0.5$; (b) Observations at $t=0.2,0.5,1$ where the samples are generated by the given SDE with drift term $\mu\\_\text{true}=x-x^{3}$ and diffusion term $\sigma=1$. The inference results are denoted by $\mu\\_{\text{chen}}$ (SOTA) and $\mu\\_{\text{wcr}}$ (ours). Parameter | $\lambda_{0}$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\sigma$
---|---|---|---|---|---
True parameters | 0 | 1 | 0 | $-1$ | 1
(a) Chen | $0.0051$ | $0.8422$ | $-0.0071$ | $-0.8994$ | $1.0347$
(a) WCR | $0$ | $\boldsymbol{1.0160}$ | $0$ | $\boldsymbol{-1.0457}$ | $\boldsymbol{1.0127}$
(b) Chen | $0.0225$ | $0.9638$ | $-0.0010$ | $-1.0035$ | $1.0138$
(b) WCR | $0$ | $\boldsymbol{0.9982}$ | $0$ | $\boldsymbol{-1.0019}$ | $\boldsymbol{1.0042}$
Table 2: The results of 1d cubic polynomial problem compared with Chen’s work
(SOTA). Reveal the unknown drift and diffusion terms with 10000 samples of
$X_{t}$ at different time snapshots: (a) Observations at $t=0.1,0.3,0.5$; (b)
Observations at $t=0.2,0.5,1$ where the samples are generated by the given SDE
with drift term $\mu\\_\text{true}=x-x^{3}$ and diffusion term $\sigma=1$. The
inference results are denoted by Chen (SOTA) and WCR (ours).

(II) Variable-dependent diffusion term problem. Much attention in the existing
literature is paid to the problem of constant diffusion in revealing the
hidden dynamics from data. However, a more general setting is the non-constant
diffusion, i.e., the variable-dependent diffusion. In this part, we will show
that the WCR method also works for the variable-dependent case. The variable-
dependent stochastic equation

$dX_{t}=(X_{t}-X_{t}^{3})dt+(1+X_{t})dW_{t}$

is considered as the true model to generate the raw data of the snapshots at
time $t=0,0.2,0.5,1$. Our task is to reveal the true model from the data by
the parameterized SDE form

$dX_{t}=(\lambda_{0}+\lambda_{1}X_{t}+\lambda_{2}X_{t}^{2}+\lambda_{3}X_{t}^{3})dt+(\sigma_{0}+\sigma_{1}X_{t})dW_{t},$
(30)

where $\lambda_{i},i=0,1,2,3$ and $\sigma_{j},j=0,1$ are the tunable
parameters. Note that the weak form of the Fokker-Planck is leveraged in our
method, where the drift term $\mu(x)$ and the $1\times 1$ diffusion matrix $D$
are revealed. By setting $D(x)=b_{0}+b_{1}x+b_{2}x^{2}$, WCR method gives the
approximation of the values of $\lambda_{i},i=0,1,2,3$ and $b_{j},j=0,1,2$.
And by the relation $D=\frac{1}{2}\sigma^{2}$, a nonlinear regression is
applied on the diffusion matrix $D$ and the parameters $\sigma_{i},i=0,1$ are
approximated, giving the explicit form of the governing equation (30). The
results of the parameters are listed in Table 3, and the WCR method still
works well in the variable-dependent problem.

Parameter | $\lambda_{0}$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\sigma_{0}$ | $\sigma_{1}$
---|---|---|---|---|---|---
True | 0 | 1 | 0 | -1 | 1 | 1
WCR | 0 | 0.9688 | 0 | -1.0264 | 0.9955 | 1.0326
Table 3: The results of 1d cubic polynomial problem when diffusion term is not
a constant. Reveal the unknown drift and diffusion terms with 10000 samples of
$X_{t}$ at $t=0,0.2,0.5,1$ where the samples are generated by the given SDE
with drift term $\mu\\_\text{true}=x-x^{3}$ and diffusion term $\sigma=1+x$.

(III) Quintic polynomial drift problem. Quintic polynomial drift has fifth-
order polynomial terms, raising the difficulty of revealing the true form of
the hidden stochastic dynamics. The problem becomes more subtle when the
coefficients of the polynomial terms are on different scales. In this
experiment, we would show that for the high order polynomial drift term with
high contrast coefficient, the WCR method can still reach a good accuracy in
modeling the data.

The raw data are generated by integrating the following SDE from $t=0$ to
$t=10$

$dX_{t}=\mu dt+dW_{t},$ (31)

where the drift term is the high order polynomial form with high contrast
coefficient as

$\mu\overset{\Delta}{=}-x(x-1)(x-2)(x-3)(x-4)=-24x+50x^{2}-35x^{3}+10x^{4}-x^{5}.$
(32)

In this experiment, WCR is applied to the aggregate data $\mathbb{X}$ to
reveal the hidden stochastic dynamics. Two kinds of the data with time
interval (a) $\Delta t=0.1$ and (b) $\Delta t=0.5$ with $N=5000$ points in
each time snapshot are considered. To model the data, we take the fifth order
of the polynomial expansion for the drift term with

$\mu(x)=\theta_{0}+\theta_{1}x+\theta_{2}x^{2}+\theta_{3}x^{3}+\theta_{4}x^{4}+\theta_{5}x^{5},$

where $\theta_{i},i=0,1,\ldots 5$ are the tunable parameters to be determined
by the data. The diffusion term is chosen as one tunable parameter $D_{0}$.

In the collocation of the kernels, 200 Gaussian functions are sampled by
randomly generating their expectations $\\{\rho_{i}\\}_{i=1}^{200}$ in the
region of the data with the default standard variance $\gamma=1$. For the
approximation of the temporal derivatives, the Milne method is applied. The
computation cost is within 2 seconds, and the results are summarized in Table
4. From the result, we can see that the WCR method achieves a good performance
in the problem with less than 3% Max Relative Error (MRE) within seconds on
the Macbook Pro.

Settings | drift coefficients $\\{\theta_{i}\\}_{i=0}^{5}$ | diffusion $D_{0}$ | MRE | Time (s)
---|---|---|---|---
True | 0 | -24 | 50 | -35 | 10 | -1 | 1 | - | -
(a) $\Delta t=0.1$ | 0 | -23.8 | 49.6 | -34.7 | 9.92 | -0.991 | 1.002 | 0.87% | 1.8
(b) $\Delta t=0.5$ | 0 | -24.9 | 51.6 | -35.9 | 10.2 | -1.018 | 1.029 | 3.82% | 0.4
Table 4: The results of one-dimensional high contrast problem with different
time snapshots. Reveal the unknown drift and diffusion terms with (a) $\Delta
t=0.1$; (b) $\Delta t=0.5$, where the samples are generated by the true SDE
with drift term $\mu(x)=-24x+50x^{2}-35x^{3}+10x^{4}-x^{5}$ and diffusion term
$\sigma=1$, gaussian functions $M=200$, samples number $N=5000$.

###  3.2 2-dimensional problem with coupled drift terms

Sombrero potential $V$ is a well-known potential form that symmetry breaking
is triggered in quantum mechanics. The gradient of the potential acts as the
force and gives the drift term of the stochastic process as $\mu=\nabla V$.
Figure 3 shows the sombrero potential $V$ and its gradient
$\mu_{1}=\partial_{x_{1}}V$ and $\mu_{2}=\partial_{x_{2}}V$ in the two
dimensional space. The drift term induced by the Sombrero potential has the
coupled terms in each dimension, i.e., $\mu_{1}$ and $\mu_{2}$ are not
independent. To illustrate the ability to reveal hidden stochastic dynamics
with coupled drift terms, we apply the WCR method to the aggregate data
generated by the 2-dimensional Brownian motion with the coupled drift term

$\boldsymbol{\mu}=-\nabla
V=\begin{pmatrix}10x_{1}-4x_{1}^{3}-4x_{1}x_{2}^{2}\\\
10x_{2}-4x_{2}x_{1}^{2}-4x_{2}^{3}\end{pmatrix},$ (33)

which is induced by the gradients of the Sombrero potential

$V=-5\|\boldsymbol{x}\|^{2}+\|\boldsymbol{x}\|^{4}=-5(x_{1}^{2}+x_{2}^{2})+(x_{1}^{2}+x_{2}^{2})^{2}.$
(34)

Figure 3:  The 2-dimensional sombrero potential (a) and its derivatives (b)
with respect to the spatial variables.

The true model decorated by the coupled drift (33) has the SDE form

$d\boldsymbol{x}=\begin{pmatrix}c_{11}x_{1}^{1}+c_{12}x_{1}^{3}+c_{13}x_{1}x_{2}^{2}\\\
c_{21}x_{2}^{2}+c_{22}x_{1}^{2}x_{2}+c_{23}x_{2}^{3}\end{pmatrix}dt+\begin{pmatrix}\sigma_{1}&0\\\
0&\sigma_{2}\end{pmatrix}dW_{t},$ (35)

where $c_{11}=c_{21}=10$, $c_{12}=c_{13}=c_{22}=c_{23}=-4$ corresponding to
the coupled drift term (33) and $\sigma_{1}=\sigma_{2}=1$ for simplicity.

The raw data are generated by integrating the SDE from $40,000$ initial
samplings from the Gaussian distributions at $t=0$ to $t=15$ using the Eular-
Maruyama method. Take $N=20,000$ samples out of the total $40,000$ points in
each snapshot and collect totally $L=151$ snapshots at $t=0,0.1,\ldots,15$, we
have the aggregate data set
$\mathbb{X}=\\{\mathbb{X}\\}_{i=1}^{L}=\\{\\{\boldsymbol{x}_{i}^{j}\\}_{j=1}^{N}\\}_{i=1}^{L}$
for the experiment. In this subsection, we apply the WCR method to the data to
reveal the true model and give a direct illustration of the performance of our
method.

To reveal the hidden dynamics of the 2-dimensional stochastic process, $200$
gaussian kernels are randomly sampled in the regime of the data, and the Milne
method (12) is used to approximate the temporal derivatives. The elements of
the unknown drift and coefficient term are approximated by the linear
combinations of the forth-order complete polynomial basis

$\Lambda=\\{1,x_{1},x_{2},x_{1}^{2},x_{1}x_{2},x_{2}^{2},x_{1}^{3},x_{1}^{2}x_{2},x_{1}x_{2}^{2},x_{2}^{3},,x_{1}^{4},x_{1}^{3}x_{2},x_{1}^{2}x_{2}^{2},x_{1}x_{2}^{3},x_{2}^{4}\\},$
(36)

with $|\Lambda|=15$. For simplicity we denote the basis vector
$\boldsymbol{\Lambda}=(1,x_{1},x_{2},\ldots,x_{2}^{4})^{T}$ which is the
vectorization of $\Lambda$. Note that the highest order of the drift term is
three in the true model, but we have used a higher order term to approximate
the unknown terms of both drift and diffusion terms. It follows that the
unknown drift and diffusion terms are approximated by the following form

$\boldsymbol{\mu}=\begin{pmatrix}\mu_{1}\\\ \mu_{2}\\\
\end{pmatrix}=\begin{pmatrix}\boldsymbol{\zeta}_{1}^{\mu}\cdot\boldsymbol{\Lambda}\\\
\boldsymbol{\zeta}_{2}^{\mu}\cdot\boldsymbol{\Lambda}\end{pmatrix},\quad
D=\begin{pmatrix}D_{11}&D_{12}\\\
D_{21}&D_{22}\end{pmatrix}=\begin{pmatrix}\boldsymbol{\zeta}_{11}^{D}\cdot\boldsymbol{\Lambda}&\boldsymbol{\zeta}_{12}^{D}\cdot\boldsymbol{\Lambda}\\\
\boldsymbol{\zeta}_{21}^{D}\cdot\boldsymbol{\Lambda}&\boldsymbol{\zeta}_{22}^{D}\cdot\boldsymbol{\Lambda}\end{pmatrix},$
(37)

where
$\boldsymbol{\zeta}^{\mu}_{i}\in\mathbb{R}^{d},\boldsymbol{\zeta}^{D}_{ij}\in\mathbb{R}^{d}$
are coefficient vectors with dimension $d=|\Lambda|=15$. Thus for each drift
term and each element of the diffusion matrix, we have 15 coefficients, making
a total of 90 parameters to be revealed. By a sparse linear regression,
instead of the 90 coefficients for a fourth-order complete polynomial basis,
only eight coefficients corresponding to those in (35) are listed for
simplicity, as the others are identified correctly as zero. All nonzero terms
of the true model are not missed, and no redundant coefficients, which should
be zeros, are not superfluous by the WCR method. The results are shown in
Table 5. The Max Relative Error (MRE) is about $2.42\%$ of the coefficients.
The experiment is done on the Macbook pro within $21$ seconds. The results in
Table 5 illustrate the abilities of the WCR method for revealing hidden
dynamics with coupled terms under high-order expansion.

Parameter | $c_{11}$ | $c_{12}$ | $c_{13}$ | $c_{21}$ | $c_{22}$ | $c_{23}$ | $\sigma_{1}$ | $\sigma_{2}$
---|---|---|---|---|---|---|---|---
True parameters | 10 | -4 | -4 | 10 | -4 | -4 | $1$ | 1
WCR | 9.7823 | -3.9032 | -3.9364 | 9.8495 | -3.9441 | -3.9409 | $0.9881$ | 1.0076
Table 5: The results of the 2-dimensional problem with Sombrebro potential.
Reveal the unknown drift and diffusion terms with the expansion of the fourth-
order complete polynomial basis. No redundant coefficients have been learned,
with only 8 out of the 90 coefficients nonzero corresponding to the true
model. The Max Relative Error of the nonzero coefficients is 2.42%.

###  3.3 Multi-dimensional problem

Many methods depending on the integration over space, always get stuck when
the dimension of the data increases because of the exponential increase of the
computational cost. Even under three or four dimensions, the integration is
computationally expensive, despite the incapacity for higher dimensions such
as 10d or 20d. In contrast, the WCR method takes a more subtle strategy to
avoid the direct computation of the integral. Thanks to the weak form, the
spatial derivatives of the probability function have been transferred to the
test function, making the integral simply computed by the summation over the
samples avoiding the curse of dimensionality. Hence WCR method can be
naturally extended to the higher dimension other than one or two dimensions.

In this subsection, we take the same setup consistent with the 3d and 4d
problems in [24] and show that the WCR method can achieve high accuracy within
seconds. The true model is a d-dimensional extension of the SDE (29) with the
form

$dX^{i}_{t}=\left(X^{i}_{t}-(X^{i}_{t})^{3}\right)dt+dW^{i}_{t},\quad
i=1,2,\cdots,d,$ (38)

where the drift term is $\mu_{i}=x_{i}-x_{i}^{3}$ in the $i$-th dimension. The
raw data are generated by the true model using the Euler-Mayaruma method. To
get the aggregate data set, $100,000$ samples each snapshot are taken and the
snapshots at (i) $t=0.1,0.3,0.5,0.7,1$ (ii) $t=0.1,0.2,0.3,0.5,0.7,0.9,1$ are
collected as the experimental data.

In the experiment, the unknown drift term in each dimension is approximated by

$\hat{\mu}_{i}=\theta_{0}+\theta_{1}x_{i}+\theta_{2}x_{i}^{2}+\theta_{3}x_{i}^{3},$
(39)

and the diffusion term approximated by a tunable parameter $D_{i}$,
$i=1,2,\ldots,d$. $100$ Gaussian kernels are sampled by the expectations using
the LHS method with the default standard derivation where $\gamma=1$. We list
the Max Relative errors and the computational time in Table 6, and the results
show that our WCR method achieves high accuracy in the 3d and 4d cases within
seconds on the MacBook Pro.

| 3D-(i) | 3D-(ii) | 4D-(i) | 4D-(ii)
---|---|---|---|---
MRE | 3.39% | 1.86% | 7.04% | 3.18%
Time(s) | 3.6 | 4.7 | 6.1 | 7.9
Table 6: The Max Relative Error of the learned coefficients and the
computational time on the MacBook Pro of revealing the unknown 3d, 4d cubic
polynomial problems. The aggregate data $\mathbb{X}$ is composed of the time
snapshots at (i) $t=0.1,0.3,0.5,0.7,0.9$ (ii) $t=0.1,0.2,0.3,0.5,0.7,0.9,1$
with $100,000$ points in each snapshot. The unknown dynamics is approximated
by the drift
$\hat{\mu}_{i}=\theta_{0}+\theta_{1}x_{i}+\theta_{2}x_{i}^{2}+\theta_{3}x_{i}^{3}$
and diffusion $D_{i}$ in each dimension with tunable parameters $\theta_{i}$
and $D_{i}$, $i=1,2,\ldots,d$. Moreover, $100$ gaussian kernels are sampled to
give the composed linear system.

Figure 4: The learned drift terms of the 3d and 4d problems. Reveal the
unknown drift and diffusion terms with $100,000$ samples each snapshot at time
$t=0.1,0.3,0.5,0.7,1$ of (a) 3-dimensional problem; (b) 4-dimensional problem,
where the samples are generated by the given SDE with drift term
$\mu\\_\text{true}=x-x^{3}$ and diffusion term $\sigma=1$ in per dimension.
The inference results are denoted by $\mu\\_\text{true}$ and
$\mu_{i}\\_\text{wcr}$ for the $i$-th dimension.

###  3.4 High-dimensional problem

For three or four-dimensional data, the WCR method shows high accuracy in
revealing the hidden stochastic dynamics. The higher dimensional problem,
especially over ten-dimensional, brings more challenges to the modeling task
because of the so-called curse of dimensionality. In this subsection, we will
show that the WCR method can be applied to high-dimensional problems thanks to
the Monte Carlo approximation of the weak form. We consider the ten and
20-dimensional aggregate data generated by the true model

$dX_{t}^{i}=\left(X_{t}^{i}-(X_{t}^{i})^{3}\right)dt+dW_{t}^{i},\qquad
i=1,2,\cdots,d,$ (40)

with $d=10,20$ using Euler-Maruyama method. Several time snapshots of the data
at time $t=0,0.1,0.2,0.3,\cdots,1$ are used as the experiment data.

In the experiment, the fourth order basis is used to approximate the drift
term in each dimension as
$\mu_{i}=\theta_{0}+\theta_{1}x_{i}+\theta_{2}x_{2}^{2}+\theta_{3}x_{3}^{3}$
and the diffusion term is approximated by a tunable parameter $\sigma_{i}$,
making totally $5d$ tunable parameters to be learned. The Milne method is used
for the temporal derivatives. The different number of data samples in each
snapshot and the Gaussian kernels are considered in the high-dimensional
problem.

In the 10-dimensional problem, we consider the four cases for samples and
kernels (a) $10,000$ samples of $X_{t}$ using $1000$ Gaussian kernels; (b)
$100,000$ samples of $X_{t}$ using $1000$ Gaussian kernels; (c)
$10,000$samples of $X_{t}$ using $10,000$ Gaussian kernels; (d) $100,000$
samples of $X_{t}$ using $10,000$ Gaussian kernels. The four cases are the
combinations of the $10,000$ and $100,000$ samples with $1000$ and $10,000$
Gaussian kernels. Figure 5 gives the display of the learned drift terms
compared with the true model in each case. From the Figure, it is easy to
check that the learned results get better when the number of the sample points
and Gaussian kernels increases. Moreover, the Max Relative Error for the best
result in case (d) achieves less than $7.8\%$ within 20 minutes.

Figure 5: The results of 10d problems with the different number of samples per
time snapshot and gaussian kernels. Reveal the unknown drift and diffusion
terms with (a) $10,000$ samples of $X_{t}$ using 1000 gaussian kernels; (b)
$100,000$ samples of $X_{t}$ using 1000 gaussian kernels; (c) $10,000$ samples
of $X_{t}$ using $10,000$ gaussian kernels; (d) $100,000$ samples of $X_{t}$
using $10,000$ gaussian kernels at time snapshots
$t=0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1$ of 10-dimensional problem, where
the samples are generated by the given SDE with drift term
$\mu\\_\text{true}=x-x^{3}$ and diffusion term $\sigma=1$ in per dimension.
The inference results are denoted by $\mu_{i}$ for the learned drift in the
$i$-th dimension.

In the 20-dimensional case, we take $10,000$ gaussian kernels and $100,000$
samples of $X_{t}$. The experiment was completed within 47 minutes on the
MacBook Pro, and the result is shown in Figure 6

Figure 6: The results of 20-dimensional problems with $100,000$ samples of
$X_{t}$ using $10,000$ gaussian kernels at time
$t=0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1$, where the samples are generated
by the given SDE with drift term $\mu\\_\text{true}=x_{i}-x_{i}^{3}$ and
diffusion term $\sigma=1$ in per dimension. The inference results are denoted
as $\mu_{i}$ for the $i$-th dimensional drift.

###  3.5 Performance on data of different qualities

In real applications, the data obtained by various means usually contain many
missing points and even flaws. The data may also be noisy, with sound signals
deeply buried. These facts make analyzing and extracting valuable models from
the data complicated and tricky. We have shown that WCR performs well with
only several snapshots with non-equally spaced data in the one-dimensional
case. In this subsection, we further investigate the performance of WCR on the
data with different qualities, such as short and long time intervals,
different numbers of samples in the snapshots, and noises.

The true model in this experiment has the following form

$\begin{split}dx_{1}&=-0.5x_{1}dt+dW_{t}\\\ dx_{2}&=-0.7x_{2}dt+dW_{t}\\\
dx_{3}&=-x_{3}dt+dW_{t},\end{split}$ (41)

where the three dimensions are not coupled for simplicity and mainly focus on
the data qualities. The raw data are generated by integrating the true model
from $t=0$ to $t=10$, and different kinds of snapshots are chosen as the
aggregate data. Noise is added for the study of the noise effect.

To reveal the hidden dynamics, the unknown drift term is expanded by the first
order complete polynomials as

$\mbox{drift}=\begin{pmatrix}\theta_{00}+\theta_{01}x_{1}+\theta_{02}x_{2}+\theta_{03}x_{3}\\\
\theta_{10}+\theta_{11}x_{1}+\theta_{12}x_{2}+\theta_{13}x_{3}\\\
\theta_{20}+\theta_{11}x_{1}+\theta_{22}x_{2}+\theta_{23}x_{3}\end{pmatrix},$
(42)

and the diffusion matrix is approximated by a diagonal matrix
$D=\begin{pmatrix}\sigma_{1}&&\\\ &\sigma_{2}&\\\ &&\sigma_{3}\\\
\end{pmatrix}$, where $[\theta_{ij}]$ and $[\sigma_{i}]$ are tunable
parameters to be revealed.

Time interval and sample number. We investigate in this part how the time
interval and the sample number of the data affect the accuracy of the results
in our WCR method. We consider four cases of the time interval $\Delta
t_{1}=0.1$, $\Delta t_{1}=0.2$, $\Delta t_{1}=0.5$ and $\Delta t_{1}=1$ of the
snapshots in the total interval $[0,10]$. To be concise, the total number of
the snapshots $L=\frac{T}{\Delta t}+1$, where $T=10$, and $\Delta t$ is one of
the time intervals. For example, we have $L=11$ snapshots in the case of
$\Delta t_{4}=1$, with each snapshot being equally spaced with $\Delta
t_{4}=1$. For the number of the samples in each snapshot, we also consider
four cases with $N_{1}=1000$, $N_{2}=2000$, $N_{3}=5000$, and $N_{4}=10,000$.
Thus we have totally tested $4\times 4=16$ cases to investigate the
performance of WCR on the time interval and sample number of the data.

In all experiments, other setups are kept the same where 100 Gaussian
functions with standard variance $\gamma=1$ are sampled as the test functions.
Moreover, the trapezoidal rule is used for the temporal derivatives if $\Delta
t=0.1,0.2$ and the Milne method is used for the big time interval $\Delta
t=0.5,1$. The Milne method with higher order accuracy is used in the big
intervals to achieve better accuracy.

We summarize all the results in Figure 7. With only $N=1000$ samples and mild
time interval $\Delta t=1$, the maximum relative error can reach a good
accuracy below $4\%$. Further, with the increasing number of samples in each
snapshot and the decreasing time interval from $\Delta t_{4}=1$ to
$\Delta_{1}=0.1$, from $N_{1}=1000$ to $N_{2}=10,000$, the results generally
get better. Moreover, all of the cases result in a low Max Relative Error,
giving proof of the robustness of our method in relatively poor data quality.

Figure 7: The heat map of the errors of the drift coefficients and diffusion
for the 3-dimensional stochastic process in $T=[0,10]$ with different
trajectory samples $N_{1}=1000$, $N_{2}=2000$, $N_{3}=5000$, $N_{4}=10,000$
and different time interval $\Delta t_{1}=0.1$, $\Delta t_{2}=0.2$, $\Delta
t_{3}=0.5$, and $\Delta t_{4}=1$. The true governing stochastic equation reads
$dx_{1}=-0.5x_{1}dt+dW_{1}$, $dx_{2}=-0.7x_{2}dt+dW_{2}$, and
$dx_{3}=-x_{3}dt+dW_{3}$, where the coefficients of $\mu_{1}$, $\mu_{2}$, and
$\mu_{3}$ are $-0.5$, $-0.7$ and $-1$, and the diffusions
$\sigma_{1}=\sigma_{2}=\sigma_{3}=1$. (a) shows the max relative error of the
coefficients for $N_{i}$ and $\Delta t_{j}$; (b) shows the relative error of
per coefficient for $N_{i}$ and $\Delta t_{j}$ with $i,j=1,2,3,4$.

Noise. In this part, we consider the performance of our WCR method on noisy
data. The raw data are obtained by the same procedure as the above, and the
random noise is added with the noise level $\delta$

$\hat{\boldsymbol{x}}_{i}^{j}=\boldsymbol{x}_{i}^{j}+\delta\mathcal{U}_{i}^{j}\boldsymbol{x}_{i}^{j},$

where $\mathcal{U}_{i}^{j}$ is an uniform random variable in $[0,1]$. Here the
data with $10,000$ samples in each snapshot are used, and two kinds of time
intervals are chosen to show the robustness of our methods on the noise.

The results are listed in Table 7, and the Max Relative Error reaches around
5% with the white noise level $\delta=10\%$ for both two cases. Further, when
we dive into the detail of the accuracy for drift and diffusion terms,
respectively, adding noise to the data did not change the accuracy of the
drift term much, unlike the diffusion term.

| $\Delta t=0.1,\delta=0\%$ | $\Delta t=0.1,\delta=10\%$ | $\Delta t=0.5,\delta=0\%$ | $\Delta t=0.5,\delta=10\%$
---|---|---|---|---
MRE | 1.87% | 5.76% | 1.46% | 6.01%
MRE in drift | 1.87% | 1.75% | 1.46% | 1.04%
MRE in diffussion | 0.72% | 5.76% | 0.77% | 6.01%
Table 7: The results of a three-dimensional linear problem with different time
snapshots and noise. Reveal the unknown drift and diffusion terms with $\Delta
t=0.1,0.5$ and $\delta=0\%,10\%$, where $10,000$ samples each snapshot are
generated by the given SDE with drift term $\mu(x,y,z)=[-0.5x,-0.7y,-z]^{T}$
and diffusion term $\sigma=I_{3\times 3}$.

##  4 Conclusion

In this work, leveraging the weak form of the Fokker-Planck equation and the
collocations of the Gaussian kernels, we proposed a framework called the Weak
Collocation Regression method (WCR) to fast reveal the hidden stochastic
dynamics from high-dimensional aggregate data. The lack of trajectory
information makes aggregate data (unpaired data) more difficult to reveal the
hidden dynamics. However, the data distribution follows the Fokker-Planck
equation under some assumptions, such as the Brownian motion. By transferring
the spatial derivatives to the test function in the weak form of the Fokker-
Planck equation, we have an integral form of the density function with
derivatives of the test functions. Thus the integral form can be easily
approximated by summing the derivatives of the test functions over the
samples. Using collocations of the Gaussian functions as the test functions
and the dictionary representation of the unknown terms, we build an extensive
linear system. Linear regression and sparse identification lead to unknown
terms, thus revealing the hidden dynamics.

Our experiments show that WCR is numerically efficient, taking only $0.02s$ on
the Macbook Pro for the $1$-dimensional problem with comparable accuracy and
seconds for multi-dimensional problems. With the LMMs of variable step size,
WCR achieves good accuracy in the non-equally spaced time data. Approximating
the integral form using the summations over samples makes the WCR method
naturally relieve the curse of dimensionality. Furthermore, the extra
computational cost is required on the collocations of the kernels but does not
exponentially increase with dimension. WCR still works well when the dimension
increases to $20$ within an hour on MacBook Pro. The variable-dependent
diffusion problem and coupled drift terms usually bring more difficulties, but
the WCR method exhibits high accuracy in these complex problems. Furthermore,
for the data with noise, the WCR method also obtains a stable performance when
adding the white noise. With different time intervals and samples, WCR shows
its robustness without losing too much accuracy when the data quality worsens.

Despite the advantages mentioned above, the Weak Collocation Regression method
still has some limitations for investigation in future work. Our approach
slows down for the ultimately high-dimensional problems such as 100 dimensions
since more Gaussian kernel collocations are needed, requiring high
computational costs. More effective sampling schemes, such as active learning
for the collocation of the Gaussian kernels, are necessary to reduce the
computational cost. Neural networks might be another choice for the test
function with a min-max optimization for reducing the cost. Further, the drift
and diffusion terms may be approximated by neural networks to give a good
approximation in the high dimension. Further discussion about the convergence
results and a general form besides the Brownian motion is needed but left for
future work.

## Acknowledgements

This work was supported by the National Key R&D Program of China (Grant No.
2021YFA0719200) and the National Science Foundation of China (Grant No.
12071244). The authors would like to thank the helpful discussions from Dr.
Liu Hong.

## References

  * [1] Roger Temam.  Navier-Stokes equations: theory and numerical analysis, volume 343\.  American Mathematical Soc., 2001.
  * [2] Leslie E Ballentine.  Quantum mechanics: a modern development.  World Scientific Publishing Company, 2014.
  * [3] Peter Klein.  Pricing black-scholes options with correlated credit risk.  Journal of Banking & Finance, 20(7):1211–1229, 1996.
  * [4] Daniel Dufresne.  The integral of geometric brownian motion.  Advances in Applied Probability, 33(1):223–241, 2001.
  * [5] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  MIT press, 2016.
  * [6] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [7] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937
  * [8] Pipi Hu, Wuyue Yang, Yi Zhu, and Liu Hong.  Revealing hidden dynamics from time-series data by odenet.  Journal of Computational Physics, 461:111203, 2022.
  * [9] Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong.  Pde-net: Learning pdes from data.  In International Conference on Machine Learning, pages 3208–3216, 2018.
  * [10] Zichao Long, Yiping Lu, and Bin Dong.  Pde-net 2.0: Learning pdes from data with a numeric-symbolic hybrid deep network.  Journal of Computational Physics, 399:108925, 2019.
  * [11] Gert-Jan Both, Subham Choudhury, Pierre Sens, and Remy Kusters.  Deepmod: Deep learning for model discovery in noisy data.  Journal of Computational Physics, 428:109985, 2021.
  * [12] Zhao Chen, Yang Liu, and Hao Sun.  Physics-informed learning of governing equations from scarce data.  Nature communications, 12(1):1–13, 2021.
  * [13] Guochang Lin, Pipi Hu, Fukai Chen, Xiang Chen, Junqing Chen, Jun Wang, and Zuoqiang Shi.  Binet: Learning to solve partial differential equations with boundary integral networks.  arXiv preprint arXiv:2110.00352, 2021.
  * [14] Qin Lou, Xuhui Meng, and George Em Karniadakis.  Physics-informed neural networks for solving forward and inverse flow problems via the boltzmann-bgk formulation.  Journal of Computational Physics, 447:110676, 2021.
  * [15] Lu Lu, Pengzhan Jin, and George Em Karniadakis.  Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.  arXiv preprint arXiv:1910.03193, 2019.
  * [16] Shaojun Ma, Shu Liu, Hongyuan Zha, and Haomin Zhou.  Learning stochastic behaviour from aggregate data.  In International Conference on Machine Learning, pages 7258–7267. PMLR, 2021.
  * [17] Liu Yang, Constantinos Daskalakis, and George E Karniadakis.  Generative ensemble regression: Learning particle dynamics from observations of ensembles with physics-informed deep generative models.  SIAM Journal on Scientific Computing, 44(1):B80–B99, 2022.
  * [18] Samuel Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian neural networks.  Advances in Neural Information Processing Systems, 32, 2019.
  * [19] Daniel Alshamaa, Aly Chkeir, Farah Mourad-Chehade, and Paul Honeine.  A hidden markov model for indoor trajectory tracking of elderly people.  In 2019 IEEE Sensors Applications Symposium (SAS), pages 1–6. IEEE, 2019.
  * [20] Maryam Baradaran Khalkhali, Abedin Vahedian, and Hadi Sadoghi Yazdi.  Vehicle tracking with kalman filter using online situation assessment.  Robotics and Autonomous Systems, 131:103596, 2020.
  * [21] Nuno Pessanha Santos, Victor Lobo, and Alexandre Bernardino.  Unmanned aerial vehicle tracking using a particle filter based approach.  In 2019 IEEE Underwater Technology (UT), pages 1–10. IEEE, 2019\.
  * [22] Yongkun Fang, Chao Wang, Wen Yao, Xijun Zhao, Huijing Zhao, and Hongbin Zha.  On-road vehicle tracking using part-based particle filter.  IEEE transactions on intelligent transportation systems, 20(12):4538–4552, 2019.
  * [23] Qiang Du, Yiqi Gu, Haizhao Yang, and Chao Zhou.  The discovery of dynamics via linear multistep methods and deep learning: Error estimation.  SIAM Journal on Numerical Analysis, 60(4):2014–2045, 2022.
  * [24] Xiaoli Chen, Liu Yang, Jinqiao Duan, and George Em Karniadakis.  Solving inverse stochastic problems from discrete particle observations using the fokker–planck equation and physics-informed neural networks.  SIAM Journal on Scientific Computing, 43(3):B811–B830, 2021.
  * [25] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In International conference on machine learning, pages 214–223. PMLR, 2017.
  * [26] Richard E Bellman.  Adaptive control processes.  In Adaptive Control Processes. Princeton university press, 2015\.
  * [27] Hannes Risken.  Fokker-planck equation.  In The Fokker-Planck Equation, pages 63–95. Springer, 1996.
  * [28] Lawrence C Evans.  Partial differential equations, volume 19.  American Mathematical Soc., 2010.
  * [29] Yaohua Zang, Gang Bao, Xiaojing Ye, and Haomin Zhou.  Weak adversarial networks for high-dimensional partial differential equations.  Journal of Computational Physics, 411:109409, 2020.
  * [30] Gang Bao, Xiaojing Ye, Yaohua Zang, and Haomin Zhou.  Numerical solution of inverse problems by weak adversarial networks.  Inverse Problems, 36(11):115003, 2020.
  * [31] Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Data-driven discovery of partial differential equations.  Science Advances, 3(4):e1602614, 2017.
  * [32] Herman Heine Goldstine.  A History of Numerical Analysis from the 16th through the 19th Century, volume 2.  Springer Science & Business Media, 2012.
  * [33] Kendall Atkinson, Weimin Han, and David E Stewart.  Numerical solution of ordinary differential equations, volume 108\.  John Wiley & Sons, 2011.
  * [34] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Multistep neural networks for data-driven discovery of nonlinear dynamical systems.  arXiv preprint arXiv:1801.01236, 2018.
  * [35] Rachael T Keller and Qiang Du.  Discovery of dynamics using linear multistep methods.  SIAM Journal on Numerical Analysis, 59(1):429–455, 2021.
  * [36] E Hairer SP Norsett and G Wanner.  Solving ordinary differential equations i: Nonsti problems, 1987.
  * [37] Michael Stein.  Large sample properties of simulations using latin hypercube sampling.  Technometrics, 29(2):143–151, 1987.
  * [38] Bernt Oksendal.  Stochastic differential equations: an introduction with applications.  Springer Science & Business Media, 2013.
  * [39] Ian Knowles and Robert J Renka.  Methods for numerical differentiation of noisy data.  Electronic Journal of Differential Equations, 21:235–246, 2014\.

## Appendix

##  Appendix A Method to compute temporal derivative

The right-hand side of equation (9) can be directly modified using the data
set $\mathbb{X}$. However, the left-hand side should be carefully considered
because of the temporal derivatives. There are several ways to handle the
temporal derivatives of $y(\mathbb{X}_{t})$ in the discrete-time snapshots,
where the most direct way is by the numerical scheme of the differentiation,
e.g., Euler’s scheme or other schemes with high order. However, due to the
high noise signal level of the data in the tasks of revealing stochastic
dynamics, these schemes may encounter failure because of the data glitch. One
practical solution is to use the schemes considering the smoothness. For
example, Tikhonov differentiation [39] given by

$\hat{y}^{\prime}=\left(A^{T}A+\lambda M^{T}M\right)^{-1}A^{T}y,$ (43)

where $\hat{y}^{\prime}$ is the approximation of the discrete form of
$y^{\prime}(t)$, $\lambda$ is a Tikhonov regularization, $A$ and $M$ come from
the approximations to the integral and the derivatives. Tikhonov
differentiation approximates the derivatives by the closeness of the integral
with the balance of the smoothness

$\hat{y}^{\prime}(t)=\operatorname{argmin}_{g}\left\|\int_{t_{0}}^{t}g(s)ds-\left(y(t)-y\left(t_{0}\right)\right)\right\|_{2}^{2}+\lambda\left\|\frac{dg}{dt}\right\|_{2}^{2}.$
(44)

Approximate the integral and the derivatives. One obtains the discrete form

$\hat{y}^{\prime}=\operatorname{argmin}_{g}\left\|Ag-(y-y_{0})\right\|_{2}^{2}+\lambda\left\|Mg\right\|_{2}^{2},$
(45)

where $A$ is a trapezoidal approximation to the integral and $M$ is a finite
difference approximation to the derivative.

Besides, the Adams-Moulton methods we used in this work is

$\left\\{\begin{gathered}y_{n}-y_{n-1}=hf_{n}\\\
y_{n+1}-y_{n}=\frac{h}{2}\bigg{(}f_{n+1}+f_{n}\bigg{)}\\\
y_{n+2}-y_{n+1}=h\bigg{(}\frac{5}{12}f_{n+2}+\frac{8}{12}f_{n+1}-\frac{1}{12}f_{n}\bigg{)}\\\
y_{n+3}-y_{n+2}=h\bigg{(}\frac{9}{24}f_{n+3}+\frac{19}{24}f_{n+2}-\frac{5}{24}f_{n+1}+\frac{1}{24}f_{n}\bigg{)}\\\
y_{n+4}-y_{n+3}=h\bigg{(}\frac{251}{720}f_{n+4}+\frac{646}{720}f_{n+3}-\frac{264}{720}f_{n+2}+\frac{106}{720}f_{n+1}-\frac{19}{720}f_{n}\bigg{)}\end{gathered}\right.$
(46)

##  Appendix B Discussion of the error

For a number of reasons, the Weak Collocation Regression method has a series
of errors. It is necessary to find out the source of the error to improve the
accuracy of our method. In general, the error of our method mainly comes from
the following aspects:

  1. 1.

The observation data of SDE

$dX_{t}=\mu dt+\sigma dW_{t}.$ (47)

Our experiment is inaccurate since the numerical error is contained in the
discretization scheme in the data generating process. In this research, we use
the Euler-Maruyama method to discrete the time $[0,T]$ with interval $\delta
t$, whose local truncation error is $O(\delta t^{2})$. In our numerical
experiments $\delta t=0.001$ or $\delta t=0.0001$ so the observation data
error here becomes $10^{-6}$ or $10^{-8}$. The initial distribution of SDE
also matters a lot because it determines whether the SDE fully evolves within
a specified time $[0,T]$.

  2. 2.

Inaccurate observation generates a random error in the measurement process,
which affects the identification of the drift and diffusion terms by
introducing an error in assembling matrix $A$ and vector $b$. For most nature
dynamics, the observation data is due to the restriction of the experimental
measurement and the measurement instrument. In this research, we show the
robustness of our methods by adding white noise to observation data. The
result shows the high insensitivity of our methods to noise. We can introduce
the adaptive regularization parameter adjustment with respect to the noise
level in the sparse regression process in the traditional Inverse Problem way
for further improvement.

  3. 3.

There exists error in computing the integration or expectation, such as
$\int_{\Omega}p(x,t)\boldsymbol{\mu}(\boldsymbol{x})\cdot\nabla\psi(\boldsymbol{x})dx$
or $\mathbb{E}_{\boldsymbol{x}\sim
p(\boldsymbol{x},t)}[\boldsymbol{\mu}(\boldsymbol{x})\cdot\nabla\psi(\boldsymbol{x})]$.
In this research, we use the Monte Carlo method to compute the integration,
whose convergence error in probability is $O(1/\sqrt{N})$:

$\lim_{N\to+\infty}\mathbb{P}\\{|I_{N}-\int_{\Omega}g(x)dx|<\frac{\varepsilon\mbox{var}(g(x))}{\sqrt{N}}\\}=\frac{1}{\sqrt{2\pi}}\int_{-\varepsilon}^{\varepsilon}e^{-\frac{t^{2}}{2}}dt$
(48)

where $N$ is the sample number, and $I_{N}$ is the average of integrand among
$N$ samples. We can find that the error here plays a decisive role in our
experimental results.

  4. 4.

The Linear Multistep Methods to deal with the temporal derivative also
produces an error. Given data of discrete time snapshots, the LMM provides
only an approximation of their temporal derivative with residual errors. Some
analyses and results already exist. Concerning the local truncated error, the
M-step Adams-Moulton methods have an order of error M+1, while the M-step BDF
has an error of order M. Moreover, the M-step BDF scheme and 0-step Adams-
Moulton method are stable, while 1 step Adams–Moulton method is marginally
stable. Hence the Adams methods of trapezoidal rule and BDF2 are both
convergent. For the most rigorous conception and theorem, we refer to [23,
35].

  5. 5.

After building linear systems, there exists an error in solving the parameter
in the drift and diffusion term. In this research, we use linear regression
with regularization to compute $A\theta=y$. If matrix $A,b$ is accurate
enough, then $\theta$ is also accurate for the equation is over-determined.
While the above three errors are objective, solving $A\theta=y$ also produces
an error that is related to the conditional number of $A$.