# Data-Driven Denoising of Accelerometer Signals

Daniel Engelsman and Itzik Klein D. Engelsman is with the Hatter Department of
Marine Technologies, University of Haifa, Israel.I. Klein is with the Hatter
Department of Marine Technologies, University of Haifa, Israel (e-mail:
kitzik@univ.haifa.ac.il).

###### Abstract

Modern navigation solutions are largely dependent on the performances of the
standalone inertial sensors, especially at times when no external sources are
available. During these outages, the inertial navigation solution is likely to
degrade over time due to instrumental noises sources, particularly when using
consumer low-cost inertial sensors. Conventionally, model-based estimation
algorithms are employed to reduce noise levels and enhance meaningful
information, thus improving the navigation solution directly. However,
guaranteeing their optimality often proves to be challenging as sensors
performance differ in manufacturing quality, process noise modeling, and
calibration precision. In the literature, most inertial denoising models are
model-based when recently several data-driven approaches were suggested
primarily for gyroscope measurements denoising. Data-driven approaches for
accelerometer denoising task are more challenging due to the unknown gravity
projection on the accelerometer axes. To fill this gap, we propose several
learning-based approaches and compare their performances with prominent
denoising algorithms, in terms of pure noise removal, followed by stationary
coarse alignment procedure. Based on the benchmarking results, obtained in
field experiments, we show that: (i) learning-based models perform better than
traditional signal processing filtering; (ii) non-parametric kNN algorithm
outperforms all state of the art deep learning models examined in this study;
(iii) denoising can be fruitful for pure inertial signal reconstruction, but
moreover for navigation-related tasks, as both errors are shown to be reduced
up to one order of magnitude. Both our code and dataset are publicly available
@ GitHub.

###### Index Terms:

Inertial sensors, MEMS IMU, denoising, deep learning, kNN regression, attitude
estimation.

##  I Introduction

Inertial Navigation System (INS) is one of the most commonly used navigation
systems. It consists of an inertial measurement units (IMU) with three
orthogonal gyroscopes and accelerometers to determine the platform position,
velocity, and orientation [1, 2]. With recent technological rise of micro-
electro-mechanical-systems (MEMS) inertial sensors, their integration became
very common in many applications, platforms, and environments (air, ground,
sea). Their popularity can be explained due to their small size, high cost-
effective, low power consumption, and low-cost prices. However, when external
position/velocity update source is not available, the INS navigation solution
becomes solely dependent on the inertial sensor performance and its error
regime. Integration of noisy measurements propagates into the navigation
solution, resulting in a solution drift [3]. This problem worsens when low
cost MEMS-IMU are used, due to their characteristic noise and inherent bias
[4].
Conventionally, INS filtering does not incorporate pure signal reconstruction.
Instead, a fusion mechanism determines the weight given to current state
estimates and external measurements, followed by integration into kinematic
quantities. But, given extreme nonlinear dynamics, high frequency noises and
non-Gaussian distribution of random noise variables, the filtering process is
contaminated. To that end, denoisers can be used for pre-filtering, thus
covering holes in the model. However, when deterministic noise sources are
introduced and a thorough calibration process cannot be taken, performance of
conventional denoisers becomes limited. This is where data-driven methods come
in handy, as trainable parameters (weights and biases) are capable of
extracting the meaningful patterns from the body of observations.

According to our literature review, among all data-driven denoising works,
only one work addresses accelerometer denoising, as the rest perform
gyroscopes denoising. It is not surprising in light of the instrumental
differences, as same physical scenario is projected inherently different over
each sensor. For example, while stationary, MEMS gyroscopes (insensitive to
Earth’s rotation rate) are expected to output a three-dimensional zero vector
$\|\boldsymbol{\omega}\|\approx 0$, regardless body orientation (in practice
we get sensor noise). Contrarily, MEMS accelerometers sense Earth gravitation
depending on a given orientation, such that their axes being regularly
subjected to a larger norm $\|\boldsymbol{f}\|\approx\text{g}$. The constant
presence of non-zero amplitudes, wider dynamic ranges and larger noise
densities [5], are ultimately reflected in the scale and variability of the
samples. Additionally, these differences also propagate into the calibration
process. While gyroscopes calibration is self-contained as output offset can
be simply subtracted from zero, accelerometers require auxiliary attitude
information to ensure parallelism to the gravity axis. When data-driven models
are required to generalize over these patterns, distributions of specific
force measurements exhibit significantly higher levels of variance and bias,
making accelerometers denoising much more challenging, thus still unexplored.

To fill this gap, we propose several data-driven models which successfully
manage to denoise accelerometer signals, improving dramatically, not only its
signal-to-noise ratio, but also subsequent navigation tasks. The main
contributions of this paper can be summarized as follows:

  1. 1.

An exhaustive literature review of MEMS-IMU denoising techniques, with respect
to sensors and dynamics.

  2. 2.

A bi-directional implementation of a recurrent neural network architecture.

  3. 3.

A naive yet competitive implementation of a k-nearest neighbor algorithm.

  4. 4.

A navigation-related metric, demonstrating how denoising improves stationary
coarse alignment (SCA).

The rest of the paper is organized as follows: Section II gives an up-to-date
literature review of MEMS-IMU denoising methods. Section III describes
auxiliary methods used to assess our strategy and section IV presents our
proposed solution. Section V elaborates the data preparation process, section
VI shows analysis and results and section VII gives conclusions.

##  II Literature Review

Study of inertial sensor denoising techniques dates back to the late 90s, and
can be inclusively divided into conventional signal processing approaches and
to recent learning-based approaches. Followed is an in-depth literature review
to conveniently map the problem domain. Some works evaluate their proposed
denoisers in terms of reconstruction of the original IMU output, while others
examine its contribution as a prefilter upon a wider endgame output (e.g.
Euler angles, navigation states). Table I presents conventional signal
processing (SP-based) approaches whose analysis and synthesis is based on the
signal structure and nature, enabling to detect components of interest.

TABLE I: Summary of SP-based denoisers Algorithm | Ref. | Gyro | Acc. | Stationary | Dynamic
---|---|---|---|---|---
Moving average | [6] | ✓ |  | ✓ | ✓
ARMA | [7, 8, 9, 10, 11] | ✓ |  | ✓ |
[12, 13, 14] | ✓ | ✓ | ✓ |
[15] | ✓ | ✓ |  | ✓
EMD | [16] | ✓ | ✓ | ✓ | ✓
[17, 18] | ✓ |  | ✓ |
[19, 20, 21] | ✓ |  | ✓ | ✓
Fuzzy logic | [22, 23] | ✓ | ✓ | ✓ | ✓
Savitsky-Golay | [24] |  | ✓ | ✓ |
[25, 26] | ✓ | ✓ | ✓ | ✓
[27] | ✓ |  | ✓ | ✓
Wavelets | [28, 29] | ✓ | ✓ |  | ✓
[30] | ✓ | ✓ | ✓ |
[31, 32, 33] | ✓ |  | ✓ |
[34, 35] | ✓ |  |  | ✓
[36, 37] | ✓ | ✓ |  | ✓

Moving average (MA) techniques can be used as efficient smoothing filter,
based on errors (residuals) from previous forecasts [6]. Other works
elaborated this by combining a weighted regression term over the lagged
values, namely auto-regressive moving-average (ARMA) [7, 8, 9, 10, 11, 12, 13,
14, 15]. Empirical mode decomposition (EMD) is a time frequency analysis which
decomposes multicomponent signals into a finite number of Intrinsic Mode
Functions (namely, its building blocks). This way, local characteristic time
scale is emphasized such that non-stationary and nonlinear signals can be
robustly handled [16, 17, 18, 19, 20, 21]. In fuzzy logic, quantified
statements are used to determine objects level of membership in values ranging
from zero to one, instead of true or false. During inference stage, the
formulated logic is used to map elements from a given input space to output
space, performing well with high levels of uncertainties in nonlinear dynamics
[22, 23]. Savitsky-Golay is a smoothing filter that fits an optimal local
curve over a moving window size, using low-degree polynomial regression [24,
25, 26, 27]. And the wavelet-based algorithms, currently most popular for
MEMS-IMU denoising, uses a window function which is scaled and shifted in the
signal time axis, projecting it into a set of basis functions named wavelets
[28, 29, 30, 31, 32, 33, 38, 34, 35, 36, 37].
Table II presents a group of learning-based algorithms, which gained much
popularity in recent years in a growing number of fields, including inertial
sensors and autonomous navigation. Their advantage lies in their ability to
identify complex patterns by learning high-level features of the data, thus
diminishing the need in domain expertise.

TABLE II: Summary of learning-based denoisers Algorithm | Ref. | Gyro | Acc. | Stationary | Dynamic
---|---|---|---|---|---
Linear regression | [39] | ✓ | ✓ | ✓ | ✓
CNN | [40] | ✓ |  | ✓ | ✓
RNN | [41, 42, 43] | ✓ |  | ✓ |
LSTM | [44, 45] | ✓ |  | ✓ |
[46] | ✓ |  | ✓ | ✓
GRU | [47] | ✓ |  | ✓ |

Gonzalez proposed a multiple linear regression (MLR) model which uses several
vector-valued variables to predict the outcome of a dependent variable [39].
Bossard proposed a convolutional neural network (CNN) which computes gyro
corrections for the undesirable noise, before integrated into orientation
increments [40]. The rest of the references utilized a common extension to
feedforward neural networks called recurrent neural networks (RNN), which were
designed to handle variable-length sequential signals. Unlike CNN, which
excels at finding spatial relations over grid-like topology, their RNN
extension has a feedback element which enables forward and backwards
connections such that complex dynamic relationships over distant time steps
are better captured [41, 42, 43]. However, despite its ability to process
temporal information of any length, it is limited with long-term dependencies
due to the exponential decay of the loss function, namely vanishing gradient
problem. To that end, an improved versions were introduced, where additional
control units (gates) were added to allow better flow of the gradient and to
maintain memory over long time periods, namely long short term memory (LSTM)
[44, 45, 46] and gated recurrent unit (GRU) [47].

As already discussed above, it can be seen that the majority of the works
(75$\%$) address gyroscopes denoising, especially in the learning-based
approaches.

##  III Problem Formulation

After describing the big picture, this section presents two functional methods
used for benchmarking (A.) and evaluating (B.) the validity of our proposed
approaches.

###  III-A Signal Processing Denoising Approaches

Conventional SP-based denoisers perform noise reduction by either spatial
smoothing, local regression, or by imposing spectral constraints to filter out
unwanted frequencies. Following are three reference methods, used for
comparison with the data-driven approaches. Their optimality was obtained by
brute-force search over the training set.

####  III-A1 Moving Average (MA)

suppresses unstable signal noises by averaging measurements inside a rolling
window [6]

$\displaystyle\hat{\boldsymbol{x}}_{\textbf{{MA}},i}=\frac{1}{T}\sum_{t=0}^{T-1}\boldsymbol{x}_{i+t-T}\quad\forall\
i\geq T$ (1)

where $\boldsymbol{x}_{i}$ is a noisy sample and $T$ is the window size.

####  III-A2 Savitsky-Golay (SG)

allows denoising by deriving observations directly from time domain, thus
avoiding spectral decomposition [24, 25, 26, 27]. By fitting successive sets
of adjacent points with a low-degree polynomial, followed by least-squares
regression, local noise is smoothed out. Noisy signal $\boldsymbol{x}_{i}$ is
replaced with a set of $m$ convolution coefficients $C_{t}$,

$\displaystyle\hat{\boldsymbol{x}}_{\textbf{{SG}},i}=\sum_{t=-t_{s}}^{t_{s}}C_{t}\boldsymbol{x}_{i+t}\
,\quad\tfrac{m+1}{2}\leq t\leq t_{f}-\tfrac{m-1}{2}$ (2)

where $t_{f}$ is the duration of data point $x_{t}$ and $t_{s}=\frac{m-1}{2}$
denotes the window margins.

####  III-A3 Discrete wavelet transform (DWT)

the most common denoising technique, where input signal $\boldsymbol{x}_{i}$
is represented in both time and frequency domains, by decomposing it into a
set of basis functions [28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. Here, a
Daubechies (db4) mother wavelet function ($\psi$) is used with different
scaling ($a=2$) and shifting ($b$) parameters, providing a progressively finer
outputs, given by coefficients matrix $\Psi$

$\displaystyle\Psi_{i}[b,a^{j}]$
$\displaystyle=\sum_{t=0}^{N-1}\boldsymbol{x}_{i}[t]\frac{1}{\sqrt{a^{j}}}\psi_{j}\left(\frac{t-b}{a^{j}}\right)$
(3)

Since small valued coefficients are dominated by noise, hard thresholding
$\text{T}_{\text{Hard}}$ is used to remove them, thus preserving only
meaningful information. Then, an inverse transform is applied on the
thresholded wavelet coefficients, to reconstruct the denoised matrix back to
time domain signal as given by

$\displaystyle\hat{\boldsymbol{x}}_{\textbf{{DWT}},i}=\Big{(}\text{T}_{\text{Hard}}(\Psi_{i})\Big{)}^{-1}$
(4)

###  III-B Inertial Navigation Coarse alignment

In addition to analyzing our proposed approach denoising capabilities on the
accelerometer signals, we examine its influence on a stationary coarse
alignment (SCA) procedure. To that end, this section gives a brief
introduction to the coarse alignment theory. The attitude (roll and pitch
angles) can be determined using accelerometer measurements. The orientation of
body coordinates frame with respect to navigation frame is given by a
transformation matrix, represented by a product of three successive rotations
about the z-y-x axes, corresponding to Euler angels: yaw ($\psi$), pitch
($\theta$) and roll ($\phi$).

$\textbf{T}_{b}^{n}=\begin{bmatrix}\text{c}_{\theta}\text{c}_{\psi}&\text{s}_{\phi}\text{s}_{\theta}\text{c}_{\psi}-\text{c}_{\phi}\text{s}_{\psi}&\text{c}_{\phi}\text{s}_{\theta}\text{c}_{\psi}+\text{s}_{\phi}\text{s}_{\psi}\\\
\text{c}_{\theta}\text{s}_{\psi}&\text{s}_{\phi}\text{s}_{\theta}\text{s}_{\psi}+\text{c}_{\phi}\text{c}_{\psi}&\text{c}_{\phi}\text{s}_{\theta}\text{s}_{\psi}-\text{s}_{\phi}\text{c}_{\psi}\\\
-\text{s}_{\theta}&\text{s}_{\phi}\text{c}_{\theta}&\text{c}_{\phi}\text{c}_{\theta}\end{bmatrix}$
(5)

In stationary conditions, accelerations in navigation frame are equal to zero
($\dot{\textbf{v}}^{n}=\textbf{0}$), thus only the gravity vector is projected
upon the accelerometer axes [3]

$\displaystyle\boldsymbol{f}_{ib}^{b}=\begin{bmatrix}f_{ib,x}^{b}\\\
f_{ib,y}^{b}\\\
f_{ib,z}^{b}\end{bmatrix}=-\textbf{T}_{n}^{b}\,\textbf{g}^{n}=\begin{bmatrix}s_{\theta}\\\
-s_{\phi}c_{\theta}\\\ -c_{\phi}c_{\theta}\end{bmatrix}\text{g}$ (6)

Using analytical coarse alignment, leveling can be computed by the axial
components of the specific force vector $\boldsymbol{f}_{ib}^{b}$ [4]

$\displaystyle\phi$
$\displaystyle=\arctan_{2}\Big{(}-f_{ib,y}^{b}\,,\,-f_{ib,z}^{b}\Big{)}$ (7)
$\displaystyle\theta$
$\displaystyle=\arctan\Big{(}\frac{-f_{ib,x}^{b}}{\sqrt{f_{ib,y}^{b\
2}+f_{ib,z}^{b\ 2}}}\Big{)}$ (8)

The accuracy of the SCA process affects the whole INS performance, especially
during pure inertial navigation. To increase the degraded signal-to-noise
ratio (SNR), instrumental errors can be reduced in proportion to a square root
factor, by either averaging over $n$ identical sensors ($1/\sqrt{n}$), or by
averaging over $T$ time steps ($1/\sqrt{T}$). Thereby, effective denoising can
reduce averaging times of the SCA procedure.

##  IV Proposed Approach

We propose implementing and modifying several learning algorithms for the
accelerometer denoising problem, which unlike conventional, model-based,
denoising filters, are independent in any analytic error model. Instead, they
learn to map noisy readings into their corresponding accurate measurements.
Figure 1 presents the proposed approach, where noisy measurements from
commercial grade inertial sensors are input to a learning model for the
denoising task. The model then outputs denoised accelerometer readings, which
are used to calculate an approximation error, with respect to a high-end,
accurate inertial sensors readings, as ground-truth (GT). The error is first
used to optimize the model during the training phase. Then, during testing
phase (VI), the well-trained models are evaluated over unseen noisy samples.

Figure 1: Denoising mapping function
$f:\boldsymbol{x}\rightarrow\hat{\boldsymbol{x}}$

Unlike SP-based methods, here extracted patterns are associated with target
functions, enabling signal denoising by generalization. Similarly to human
processing, feedforward neural networks (FNN) excel at pattern recognition, as
their layered structure perform hierarchical representations, allowing
extraction of spatial features [48]. However given sequential relations, their
node-specific weights eventually fail to learn, as lengthy inputs and complex
ordering impose an exponentially growing number of parameters. To that end, a
recurrent feedback mechanism was proposed [49]

Figure 2: One-layer recurrent neural network (RNN)

Three dimensional measurements are processed inside an encoding cell called
hidden state, using a weight vector of $m$ learnable parameters, such that
$\overrightarrow{\textbf{{h}}}\in\mathbb{R}^{3\times m}$. Unlike FFN, here
same weights are shared across the entire sequence, and current inputs are
taken from previous outputs. This way, the model avoids spatial memorization
and focuses on intertemporal dependencies. Figure 2 illustrates a sequence-to-
sequence mapping, consists of single 3D noisy measurements, estimated
(denoised) every time step $\hat{x}_{t}\ \in\mathbb{R}^{3}$.
Following, we elaborate on our proposed models: three variations of deep
recurrent neural networks (RNN) and one naive machine learning model (kNN).

###  IV-A Unidirectional bi-layer LSTM

The architecture of the first model is similarly arranged as shown in Figure
2, but consists of two (bi) stacked layers. The hidden states are an RNN
variant called long-short term memory (LSTM) [50], originally developed to
handle the decay of the loss function:

Figure 3: Hidden state unit of an LSTM

As illustrated in Figure 3, the LSTM hidden state contains three non-linear
gates, used to regulate the signal flow. The forget gate $\textbf{{f}}_{t}$
decides whether ignoring or adding new (short-term) information to the cell
state $\textbf{{C}}_{t}$ memory. The input gate $\textbf{{i}}_{t}$, is then
used to estimate the relevance of current input to past hidden states. This
way, the cell state acts as a global memory unit, where distant dependencies
(long-term) are maintained across time

$\displaystyle\textbf{{f}}_{t}$
$\displaystyle=\sigma(\textbf{W}_{f}[x_{t},\,\overrightarrow{h}_{t-1}]+\textbf{b}_{f})$
(9) $\displaystyle\textbf{{i}}_{t}$ $\displaystyle=\sigma(\textbf{W}_{i}\hskip
1.99168pt[x_{t},\,\overrightarrow{h}_{t-1}]+\textbf{b}_{i})$ (10)
$\displaystyle\tilde{\textbf{{C}}}_{t}$
$\displaystyle=\tanh(\textbf{W}_{c}[x_{t},\,\overrightarrow{h}_{t-1}]+\textbf{b}_{c})$
(11) $\displaystyle\textbf{{C}}_{t}$
$\displaystyle=\textbf{{f}}_{t}\otimes\textbf{{C}}_{t-1}+\textbf{{i}}_{t}\otimes\tilde{\textbf{{C}}}_{t}$
(12)

Operator $\otimes$ denotes an element-wise multiplication, $\textbf{W}_{j}$
and $\textbf{b}_{j}$ are the weights and biases respectively, and the sigmoid
function scales encoded vectors as follows
$\sigma:\mathbb{R}^{m}\rightarrow(0,1)^{m}$. Next, the output gate
$\textbf{{o}}_{t}$ is fused with the updated cell state

$\displaystyle\textbf{{o}}_{t}$ $\displaystyle=\sigma(\textbf{W}_{o}\hskip
1.42262pt[x_{t},\,\overrightarrow{h}_{t-1}]+\textbf{b}_{o})$ (13)

such that the next hidden state is determined as follows:

$\displaystyle\overrightarrow{\textbf{{h}}_{t}}$
$\displaystyle=\textbf{{o}}_{t}\otimes\tanh(\textbf{{C}}_{t})$ (14)

Finally, the hidden state signal is taken for state prediction, using a linear
activation function $\boldsymbol{\phi}$ (neuron) which scales it back into a
continuous and unbounded range, given by

$\displaystyle\hat{\boldsymbol{x}}_{t}$
$\displaystyle=\boldsymbol{\phi}(\overrightarrow{\textbf{{h}}_{t}})=\textbf{W}_{\hat{y}}^{\mathsf{T}}\overrightarrow{\textbf{{h}}_{t}}+\textbf{b}_{\hat{y}}\
\in\ \mathbb{R}^{3}$ (15)

###  IV-B Bi-directional one-layer RNN

Conventionally, sequential inputs are analyzed along the positive time
direction (chronologically), thus limiting the learning with respect to only
past states. Figure 4 illustrates a two opposing cells architecture, where
future contexts are learned from the opposite time direction (anti-
chronological):

Figure 4: Bi-directional RNN [51]

Although RNNs lack the cell state concept, their simple parametric relations
managed to exhibit competitive results when applied in a shallow (one-layer)
arrangement. Figure 5 presents a forward directed hidden unit, at time step
$t$:

Figure 5: RNN hidden state and its denoised output

Here, two different hidden layers are allocated to process data both in
forward $\overrightarrow{\textbf{{h}}_{t}}$ and backward
$\overleftarrow{\textbf{{h}}_{t}}$ directions

$\displaystyle\overrightarrow{\textbf{{h}}_{t}}$
$\displaystyle=\sigma(\textbf{W}_{\overrightarrow{h}}[x_{t},\,\overrightarrow{h}_{t-1}]+\textbf{b}_{\overrightarrow{h}})$
(16) $\displaystyle\overleftarrow{\textbf{{h}}_{t}}$
$\displaystyle=\sigma(\textbf{W}_{\overleftarrow{h}}[x_{t},\,\overleftarrow{h}_{t-1}]+\textbf{b}_{\overleftarrow{h}})$
(17)

Then, parameters of the opposing hidden states are stacked
$\overline{\textbf{{h}}}_{t}$, and fed into the output layer to provide state
denoising

$\displaystyle\overline{\textbf{{h}}}_{t}$
$\displaystyle=[\overrightarrow{\textbf{{h}}_{t}},\overleftarrow{\textbf{{h}}_{t}}]\
\in\ \mathbb{R}^{3\times 2m}$ (18) $\displaystyle\hat{\boldsymbol{x}}_{t}$
$\displaystyle=\boldsymbol{\phi}(\overline{\textbf{{h}}}_{t})=\textbf{W}_{\hat{y}}\overline{\textbf{{h}}}_{t}+\textbf{b}_{\hat{y}}\
\in\ \mathbb{R}^{3}$ (19)

###  IV-C Bi-directional one-layer GRU

The third model is also arranged in a bi-directional form (Figure 4), but
utilizes a different gating mechanism called gated recurrent unit (GRU), which
merges both input and forget gates into a single update gate [52]. Using only
two-gates, it requires shorter training epochs, less computational efforts and
is thus more robust to the vanishing gradient problem. Using a sigmoid
function, the gates regulate the trade-off between previous hidden states and
new input information. Zero output means closed gate, such that historical
data cannot pass beyond, and only the new input is emphasized, as presented in
Figure 6.

Figure 6: Hidden state unit of a GRU model

The reset gate $\textbf{{r}}_{t}$ controls the dominance of past states with
respect to new input, and the update gate $\textbf{{z}}_{t}$ specifies how
much memory will pass to the next hidden state.

Since both hidden states are symmetrically opposite, only the chronological
direction is formulated to simplify notation

$\displaystyle\textbf{{r}}_{t}$
$\displaystyle=\sigma(\textbf{W}_{r}[x_{t},\,\overrightarrow{h}_{t-1}]+\textbf{b}_{r})$
(20) $\displaystyle\textbf{{z}}_{t}$
$\displaystyle=\sigma(\textbf{W}_{z}[x_{t},\,\overrightarrow{h}_{t-1}]+\textbf{b}_{z})$
(21)

where the chronological hidden state is defined as

$\displaystyle\tilde{\textbf{{h}}}_{t}$
$\displaystyle=\tanh(\textbf{W}_{\overrightarrow{h}}[x_{t},\
\textbf{{r}}_{t}\otimes\overrightarrow{h}_{t-1}]+\textbf{b}_{\overrightarrow{h}})$
(22) $\displaystyle\overrightarrow{\textbf{{h}}}_{t}$
$\displaystyle=(1-\textbf{{z}}_{t})\otimes\overrightarrow{\textit{h}}_{t-1}+\textbf{{z}}_{t}\otimes\tilde{\textbf{{h}}}_{t}$
(23)

Similarly to (19), the opposite cells are horizontally stacked and
approximated by the external output layer as

$\displaystyle\overline{\textbf{{h}}}_{t}$
$\displaystyle=[\overrightarrow{\textbf{{h}}_{t}},\overleftarrow{\textbf{{h}}_{t}}]\
\in\ \mathbb{R}^{3\times 2m}$ (24) $\displaystyle\hat{\boldsymbol{x}}_{t}$
$\displaystyle=\boldsymbol{\phi}(\overline{\textbf{{h}}}_{t})=\textbf{W}_{\hat{y}}\overline{\textbf{{h}}}_{t}+\textbf{b}_{\hat{y}}\
\in\ \mathbb{R}^{3}$ (25)

Note that the optimality of the above models was obtained heuristically, with
respect to hidden state architecture, number of parameters, directionality of
signal flow and model depth.

###  IV-D k-nearest neighbors algorithm (kNN)

Our last proposed model does not assume priors on the underlying distribution,
but is only determined by the number of $k$ nearest neighbors. Let the dataset
be a set of $n$ pairs

$\displaystyle\mathcal{D}=\\{(\boldsymbol{x}_{i},\boldsymbol{x}_{{}_{GT},i})\\}_{i=1}^{n}$
(26)

where denoised outputs are determined by the mean GT value of $k$ closest
(noisy) samples to the new data point

$\displaystyle\hat{\boldsymbol{x}}(k)=\frac{1}{k}\sum_{\boldsymbol{x}_{i}\in\boldsymbol{N}(k)}\boldsymbol{x}_{{}_{GT},i}$
(27)

Given r-dimensional data points $\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{2}$,
neighborhood $\boldsymbol{N}$ will be consisted of $k$ closest neighbors, as
proximity is calculated using an Euclidean distance function

$\displaystyle
d(\boldsymbol{x}_{1},\boldsymbol{x}_{2})=\sqrt{\sum_{j=1}^{r}(x_{1,j}-x_{2,j})^{2}}$
(28)

To simplify demonstration, Figure 7 illustrates a kNN denoising (k=6), between
one-dimensional data points, but applies similarly to higher dimensional
vectors as in our experiment

Figure 7: Non-parametric kNN denoising

Each coordinate in the acceleration plane, refers to a different sensor grade.
The horizontal and vertical axes refer to the noisy and GT measurements
respectively, where the colorbar indicates only the GT values. Estimating a
new data point $\hat{\boldsymbol{x}}$, is obtained by averaging all GT data
points inside $\boldsymbol{N}(k)$. In the absence of pre-defined statistical
methods to find optimality, an heuristic search brought up the following
optimal $k$

$\displaystyle
k^{*}=\arg\min_{k}\|\hat{\boldsymbol{X}}(k)-\boldsymbol{X}_{GT}\|^{2}\approx\sqrt{n/5}$
(29)

where $\hat{\boldsymbol{X}}$ is the kNN estimates and $\boldsymbol{X}_{GT}$ is
their corresponding GT samples, both taken from the test-set, which is one
fifth of the dataset ($n=90,000$).

###  IV-E Loss function

Unlike the kNN optimization, deep learning architectures use a loss function
to assess deviations between their denoised predictions and the actual ground
truth. The axial error of a single sample containing $H$ time steps is given
by

$\displaystyle\boldsymbol{e}_{j}=\hat{\boldsymbol{x}}_{j}-\boldsymbol{x}_{j,_{GT}}\
\in\ \mathbb{R}^{H}$ (30)

where $j$ denotes an axis index as each sensor measurement consists of three
channels - x,y and z. Then, during training, performance is assessed using
mean squared error (MSE), such that loss surface is continuously
differentiable, and errors can be minimized progressively as a function of the
model weights

$\displaystyle\text{MSE}=\frac{1}{H\times
3}\sum_{j\in\\{\text{x,y,z}\\}}\sum_{i=1}^{H}(e_{i,j})^{2}$ (31)

##  V Datasets and Error Measures

To evaluate our methodology, models from both signal processing-based and
learning-based approaches are evaluated and compared over simulated and
experimental datasets, using unique error measures for the evaluation process.

###  V-A Simulated Dataset

In stationary conditions, the gravity vector is projected into the specific
force measurements of the accelerometers. This projection can be represented
using Euler angles.
In our simulative setup, Euler angles are drawn within a bounded domain, and
using small intervals, a rich combination of data points (orientations) is
guaranteed. Using the transformation matrix (5), Euler angles dictate how
gravity is projected on the sensor axes (6). Without loss of generality, the
yaw angle has no influence on the gravity projection, thus it remains constant
($\psi_{sim.}=0$), and we address only the roll and pitch angles in the range
of

$\displaystyle-15^{\circ}\,\leq\,$ $\displaystyle\phi_{sim.}\leq+15^{\circ}$
(32) $\displaystyle-15^{\circ}\,\leq\,$
$\displaystyle\theta_{sim.}\leq+15^{\circ}$ (33)

Both angles are divided into intervals of $0.1$ deg, such that
$\frac{15-(-15)}{0.1}=300$ increments are obtained. Then, the angles are
combined together, and the overall number of simulated instances becomes
3002=90,000 combinations. Notice, that in this study we limit the angles range
to $\pm 15^{\circ}$ to enable model training in fair times on our hardware
(Intel i5-9600K CPU @ 3.70 GHz and NVIDIA GTX2080 GPU). Yet, of course, the
proposed methods can be applied in any angle range, influencing only the
training time. The inference mode (real-time) operation is not influenced by
the dataset size and the training time. Next, synthetic noise and error terms
are added to each accelerometer channel using the following error terms:

  1. 1.

Velocity random walk (VRW) - accumulative error due to white noise in the
measurement.

  2. 2.

Bias instability (BI) - stationary stochastic process which act as a low-order
Gauss-Markov process.

  3. 3.

Bias offset (BO) - a constant offset caused by misprojection of the true
accelerations.

Figure 8 illustrates the simulation, showing that the constructed dataset
contains both simulated GT and measured accelerometer readings. Table III
specifies the synthetic additive noise values (Simulation column).

Figure 8: Simulated dataset generation procedure

###  V-B Field Experiment Dataset

The stationary conditions field experiment dataset is based on a dataset
published in [53]. There, a unique device, as shown in Figure 9 was built to
align between a Huawei P40 smartphone and an Inertial Lab MRU-P unit [54].

Figure 9: Unique device aligning between a Huawei P40 smartphone and an
Inertial Lab MRU-P unit.

In this setup, the smartphone accelerometer readings are used as unit under
test while the MRU’s accelerometer readings serve as GT measurements. Both
accelerometers provided measurements at 100Hz sampling rate.
Due to practical constraints, smaller amount of raw measurements were taken in
much bigger intervals of $\approx 3.0^{\circ}$ deg

$\displaystyle-15^{\circ}\,\leq\,$ $\displaystyle\phi_{exp.}\leq+15^{\circ}$
(34) $\displaystyle-15^{\circ}\,\leq\,$
$\displaystyle\theta_{exp.}\leq+15^{\circ}$ (35)

The experimental dataset consists of $10^{2}=100$ recordings at different
orientations and duration, divided into sub-samples depending on the window
size. Unlike the simulated dataset, yaw angles here are non-zeros
($\psi_{exp.}\neq 0$), as they were taken from a wide variety of platform
orientations. Table III specifies the sensors noise specs, as stated by
manufacturers (Experiment column).

TABLE III: Error sources used in the simulation and experiments | Simulation | Experiment
---|---|---
Error | Units | Noisy | GT | P40 | MRU-P (GT)
VRW | [m/s/$\sqrt{s}$] | 0.005 | 0.00001 | 0.003 | 0.00025
BI | [m/s2] | 0.001 | 0.00001 | 0.001 | 0.00005
BO | [m/s2] | 0.05 | 0.00001 | 0.067 | 0.0001

Lastly, to increase variability of the experimental dataset and avoid model
overfitting, each sample underwent two consecutive transformations suitable
for time-series [55, 56, 57]: i) Angular augmentation - applying angular
transformations over a given orientation, to obtain more different samples and
densify the sparse distribution. ii) Noise augmentation - applying stochastic
error in form of additive white Gaussian noise. This way, the augmented
dataset also contains a total of 90,000 samples, improving generalizability
and noise robustness of the learning-based models.

###  V-C Performance Metrics

During the testing phase, several evaluation functions are used to provide
relevant measures of the models performances from the perspective of signal
reconstruction.

  * •

Root mean squared error (RMSE) - is defined as the square root of the MSE
(31), given by

$\text{RMSE}=\sqrt{\text{MSE}}$ (36)

  * •

Mean absolute error (MAE) - returns an average magnitude of the residuals
(30), calculated over $n$ time steps, disregarding their direction, given by

$\text{MAE}=\frac{1}{n}\sum_{i=1}^{n}|e_{i}|$ (37)

  * •

Peak signal-to-noise ratio (PSNR) - returns the ratio between maximum signal
value (MAX) and noise levels. However here, PSNR expresses ratio between two
amplitudes (root-power quantity), thus a measure of order of magnitude is
equal to 20 [dB] difference, given by

$\text{PSNR}=20\log_{10}\left(\frac{\text{MAX}}{\text{RMSE}}\right)$ (38)

  * •

Relative absolute error (RAE) - returns a relative measure of the discrepancy
between the residuals and their corresponding ground truth mean value
$\mu_{{}_{GT}}$, given by

$\text{RAE}=\frac{\sum_{i=1}^{n}|e_{i}|}{\sum_{i=1}^{n}|x_{{}_{GT},i}-\mu_{{}_{GT}}|}$
(39)

For example, Table IV presents several error metrics, where outputs of a
denoising model are compared to raw accelerometer measurements, in terms of
dissimilarity to the accurate GT measurements. Therefore, the optimization
task is $\hat{x}\approx x_{GT}$.

TABLE IV: Toy Example: Evaluation | RMSE [m/s2] | MAE [m/s2] | PSNR [dB] | RAE [$\%$]
---|---|---|---|---
Model | 0.0251 | 0.0198 | 51.8317 | 0.4436
Noisy | 0.0626 | 0.0464 | 43.8681 | 1.1096

It can be seen that noise levels are reduced by several aspects and peak
signal-to-noise ratio is improved by 8 [dB].

##  VI Analysis and Results

This section presents a comparative performance analysis between all models
introduced above: signal processing-based (III-A) and our proposed learning-
based (IV). Both simulated and the experimental datasets are assessed, first
in terms of pure signal reconstruction, followed by measuring its influence on
the stationary coarse alignment procedure.

###  VI-A Simulation Assessment

After the simulated dataset is generated, models are trained on the training
set before comparing performances. Then, the simulated test-set utilizes to
examine the models capability of reducing noise sources, in terms of error
minimization. Figure 10 visualizes model inference on a single noisy sample
(bubbles), where GT measurements are denoted by solid blue line and models
estimates $\hat{\boldsymbol{x}}$ are marked with dashed lines.

Figure 10: Reconstruction comparison: Simulation

Since this dataset is simulated, different tradeoffs between deterministic and
stochastic noise sources were tested to examine models robustness. In the
absence of a calibration process or knowledge about inherent bias, SP-based
performances remained bounded by the sensors offset errors. In contrast, the
learning-based estimators are optimized during training phase, capable of
compensating wide range of error sources.

Table V summarizes the models denoising performances across the simulated
training set, providing dissimilarity measures between reconstructed signals
and the GT observations.

TABLE V: Reconstruction comparison: Simulation Model | RMSE [m/s2] | MAE [m/s2] | PSNR [dB] | RAE [$\%$]
---|---|---|---|---
kNN | 0.02481 | 0.01912 | 52.10719 | 0.44484
GRU | 0.02635 | 0.01950 | 51.39078 | 0.46668
LSTM | 0.02701 | 0.01999 | 50.10601 | 0.47835
RNN | 0.02898 | 0.02145 | 46.25170 | 0.51335
DWT | 0.06278 | 0.04623 | 43.81552 | 1.11545
SG | 0.06298 | 0.04678 | 43.62187 | 1.11605
MA | 0.06344 | 0.04872 | 43.42730 | 1.12842
Noisy | 0.06614 | 0.04913 | 41.73747 | 1.18383

Models are ranked according to RMSE results, where learning-based denoisers
have a clear advantage over SP-based algorithms. Since the kNN model leads the
table, let us define a suppression ratio $\gamma$ to quantify its ability to
minimize the reconstruction error, with respect to raw noisy measurements

$\displaystyle\gamma_{\text{sim.}}=\frac{\text{RAE}_{\text{kNN}}}{\text{RAE}_{\text{Noisy}}}\approx
37.5\%$ (40)

As seen, the relative absolute error (RAE) of the kNN estimates, cuts
dissimilarity error in more than half when compared to noisy measurements.
Generalization error is minimized the closer the output reaches the GT, but
increases proportionally with respect to sudden noise bursts.
Next, the denoised signals in Figure 10 are computed through the analytic SCA
equations (7, 8), to examine whether the reconstructed signals improve the
roll and pitch accuracy.

Figure 11: SCA comparison: Simulation

Figure 11 demonstrates how the computed angles obtained from models outputs,
improve the SCA procedure, as computed roll and pitch angles lay closer to the
ground-truth. Since the kNN is leading the accuracy, the next definition is
used to assess angular error between computed angle and GT angle

$\displaystyle\varepsilon=\alpha-\alpha_{GT}\ \in\ \mathbb{R}$ (41)

Figure 12 demonstrates the angular errors with respect to GT, as kNN estimated
angles (brown) are compared with noisy raw measurements (blue), over the
entire simulated test-set.

Figure 12: SCA results of kNN: Simulation

To understand the meaning of these results, Table VI utilizes RMSE as a
performance indicator, to summarize the kNN contribution to lowering the
overall angular error.

TABLE VI: Overall improvement rate: Simulation | $\varepsilon_{\text{Noisy}}\ [^{\circ}]$ | $\varepsilon_{\text{kNN}}\ [^{\circ}]$ | ratio [$\%$]
---|---|---|---
RMSE($\phi$) | 0.18452 | 0.14895 | 80.695
RMSE($\theta$) | 0.18147 | 0.14167 | 78.068

To conclude, synthetic noises and biases were generated to simulate
contaminated measurements and to examine the models capabilities to remove
them. Empirical experiments showed that each architecture exhibited different
bias-variance tradeoffs. Some estimators performed better using shallow
architectures such that variance error (overfitting) was reduced. Others
improved when estimated parameters were more noise-sensitive, thus avoiding
underfitting due to sensors biases. Over the entire test-set, it was shown
that learning-based estimates reduce reconstruction error by up to 37.5$\%$,
followed by 20$\%$ accuracy improvement of the SCA procedure.

###  VI-B Experimental Assessment

Unlike the simulated dataset, error sources in the experimental setup are
sensor-specific, and tradeoff between them cannot be modified. Here, noisy
samples are given by a consumer-grade smartphone sensor, whose readings are
submerged in lower levels of stochastic noise but higher bias levels, as
smartphones can be prone to mechanical shocks that impair the orthogonality of
the sensor axes.
In contrast, the GT references are taken from a high-end sensor, whose
accurate measurements are characterized by significantly lower noise levels
[53]. In other words, the generalization task here is to find an approximation
function that mimics best the GT sensor.

Similarly to before, Figure 13 compares the denoising capabilities of each
model, as the optimization task is to minimize dissimilarity between the model
estimates and GT.

Figure 13: Reconstruction comparison: Experiment

When compared to the simulated dataset, noise regime (bubbles) exhibits more
bias but less stochastic noise, as fluctuations are smaller. Here as well, the
learning-based algorithms lay significantly closer to GT measurements,
outperforming the SP-based models which lay in the vicinity of the noisy
measurements. It is not unusual however, as conventional filters may excel at
spatial filtering, but in the absence of a calibration process, constant
biases remain unfiltered.
Table VII summarizes the models reconstruction errors when computed across the
entire experimental test-set.

TABLE VII: Reconstruction comparison: Experiment Model | RMSE [m/s2] | MAE [m/s2] | PSNR [dB] | RAE [$\%$]
---|---|---|---|---
kNN | 0.00969 | 0.00658 | 60.09720 | 0.171414
GRU | 0.01373 | 0.01016 | 57.07227 | 0.242834
LSTM | 0.01414 | 0.01082 | 56.82121 | 0.249957
RNN | 0.01505 | 0.01231 | 56.27788 | 0.266084
SG | 0.10317 | 0.09124 | 39.81290 | 1.727121
DWT | 0.10326 | 0.09126 | 39.55169 | 1.823734
MA | 0.10588 | 0.09396 | 39.52071 | 1.825294
Noisy | 0.10843 | 0.09427 | 39.53768 | 1.828244

Once again, the kNN model demonstrates superiority among all models, as its
relative absolute error reduces error by an order of magnitude, as given by
the suppression ratio

$\displaystyle\gamma_{exp.}=\frac{\text{RAE}_{\text{kNN}}}{\text{RAE}_{\text{Noisy}}}\approx
9.38\%$ (42)

The strong differences between $\gamma_{sim.}$ and $\gamma_{exp.}$ can be
explained by unknown error sources that are present in real-world devices, but
underestimated or mismodeled during the simulation phase, thus beneficial to
the experimental results.
Next, the reconstructed signals shown in Figure 13, are used to estimate the
roll and pitch angles of the SCA procedure, as shown in Figure 14.

Figure 14: SCA comparison: Experiment

The learning-based denoisers exhibit better performances than the SP-based, as
their outputs manage to minimize the distance error with respect to the GT
angles. Focusing again on the best estimator, Figure 15 shows angular errors
of computed roll and pitch angles, comparing noisy measurements (blue) with
kNN estimates (brown), across the experimental test-set.

Figure 15: SCA results of kNN: Experiment

Despite greater noise amplitude of the instruments, the kNN model manages to
reduce reconstruction errors drastically. Table VIII shows the overall
improvement rate of the computed angles, using RMSE performance indicators to
compare the kNN contribution with respect to noisy sensor readings.

TABLE VIII: Overall improvement rate: Experiment | $\varepsilon_{\text{Noisy}}\ [^{\circ}]$ | $\varepsilon_{\text{kNN}}\ [^{\circ}]$ | ratio [$\%$]
---|---|---|---
RMSE($\phi$) | 0.30042 | 0.03398 | 11.311
RMSE($\theta$) | 0.50251 | 0.06676 | 13.285

Similarly to the reconstruction results (42), the angular errors are also
reduced by one order of magnitude, confirming the contribution of learning-
based denoising techniques.
To conclude the experimental section, accelerometer readings from a
commercial-grade smartphone sensor were used as noisy measurements, whereas
accurate readings from an aligned high-end sensor were used as a GT reference.
The denoising capabilities of different models were investigated, followed by
examining their contribution to the SCA procedure. Similarly to the simulated
scenario, our proposed kNN algorithm happened to outperform all other models,
managing to reduce the angular errors by one order of magnitude. We explain
its superiority by the simplistic non-parametric approach, which relies on
semantic similarity between independent variables, here specific force
measurements, with respect to $k$ nearest data points. It is noteworthy that
all visualizations were deliberately focused on short time periods, to
emphasize local fluctuations and responses.

##  VII Conclusions

In this work we addressed the challenging and previously unresearched task of
denoising accelerometer signals using data-driven methods, over simulated and
experimental datasets. By implementing and adjusting a wide range of denoising
algorithms from two different approaches, we assessed their denoising
capabilities in terms of pure signal reconstruction, followed by improvement
of an SCA procedure. The results showed a clear advantage of the learning-
based methods over conventional signal processing algorithms, presumably due
to their generalization ability to compensate a wide range of error sources.
Improvement rates seemed to vary between both datasets, depending on the
trade-off between noise components, as instruments exhibited a more intensive
noise regime. While errors in the simulated assessment dropped to half, in the
experimental assessment they were reduced by more than one order of magnitude.
That is to say that denoising accelerometers by learning-based approaches, can
be highly fruitful not only for pure noise reduction, but especially for
subsequent INS-related tasks.

## References

  * [1] K. R. Britting, _Inertial Navigation Systems Analysis_. John Wiley & Sons, Inc., 1971.
  * [2] C. Jekeli, _Inertial navigation systems with geodetic applications_. Walter de Gruyter Berlin, Germany, 2000.
  * [3] D. Titterton and J. L. Weston, _Strapdown Inertial Navigation Technology_. American Institute of Aeronautics and Astronautics and the Institution of Electrical Engineers, 2004\.
  * [4] P. D. Groves, _Principles of GNSS, Inertial and Multisensor Integrated Navigation Systems_. Artech House, 2013\.
  * [5] M. Dadafshar, “Accelerometer and gyroscopes sensors: operation, sensing, and applications,” _Maxim Integrated [online]_ , 2014.
  * [6] R. Gonzalez and C. A. Catania, “A statistical approach for optimal order adjustment of a moving average filter,” in _2018 IEEE/ION Position, Location and Navigation Symposium (PLANS)_. IEEE, 2018, pp. 1542–1546.
  * [7] A. Waegli, J. Skaloud, S. Guerrier, M. E. Parés, and I. Colomina, “Noise reduction and estimation in multiple micro-electro-mechanical inertial systems,” _Measurement Science and Technology_ , vol. 21, no. 6, p. 065201, 2010.
  * [8] Z. Diao, H. Quan, L. Lan, and Y. Han, “Analysis and compensation of mems gyroscope drift,” in _2013 Seventh International Conference on Sensing Technology (ICST)_. IEEE, 2013, pp. 592–596.
  * [9] S. Yong, C. Jiabin, S. Chunlei, and H. Yongqiang, “Research on the compensation in mems gyroscope random drift based on time-series analysis and kalman filtering,” in _2015 34th Chinese Control Conference (CCC)_. IEEE, 2015, pp. 2078–2082.
  * [10] Y.-H. Tu and C.-C. Peng, “An arma-based digital twin for mems gyroscope drift dynamics modeling and real-time compensation,” _IEEE Sensors Journal_ , vol. 21, no. 3, pp. 2712–2724, 2020.
  * [11] M. Ding, Z. Shi, B. Du, H. Wang, L. Han, and Z. Cheng, “The method of mems gyroscope random error compensation based on arma,” _Measurement Science and Technology_ , vol. 32, no. 12, p. 125109, 2021.
  * [12] S. Nassar, K.-P. Schwarz, A. Noureldin, _et al._ , “Modeling inertial sensor errors using autoregressive (ar) models,” _NAVIGATION, Journal of the Institute of Navigation_ , vol. 51, no. 4, pp. 259–268, 2004.
  * [13] S. Nassar, “Accurate insidgps positioning using ins data de-noising and autoregressive (ar) modeling of inertial sensor errors,” _Geomatica_ , vol. 59, no. 3, pp. 283–294, 2005.
  * [14] X. Wang and L. Wang, “Time-serial modeling and kalman filter of mems gyroscope random drift compensation,” in _2018 IEEE CSAA Guidance, Navigation and Control Conference (CGNCC)_. IEEE, 2018, pp. 1–5.
  * [15] G. Yuan, H. Liang, K. He, and Y. Xie, “Research on signal de-noising technique for mems gyro,” in _2010 3rd International Symposium on Systems and Control in Aeronautics and Astronautics_. IEEE, 2010, pp. 1288–1291.
  * [16] Y. Gan, L. Sui, J. Wu, B. Wang, Q. Zhang, and G. Xiao, “An emd threshold de-noising method for inertial sensors,” _Measurement_ , vol. 49, pp. 34–41, 2014.
  * [17] C. Liu, Z. Yang, Z. Shi, J. Ma, and J. Cao, “A gyroscope signal denoising method based on empirical mode decomposition and signal reconstruction,” _Sensors_ , vol. 19, no. 23, p. 5064, 2019.
  * [18] X. Wang, H. Cao, Y. Jiao, T. Lou, G. Ding, H. Zhao, and X. Duan, “Research on novel denoising method of variational mode decomposition in mems gyroscope,” _Measurement Science Review_ , vol. 21, no. 1, pp. 19–24, 2021.
  * [19] C. Shen, J. Li, X. Zhang, Y. Shi, J. Tang, H. Cao, and J. Liu, “A noise reduction method for dual-mass micro-electromechanical gyroscopes based on sample entropy empirical mode decomposition and time-frequency peak filtering,” _Sensors_ , vol. 16, no. 6, p. 796, 2016.
  * [20] X. Guo, C. Sun, P. Wang, and L. Huang, “Hybrid methods for mems gyro signal noise reduction with fast convergence rate and small steady-state error,” _Sensors and Actuators A: Physical_ , vol. 269, pp. 145–159, 2018.
  * [21] Y. Liu, G. Chen, Z. Wei, J. Yang, and D. Xing, “Denoising method of mems gyroscope based on interval empirical mode decomposition,” _Mathematical Problems in Engineering_ , vol. 2020, 2020.
  * [22] T. L. Grigorie, I. R. Edu, and J. I. Corcau, “Fuzzy logic denoising of the miniaturized inertial sensors in redundant configurations,” in _Proceedings of the ITI 2011, 33rd International Conference on Information Technology Interfaces_. IEEE, 2011, pp. 521–526.
  * [23] T. L. Grigorie and R. M. Botez, “Positioning monitoring improvement in a horizontal plane ins by using fuzzy logic data fusion for denoising of inertial sensors in redundant clusters,” _International Journal of Fuzzy Systems And Advanced Applications_ , vol. 2, pp. 33–40, 2015.
  * [24] Q. Li, X. Chen, and W. Xu, “Noise reduction of accelerometer signal with singular value decomposition and savitzky-golay filter,” _JOURNAL OF INFORMATION &COMPUTATIONAL SCIENCE_, vol. 10, no. 15, pp. 4783–4793, 2013.
  * [25] K. Nirmal, A. Sreejith, J. Mathew, M. Sarpotdar, A. Suresh, A. Prakash, M. Safonova, and J. Murthy, “Noise modeling and analysis of an imu-based attitude sensor: improvement of performance by filtering and sensor fusion,” in _Advances in Optical and Mechanical Technologies for Telescopes and Instrumentation II_ , vol. 9912. International Society for Optics and Photonics, 2016, p. 99126W.
  * [26] M. Karaim, A. Noureldin, and T. B. Karamat, “Low-cost imu data denoising using savitzky-golay filters,” in _2019 International Conference on Communications, Signal Processing, and their Applications (ICCSPA)_. IEEE, 2019, pp. 1–5.
  * [27] J. He, C. Sun, and P. Wang, “Noise reduction for mems gyroscope signal: a novel method combining acmp with adaptive multiscale sg filter based on ama,” _Sensors_ , vol. 19, no. 20, p. 4382, 2019.
  * [28] C. W. Kang, C. H. Kang, and C. G. Park, “Wavelet denoising technique for improvement of the low cost mems-gps integrated system,” _Wiley New York_ , 2010.
  * [29] C.-H. Kang, S.-Y. Kim, and C.-G. Park, “Improvement of a low cost mems inertial-gps integrated system using wavelet denoising techniques,” _International Journal of Aeronautical and Space Sciences_ , vol. 12, no. 4, pp. 371–378, 2011.
  * [30] N. El-Sheimy, S. Nassar, and A. Noureldin, “Wavelet de-noising for imu alignment,” _IEEE Aerospace and Electronic Systems Magazine_ , vol. 19, no. 10, pp. 32–39, 2004.
  * [31] F. Liu, F. Liu, W. Wang, and B. Xu, “Mems gyro’s output signal de-noising based on wavelet analysis,” in _2007 International Conference on Mechatronics and Automation_. IEEE, 2007, pp. 1288–1293.
  * [32] Z.-p. Li, Q.-j. Fan, L.-m. Chang, and X.-h. Yang, “Improved wavelet threshold denoising method for mems gyroscope,” in _11th IEEE International Conference on Control & Automation (ICCA)_. IEEE, 2014, pp. 530–534.
  * [33] J. Song, Z. Shi, B. Du, L. Han, H. Wang, and Z. Wang, “Mems gyroscope wavelet de-noising method based on redundancy and sparse representation,” _Microelectronic Engineering_ , vol. 217, p. 111112, 2019.
  * [34] Y. S. Shi and Z. F. Gao, “Study on mems gyro signal de-noising based on improved wavelet threshold method,” in _Applied Mechanics and Materials_ , vol. 433. Trans Tech Publ, 2013, pp. 1558–1562.
  * [35] J. Yuan, Y. Yuan, F. Liu, Y. Pang, and J. Lin, “An improved noise reduction algorithm based on wavelet transformation for mems gyroscope,” _Frontiers of Optoelectronics_ , vol. 8, no. 4, pp. 413–418, 2015.
  * [36] A. S. El-Wakeel, A. Noureldin, H. S. Hassanein, and N. Zorba, “Utilization of wavelet packet sensor de-noising for accurate positioning in intelligent road services,” in _2018 14th International Wireless Communications & Mobile Computing Conference (IWCMC)_. IEEE, 2018, pp. 1231–1236.
  * [37] T. A. Ali and A. M. Hasan, “A wavelet-NARX model for SDINS/GPS integration system,” _IOP Conference Series: Materials Science and Engineering_ , vol. 1094, no. 1, p. 012066, feb 2021. [Online]. Available: https://doi.org/10.1088/1757-899x/1094/1/012066
  * [38] M. Ben, W. J. Wei, W. J. Tong, and Z. X. Mei, “Mems gyro denoising based on second generation wavelet transform,” in _2010 First International Conference on Pervasive Computing, Signal Processing and Applications_. IEEE, 2010, pp. 255–258.
  * [39] R. Gonzalez and C. A. Catania, “Time-delayed multiple linear regression for de-noising mems inertial sensors,” _Computers & Electrical Engineering_, vol. 76, pp. 1–12, 2019.
  * [40] M. Brossard, S. Bonnabel, and A. Barrau, “Denoising imu gyroscopes with deep learning for open-loop attitude estimation,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 4796–4803, 2020.
  * [41] C. Jiang, S. Chen, Y. Chen, B. Zhang, Z. Feng, H. Zhou, and Y. Bo, “A mems imu de-noising method using long short term memory recurrent neural networks (lstm-rnn),” _Sensors_ , vol. 18, no. 10, p. 3470, 2018.
  * [42] C. Jiang, S. Chen, Y. Chen, Y. Bo, L. Han, J. Guo, Z. Feng, and H. Zhou, “Performance analysis of a deep simple recurrent unit recurrent neural network (sru-rnn) in mems gyroscope de-noising,” _Sensors_ , vol. 18, no. 12, p. 4471, 2018.
  * [43] Z. Ruoyu, G. Shuang, and C. Xiaowen, “Modeling of mems gyro drift based on wavelet threshold denoising and improved elman neural network,” in _2019 14th IEEE International Conference on Electronic Measurement & Instruments (ICEMI)_. IEEE, 2019, pp. 1754–1761.
  * [44] Z. Zhu, Y. Bo, and C. Jiang, “A mems gyroscope noise suppressing method using neural architecture search neural network,” _Mathematical Problems in Engineering_ , vol. 2019, 2019.
  * [45] C. Zhu, S. Cai, Y. Yang, W. Xu, H. Shen, and H. Chu, “A combined method for mems gyroscope error compensation using a long short-term memory network and kalman filter in random vibration environments,” _Sensors_ , vol. 21, no. 4, p. 1181, 2021.
  * [46] S. Han, Z. Meng, X. Zhang, and Y. Yan, “Hybrid deep recurrent neural networks for noise reduction of mems-imu with static and dynamic conditions,” _Micromachines_ , vol. 12, no. 2, p. 214, 2021.
  * [47] C. Jiang, Y. Chen, S. Chen, Y. Bo, W. Li, W. Tian, and J. Guo, “A mixed deep recurrent neural network for mems gyroscope noise suppressing,” _Electronics_ , vol. 8, no. 2, p. 181, 2019.
  * [48] J. Schmidhuber, “Deep learning in neural networks: An overview,” _Neural networks_ , vol. 61, pp. 85–117, 2015.
  * [49] J. L. Elman, “Finding structure in time,” _Cognitive science_ , vol. 14, no. 2, pp. 179–211, 1990.
  * [50] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [51] M. Schuster and K. K. Paliwal, “Bidirectional recurrent neural networks,” _IEEE transactions on Signal Processing_ , vol. 45, no. 11, pp. 2673–2681, 1997.
  * [52] K. Cho, B. Van Merriënboer, D. Bahdanau, and Y. Bengio, “On the properties of neural machine translation: Encoder-decoder approaches,” _arXiv preprint arXiv:1409.1259_ , 2014.
  * [53] A. Shurin, A. Saraev, M. Yona, Y. Gutnik, S. Faber, A. Etzion, and I. Klein, “The autonomous platforms inertial dataset,” _IEEE Access_ , vol. 10, pp. 10 191–10 201, 2022.
  * [54] I. Labs, “Motion Reference Unit,” https://www.inertiallabs.com/mru-datasheet, 2021, [Online; accessed 19-July-2021].
  * [55] T. T. Um, F. M. Pfister, D. Pichler, S. Endo, M. Lang, S. Hirche, U. Fietzek, and D. Kulić, “Data augmentation of wearable sensor data for parkinson’s disease monitoring using convolutional neural networks,” in _Proceedings of the 19th ACM international conference on multimodal interaction_ , 2017, pp. 216–220.
  * [56] H. Ohashi, M. Al-Nasser, S. Ahmed, T. Akiyama, T. Sato, P. Nguyen, K. Nakamura, and A. Dengel, “Augmenting wearable sensor data with physical constraint for dnn-based human-action recognition,” in _ICML 2017 times series workshop_ , 2017, pp. 6–11.
  * [57] O. Steven Eyobu and D. S. Han, “Feature representation and data augmentation for human activity classification based on wearable imu sensor data using a deep lstm neural network,” _Sensors_ , vol. 18, no. 9, p. 2892, 2018.

|  Daniel Engelsman received B.Sc. in Mechanical Engineering from Afeka
College of Engineering (2017) and M.Eng. in Aerospace Engineering from the
Technion-Israel Institute of Technology (2020). He is currently working as an
algorithm researcher in the Autonomous Navigation and Sensor Fusion Laboratory
at the University of Haifa, Department of Marine Technologies. His research
interest is focused on the thin line between conventional filtering techniques
and learning-based approaches for improving performance of inertial navigation
systems.
---|---
|  Itzik Klein (Senior Member, IEEE) received the B.Sc. and M.Sc. degrees in
Aerospace Engineering from the Technion - Israel Institute of Technology,
Haifa, Israel, in 2004 and 2007, respectively, and a Ph.D. degree in Geo-
information Engineering from the Technion - Israel Institute of Technology, in
2011. He is currently an Assistant Professor, heading the Autonomous
Navigation and Sensor Fusion Lab, at the Hatter Department of Marine
Technologies, University of Haifa. His research interests include data driven
based navigation, novel inertial navigation architectures, autonomous
underwater vehicles, sensor fusion, and estimation theory.
---|---


# A Constructive Approach to Function Realization
by Neural Stochastic Differential Equations

Tanya Veeravalli and Maxim Raginsky This work was supported in part by the NSF
under award CCF-2106358 (“Analysis and Geometry of Neural Dynamical Systems”)
and in part by the Illinois Institute for Data Science and Dynamical Systems
(iDS2), an NSF HDR TRIPODS institute, under award CCF-1934986.The authors are
with the Department of Electrical and Computer Engineering and the Coordinated
Science Laboratory, University of Illinois, Urbana, IL, USA. Emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The problem of function approximation by neural dynamical systems has
typically been approached in a top-down manner: Any continuous function can be
approximated to an arbitrary accuracy by a sufficiently complex model with a
given architecture. This can lead to high-complexity controls which are
impractical in applications. In this paper, we take the opposite, constructive
approach: We impose various structural restrictions on system dynamics and
consequently characterize the class of functions that can be realized by such
a system. The systems are implemented as a cascade interconnection of a neural
stochastic differential equation (Neural SDE), a deterministic dynamical
system, and a readout map. Both probabilistic and geometric (Lie-theoretic)
methods are used to characterize the classes of functions realized by such
systems.

###### Index Terms:

neural dynamical systems; stochastic differential equations; geometric control
theory

##  I Introduction

There is an extensive literature on function approximation by neural nets. For
instance, a classical result due to Leshno et al. [1] states the following:
Let $\sigma\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}\to{\mathbb{R}}$
be a continuous function which is not a polynomial. Then for any continuous
$f\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}^{n}\to{\mathbb{R}}$, any
compact $K\subset{\mathbb{R}}^{n}$, and any $\varepsilon>0$, there exist $N$
tuples
$(c_{i},a_{i},b_{i})\in{\mathbb{R}}\times{\mathbb{R}}^{n}\times{\mathbb{R}}$
such that

$\displaystyle\sup_{x\in
K}\Bigg{|}f(x)-\sum^{N}_{i=1}c_{i}\sigma(a_{i}^{\hbox{\it\tiny
T}}x+b_{i})\Bigg{|}\leq\varepsilon$ (1)

(in fact, the universal approximation property holds only if $\sigma$ is not a
polynomial). The approximating function in (1) is an instance of a neural net
with one hidden layer, $N$ hidden units (or neurons), and activation function
$\sigma$. There are also universal approximation results for multilayer neural
nets [2], as well as depth separation results [3] which show (constructively)
that there exist functions that admit efficient representation by multilayer
neural nets, yet require exponentially more hidden units for representation by
“shallower” nets.

Recently, motivated in large part by the widespread use of deep learning,
there has been a lot of interest in continuous-time abstractions of neural
nets [4, 5, 6, 7], aptly termed neural ODEs. These are modeled by dynamical
systems of the form

$\displaystyle\dot{z}(t)$ $\displaystyle=f(z(t),w(t)),\qquad z(0)=\alpha(x)$
(2a) $\displaystyle y$ $\displaystyle=\beta(z(1))$ (2b)

where the $q$-dimensional input $x$ is mapped to an initial $n$-dimensional
state $z(0)$ by a read-in map $\alpha\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{q}\to{\mathbb{R}}^{n}$, the $p$-dimensional output $y$
is produced from the state $z(1)$ at time $t=1$ by a read-out map
$\beta\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{n}\to{\mathbb{R}}^{p}$, and $f(\cdot,w)$ is a family
of vector fields on ${\mathbb{R}}^{n}$ parametrized by $w\in{\mathbb{R}}^{m}$.
For example, we could have something like $w=(C,A,b)\in{\mathbb{R}}^{n\times
k}\times{\mathbb{R}}^{k\times n}\times{\mathbb{R}}^{k}$ and

$\displaystyle f(x,w)=C\boldsymbol{\sigma}(Ax+b),$

where $\boldsymbol{\sigma}$ is the diagonal map

$\displaystyle\boldsymbol{\sigma}(z)\mathrel{\mathop{\mathchar
58\relax}}=\big{(}\sigma(z_{1}),\dots,\sigma(z_{k})\big{)}^{\hbox{\it\tiny
T}},\,z\in{\mathbb{R}}^{k}.$ (3)

The problem of universal approximation is then to establish the conditions
under which any continuous function $f\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{q}\to{\mathbb{R}}^{p}$ could be approximated, to any
given accuracy $\varepsilon$, on a given compact set
$K\subset{\mathbb{R}}^{q}$ by a system of the form (2) for a suitable choice
of the control law $w\mathrel{\mathop{\mathchar
58\relax}}[0,1]\to{\mathbb{R}}^{m}$.

This problem has been addressed in a number of recent works [8, 9, 10, 11].
While the methods and techniques differ, the overall philosophy of the results
is top-down—under suitable structural assumptions on $\alpha$, $\beta$, and
$f(\cdot,w)$, any continuous function can be approximated to any desired
accuracy by choosing a suitable control $w(\cdot)$. The resulting controls,
however, tend to exhibit fairly high complexity unlikely to be tolerated in
practical applications. For example, they may be piecewise constant, but the
number of pieces will grow like $(1/\varepsilon)^{d}$, where $d$ depends
polynomially on $q$ [8, 9, 10]; other constructions make use of high-gain,
high-frequency controls [11], with similar complexity issues.

###  I-A A constructive approach

In this paper, we approach the problem of function realization by neural
dynamical systems from a constructive, bottom-up perspective. That is, instead
of fixing a priori the class of functions to be approximated, we impose
various structural restrictions on the system dynamics and then characterize
the class of functions that can be realized by the system. To fix ideas, let
us consider a cascade interconnection [12] of two dynamical systems of the
following form:

$\displaystyle\dot{w}(t)$ $\displaystyle=g(w(t),\theta)$ (4a)
$\displaystyle\dot{z}(t)$ $\displaystyle=f(z(t),w(t)),$ (4b)

where $g(\cdot,\theta)$ is a family of vector fields on ${\mathbb{R}}^{m}$
parametrized by $\theta\in{\mathbb{R}}^{k}$. This system differs from (2) in
one key respect: The trajectory $w(t)$ is not prescribed externally, but is
instead generated internally according to (4a). Given the read-in map $\alpha$
and the read-out map $\beta$, as in (2), we can realize different functions of
$x\in{\mathbb{R}}^{q}$ by tuning the parameter $\theta$ [and, possibly, the
initial condition $w(0)$].

Neural ODE models of this type have been considered in the literature. For
example, Choromanski et al. [13] analyzed the case when $w$ takes values in
the orthogonal group $O(n)$, the vector fields $g(\cdot,\theta)$ are of the
form

$\displaystyle g(w,\theta)=w\tilde{g}(w,\theta)$

for some finitely parametrized map
$\tilde{g}(\cdot,\theta)\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{n\times n}\to{\rm Skew}(n)$, where ${\rm Skew}(n)$ is
the set of all $n\times n$ skew-symmetric matrices with real entries, and

$\displaystyle f(z,w)=\boldsymbol{\sigma}(wz),$

where $\boldsymbol{\sigma}$ is the coordinatewise application of some scalar
nonlinearity $\sigma\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}\to{\mathbb{R}}$, as in (3). While Choromanski et al.
did not explicitly characterize the class of functions that can be realized by
this model, they showed that they enjoy certain stability properties that are
advantageous in applications.

We will use cascade models like (4) as a starting point, but will replace the
deterministic dynamics (4a) with a stochastic differential equation (SDE)
driven by a multidimensional Brownian motion. The rationale for using such
neural SDEs [14, 15, 16, 17] is twofold: First, by using Brownian motion
processes as generalized inputs [18], we will be able to generate all the
internal complexity we need while only varying the global parameters of the
SDE. Second, it provides us with a sampling-based mechanism for constructing
finite approximations [19]. Indeed, if $F(x;\omega)$ is a random function of
$x$, then the expected value $f(x)\mathrel{\mathop{\mathchar
58\relax}}={\mathbf{E}}_{\omega}[F(x;\omega)]$ can be closely approximated by
finite sums of the form

$\displaystyle\hat{f}_{N}(x)=\frac{1}{N}\sum^{N}_{i=1}F(x;\omega_{i}),$

where $\omega_{1},\dots,\omega_{N}$ are independent copies of $\omega$, with
the $L^{2}$ approximation error ${\mathbf{E}}[(f(x)-\hat{f}_{N}(x))^{2}]$
scaling as $O(1/N)$ under suitable assumptions on the second moment of
$F(x;\omega)$. (While finite approximation is not the focus of this paper, our
earlier work [20] contains a discussion of sampling in the context of neural
SDEs.)

###  I-B Previous results

In an earlier work [20], we have obtained some constructive realizability
results for a certain class of neural SDEs. Briefly, we have considered
$n$-dimensional Itô SDEs

$\displaystyle\operatorname{d\\!}X_{t}$
$\displaystyle=a(X_{t};\theta)\operatorname{d\\!}t+b(X_{t};\theta)\operatorname{d\\!}V_{t},\qquad
X_{0}=x$ (5)

with finitely parametrized drift and diffusion coefficients satisfying the
uniform ellipticity condition—the eigenvalues of the matrix
$b(x;\theta)b(x;\theta)^{\hbox{\it\tiny T}}$ are uniformly bounded away from
$0$. Let $v$ be a fixed (nonrandom) vector in ${\mathbb{R}}^{n}$. We say that,
for each choice of $\theta$, the above system realizes the function
$f(x;\theta)\mathrel{\mathop{\mathchar
58\relax}}={\mathbf{E}}[Y]={\mathbf{E}}[v^{\hbox{\it\tiny T}}X_{1}]$ of the
initial condition $x$. We have shown that $f(x;\theta)$ can be approximated by

$\displaystyle f(x;\theta)\approx
c_{1}\int_{{\mathbb{R}}^{n}}v^{\hbox{\it\tiny
T}}z\cdot\exp\big{(}-c_{2}I(x,z;\theta)\big{)}\operatorname{d\\!}z,$

where $c_{1},c_{2}$ are some constants and where $I(x,z;\theta)$ is the
minimum of the control energy
$\frac{1}{2}\int^{1}_{0}|v(t)|^{2}\operatorname{d\\!}t$ among all sufficiently
regular (e.g., $L^{2}$) controls $v\mathrel{\mathop{\mathchar
58\relax}}[0,1]\to{\mathbb{R}}^{n}$ that transfer the state of the
deterministic control-affine system

$\displaystyle\dot{x}(t)=a(x(t);\theta)+b(x(t);\theta)v(t)$ (6)

from $x(0)=x$ to $x(1)=z$. Thus, the problem of characterizing the class of
functions realized in this way reduces to analyzing deterministic minimum-
energy control problems, and some upper and lower bounds on $I(x,z;\theta)$
are given in [20]. However, the assumption of uniform ellipticity is rather
restrictive and, in particular, rules out the interesting cases when the
diffusion process in (5) is degenerate, yet the deterministic system (6) is
nevertheless completely controllable [21, 22].

##  II The basic model and its properties

We start by presenting a relatively simple model which will allow for a clean
analysis. It has the following form:

$\displaystyle\operatorname{d\\!}W_{t}B$
$\displaystyle=a(W_{t};\theta)\operatorname{d\\!}t+b(W_{t};\theta)\circ\operatorname{d\\!}V_{t}$
(7a) $\displaystyle\operatorname{d\\!}Z_{t}$
$\displaystyle=h(W_{t},t)Z_{t}\operatorname{d\\!}t$ (7b) $\displaystyle Y$
$\displaystyle=Z_{1}\sigma(W_{1}^{\hbox{\it\tiny T}}x)$ (7c)

Here, (7a) is a Stratonovich SDE for an $n$-dimensional diffusion process
$(W_{t})$ with a deterministic initial condition $W_{0}=w_{0}$. It is driven
by an $m$-dimensional Brownian motion $(V_{t})$, and its drift and diffusion
coefficients are parametrized by a $k$-dimensional vector of parameters
$\theta$. We assume throughout that the drift and the diffusion coefficients
are sufficiently well-behaved (e.g., Lipschitz) to guarantee existence and
uniqueness of strong solutions of the SDE. The process $(Z_{t})$ is scalar,
initialized at $Z_{0}=1$, and $h\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{n}\times[0,1]\to{\mathbb{R}}$ is a continuous
function. The real-valued output $Y$ is obtained by multiplying $Z_{1}$ by
$\sigma(W_{1}^{\hbox{\it\tiny T}}x)$.

######  Remark 1.

By adding an equation $\operatorname{d\\!}X_{t}=0$ initialized at $X_{0}=x$,
we can view the above system as a stochastic variant of (4) with internal
state $(W_{t},Z_{t},X_{t})$, read-in map

$\displaystyle\alpha(x)=(w_{0},1,x),$

and read-out map

$\displaystyle\beta(w,z,x)=z\sigma(w^{\hbox{\it\tiny T}}x).$

The output $Y$ of (7) is a random function of $x$. We take its expected value
as the function realized by this system:

$\displaystyle f(x;\theta)\mathrel{\mathop{\mathchar
58\relax}}={\mathbf{E}}[Z_{1}\sigma(W^{\hbox{\it\tiny T}}_{1}x)],$ (8)

where in the left-hand side we have explicitly indicated the dependence of
this function on the parameters $\theta$ of the SDE (7a). The class of
functions realized in this way can be characterized as follows:

######  Theorem 1.

The function (8) realized by the system (7) can be expressed as

$\displaystyle
f(x;\theta)={\mathbf{E}}\Bigg{[}\exp\Bigg{(}\int^{1}_{0}h(W_{t},t)\operatorname{d\\!}t\Bigg{)}\sigma(W^{\hbox{\it\tiny
T}}_{1}x)\Bigg{]}.$ (9)

Moreover, let $\mathcal{L}^{\theta}$ be the second-order linear differential
operator

$\displaystyle\mathcal{L}^{\theta}\varphi(w)\mathrel{\mathop{\mathchar
58\relax}}=\tilde{a}(w;\theta)^{\hbox{\it\tiny
T}}\nabla\varphi(w)+\frac{1}{2}{\rm
tr}\\{b(w;\theta)b(w;\theta)^{\hbox{\it\tiny T}}\nabla^{2}\varphi(w)\\}$

where

$\displaystyle\tilde{a}(w;\theta)=a(w;\theta)+\frac{1}{2}\sum^{m}_{i=1}\frac{\partial
b_{i}}{\partial w}(w;\theta)b_{i}(w;\theta)$

is the Itô-corrected drift, and where $b_{i}(w;\theta)$ is the $i$th column of
$b(w;\theta)$. Consider the PDE

$\displaystyle\frac{\partial u}{\partial
w}(w,t;x,\theta)=\mathcal{L}^{\theta}u(w,t;x,\theta)+h(w,t)u(w,t;x,\theta)$
(10)

with the initial condition $u(w,0;x,\theta)=\sigma(w^{\hbox{\it\tiny T}}x)$.
If the function $w\mapsto\sigma(w^{\hbox{\it\tiny T}}x)$ is in the domain of
$\mathcal{L}^{\theta}$, then

$\displaystyle f(x;\theta)=u(w_{0},1;x,\theta).$

###### Proof.

The equation (7b) for $Z_{t}$ can be solved for each trajectory $(W_{t})$:

$\displaystyle Z_{t}$
$\displaystyle=Z_{0}\exp\Bigg{(}\int^{t}_{0}h(W_{s},s)\operatorname{d\\!}s\Bigg{)}$
$\displaystyle=\exp\Bigg{(}\int^{t}_{0}h(W_{s},s)\operatorname{d\\!}s\Bigg{)}.$

Substituting this into the expression for $Y$ and taking expectations, we
obtain (9).

The operator $\mathcal{L}^{\theta}$ is the infinitesimal generator of the
diffusion process $(W_{t})$ governed by the SDE (7a). Since
$w\mapsto\sigma(w^{\hbox{\it\tiny T}}x)$ is in the domain of
$\mathcal{L}^{\theta}$, the function

$\displaystyle u(w,t;x,\theta)$
$\displaystyle={\mathbf{E}}\Bigg{[}\exp\Bigg{(}\int^{t}_{0}h(W_{s},s)\operatorname{d\\!}s\Bigg{)}\sigma(W^{\hbox{\it\tiny
T}}_{t}x)\Bigg{|}W_{0}=w\Bigg{]}$

is a solution of (10) subject to the initial condition
$u(w,0;x,\theta)=\sigma(w^{\hbox{\it\tiny T}}x)$ by the (converse of) the
Feynman–Kac theorem [23, Thm. 21.1]. ∎

A more elaborate set-up, in the spirit of [13], is as follows:

$\displaystyle\operatorname{d\\!}W_{t}$
$\displaystyle=a(W_{t};\theta)\operatorname{d\\!}t+b(W_{t};\theta)\circ\operatorname{d\\!}V_{t}$
(11a) $\displaystyle\operatorname{d\\!}Z_{t}$
$\displaystyle=\boldsymbol{\sigma}(W_{t}Z_{t})\operatorname{d\\!}t$ (11b)
$\displaystyle Y$ $\displaystyle=v^{\hbox{\it\tiny T}}Z_{1}.$ (11c)

The processes $(W_{t})$ and $(Z_{t})$ in Eqs. (11) and (11b) are now evolving
in ${\mathbb{R}}^{n\times n}$ and ${\mathbb{R}}^{n}$, respectively, and $v$ in
(11c) is a fixed vector in ${\mathbb{R}}^{n}$. With the initial condition
$W_{0}=w_{0}$ fixed, we introduce the read-in map $\alpha(x)=(w_{0},x)$ and
the read-out map $\beta(w,z)=v^{\hbox{\it\tiny T}}z$. The stochastic dynamics
in (11a) generates the matrices of weights $W_{t}$, which are used to control
the dynamics of internal activations $Z_{t}$ in (11b). The cascade form of the
system allows us to first generate the random trajectory $(W_{t})$ and then,
conditionally on that, solve the ODE (11b). As before, we say that the model
(11) realizes the function
$f(x;\theta)={\mathbf{E}}[Y]={\mathbf{E}}[v^{\hbox{\it\tiny T}}Z_{1}]$. In
contrast to (7), where we were able to characterize the class of realized
functions in a relatively clean manner, the determination of $(Z_{t})$, given
a realization $(W_{t})$, involves solving a time-inhomogeneous ODE

$\displaystyle\frac{\operatorname{d\\!}Z_{t}}{\operatorname{d\\!}t}=g_{W_{t}}(Z_{t}),\,Z_{0}=x$
(12)

where we have defined

$\displaystyle g_{W_{t}}(z)\mathrel{\mathop{\mathchar
58\relax}}=\boldsymbol{\sigma}(W_{t}z).$ (13)

We can compactly express the function realized by (11) using the chronological
exponential notation [24] as

$\displaystyle f(x;\theta)={\mathbf{E}}\Bigg{[}v^{\hbox{\it\tiny
T}}\overrightarrow{\exp}\left(\int^{1}_{0}g_{W_{t}}\operatorname{d\\!}t\right)x\Bigg{]},$

which is just shorthand for the dependence of the solution of (12) on the
initial condition $x$ and on the trajectory $(W_{t})$, but any further
analysis would involve asymptotic expansions in terms of iterated integrals
and Lie derivatives.

###  II-A The role of Lie theory

The above construction is quite general as it does not specify the form of the
drift and the diffusion coefficients $a(w;\theta)$ and $b(w;\theta)$ in (7a)
and in (11a). A priori we would expect that, in order to make the model (7)
sufficiently expressive (i.e., to realize a sufficiently rich class of
functions), we would need $a$ and $b$ to depend nonlinearly on both $w$ and
$\theta$. This intuition, however, is misleading since we can let the
functions $h$ and $\sigma$ do all the “nonlinear work” when the Lie algebra
generated by the vector fields
$a(\cdot;\theta),b_{1}(\cdot;\theta),\dots,b_{m}(\cdot;\theta)$ is finite-
dimensional.

Thus, let $\mathcal{A}$ be a finite-dimensional Lie algebra of vector fields
on ${\mathbb{R}}^{n}$ and let $g_{1},\dots,g_{d}$ be a basis of $\mathcal{A}$,
$d=\dim\mathcal{A}$. For $\theta=(\theta_{ij})_{i=0,\dots,m,j=1,\dots d}$,
take

$\displaystyle a(w;\theta)=\sum^{d}_{j=1}\theta_{0j}g_{j}(w)$

and

$\displaystyle b_{i}(w;\theta)=\sum^{d}_{j=1}\theta_{ij}g_{j}(w),\qquad
i=1,\dots,m.$

Since $\mathcal{A}$ is finite-dimensional, it is isomorphic to a subalgebra of
${\rm gl}(\ell,{\mathbb{R}})$ (the Lie algebra of $\ell\times\ell$ real
matrices with the commutator bracket $[A,B]=AB-BA$) for some finite $\ell$ by
Ado’s theorem [25]. Thus, without loss of generality (and by changing $n$ if
necessary) we may take $g_{i}(w)=G_{i}w$ for some
$G_{i}\in{\mathbb{R}}^{n\times n}$, so (7a) will take the form

$\displaystyle\operatorname{d\\!}W_{t}=A(\theta)W_{t}\operatorname{d\\!}t+\sum^{m}_{i=1}B_{i}(\theta)W_{t}\circ\operatorname{d\\!}V^{i}_{t},$
(14)

where $V^{1}_{t},\dots,V^{m}_{t}$ are $m$ independent scalar Brownian motion
processes, and where

$\displaystyle A(\theta)$ $\displaystyle=\sum^{d}_{j=1}\theta_{0j}G_{j},$
$\displaystyle B_{i}(\theta)$
$\displaystyle=\sum^{d}_{j=1}\theta_{ij}G_{j},\quad i=1,\dots,m.$

As an illustration of the expressive capabilities of this type of linear
parametrization, we can consider Brockett’s construction of a diffusion
process on the sphere [22]: Let $A,B_{1},\dots,B_{m}$ be $n\times n$ skew-
symmetric matrices and consider the $n$-dimensional Stratonovich SDE

$\displaystyle\operatorname{d\\!}W_{t}=AW_{t}\operatorname{d\\!}t+\sum^{m}_{i=1}B_{i}W_{t}\circ\operatorname{d\\!}V^{i}_{t}.$
(15)

or the equivalent Itô SDE

$\displaystyle\operatorname{d\\!}W_{t}=\Bigg{(}A+\frac{1}{2}\sum^{m}_{i=1}B^{2}_{i}\Bigg{)}W_{t}\operatorname{d\\!}t+\sum^{m}_{i=1}B_{i}W_{t}\operatorname{d\\!}V^{i}_{t}.$

Applying Itô’s rule gives

$\displaystyle\operatorname{d\\!}\,(W^{\hbox{\it\tiny T}}_{t}W_{t})$
$\displaystyle=W^{\hbox{\it\tiny T}}_{t}\Bigg{(}A+A^{\hbox{\it\tiny
T}}+\sum^{m}_{i=1}B^{2}_{i}\Bigg{)}W_{\hbox{\it\tiny T}}\operatorname{d\\!}t$
$\displaystyle\qquad+\sum^{m}_{i=1}W^{\hbox{\it\tiny
T}}_{t}(B_{i}+B^{\hbox{\it\tiny T}}_{i})W_{t}\operatorname{d\\!}V^{i}_{t}$
$\displaystyle\qquad+\sum^{m}_{i=1}B^{\hbox{\it\tiny
T}}_{i}B_{i}\operatorname{d\\!}t$ $\displaystyle=0$

so the Euclidean norm $|W_{t}|$ stays constant for all $t$. Thus, if we choose
$w_{0}\in S^{n-1}$, then the random trajectory $(W_{t})_{t\geq 0}$ will be
confined to $S^{n-1}$. The algebra ${\rm Skew}(n)$ of $n\times n$ skew-
symmetric matrices has dimension $d=n(n-1)/2$, so we can express (15) in the
form (14) with a vector $\theta$ of $k=(m+1)n(n-1)/2$ parameters.

Similar considerations apply to the matrix-valued case. For example, if
$A,B_{1},\dots,B_{m}$ are $n\times n$ skew-symmetric matrices and the matrix
Stratonovich SDE

$\displaystyle\operatorname{d\\!}W_{t}=AW_{t}\operatorname{d\\!}t+\sum^{m}_{i=1}B_{i}W_{t}\circ\operatorname{d\\!}V^{i}_{t},$
(16)

is initialized with $W_{0}=w_{0}\in O(n)$, then all $W_{t}$ will take values
in $O(n)$ as well [22].

With regards to the role of $\sigma$, let us consider the family $\mathcal{F}$
of the vector fields $g_{W}(z)$, defined in (13). With typical choices of the
nonlinearity $\sigma\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}\to{\mathbb{R}}$ (e.g., the hyperbolic tangent
$\sigma(r)=\tanh r$), the Lie algebra generated by $\mathcal{F}$ will be
infinite-dimensional. To see this, suppose that $\sigma$ is $C^{\infty}$. Let
$g,g^{\prime}$ denote the vector fields $g_{W},g_{W^{\prime}}$ for two
matrices $W,W^{\prime}\in{\mathbb{R}}^{n}$. Then, using the formula

$\displaystyle\frac{\partial}{\partial
z}\boldsymbol{\sigma}(Wz)=\boldsymbol{\sigma}^{(1)}(Wz)W$

for the Jacobian of $g_{W}$, where

$\boldsymbol{\sigma}^{(1)}(z)\mathrel{\mathop{\mathchar 58\relax}}={\rm
diag}\big{(}\sigma^{\prime}(z_{1}),\dots,\sigma^{\prime}(z_{n})\big{)},$

a straightforward but tedious computation of the iterated Lie brackets

$\displaystyle{\rm ad}^{k+1}_{g}g^{\prime}={\rm ad}_{g}({\rm
ad}^{k}_{g}g^{\prime}),\quad k=0,1,\dots$

where ${\rm ad}^{0}_{g}g^{\prime}\mathrel{\mathop{\mathchar
58\relax}}=g^{\prime}$ and ${\rm ad}_{g}g^{\prime}\mathrel{\mathop{\mathchar
58\relax}}=[g,g^{\prime}],$ yields vector fields whose coordinates involve
derivatives of $\sigma$ of arbitrary orders. As an example, consider the case
$\sigma(r)=\tanh r$. Since the derivative of $\sigma$ satisfies the relation
$\sigma^{\prime}(r)=1-\sigma^{2}(r)$, the $k$th iterated Lie bracket ${\rm
ad}^{k}_{g}g^{\prime}$ will involve polynomials of degree $k+1$ in the entries
of $\boldsymbol{\sigma}(Wz)$ and $\boldsymbol{\sigma}(W^{\prime}z)$. Hence,
using the linear independence of the monomials $1,r,r^{2},\dots$, it is easy
to see that the Lie algebra generated by $\mathcal{F}$ will be infinite-
dimensional. Moreover, the Lie algebra generated by $\mathcal{F}$ may be
infinite-dimensional even if $\sigma$ is a polynomial. Consider, for example,
$\sigma(r)=1+r^{3}$; then even for $n=1$ the Lie algebra generated by
$\mathcal{F}$ will contain polynomial vector fields of arbitrary degree (and
thus any continuous function can be approximated arbitrarily well on any given
compact set by some element of this Lie algebra [11]). This is in sharp
contrast with the result of Leshno et al. [1] on the necessity of
nonpolynomial continuous activation functions for universal approximation of
continuous functions by neural nets with one hidden layer.

###  II-B Elaborations and extensions

The basic model in (7) can be extended in various ways. For instance, we can
consider matrix-valued processes $W_{t}$ and replace the output equation (7c)
with

$\displaystyle Y=Z_{1}v^{\hbox{\it\tiny T}}\boldsymbol{\sigma}(W_{1}x),$

where $W_{t}$ takes values in ${\mathbb{R}}^{p\times n}$ and
$v\in{\mathbb{R}}^{p}$ is a fixed vector. Another possibility is to let
$W_{t}=(U_{t},\tilde{W}_{t})$ take values in
${\mathbb{R}}^{p}\times{\mathbb{R}}^{p\times n}$ and let

$\displaystyle Y=Z_{1}U_{1}^{\hbox{\it\tiny
T}}\boldsymbol{\sigma}(\tilde{W}_{1}x).$

If $h\equiv 0$ in (7b), then the class of functions realized in this way is
the closed convex hull of neural nets with one hidden layer and with
arbitrarily many hidden units [26]. To see this, let
$(U^{\nu},\tilde{W}^{\nu})$, $\nu=1,\dots,N$, be independent copies of
$(U_{1},\tilde{W}_{1})$. (In other words, we just run $N$ independent copies
of the stochastic dynamics (7a) in parallel.) Then, assuming $U_{t}$ evolves
on the unit sphere $S^{n-1}$, $\tilde{W}_{t}$ evolves on the orthogonal group
$O(n)$, and the nonlinearity $\sigma\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}\to{\mathbb{R}}$ is bounded, the function

$\displaystyle\hat{f}_{N}(x)\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{N}\sum^{N}_{\nu=1}(U^{\nu})^{\hbox{\it\tiny
T}}\boldsymbol{\sigma}(W^{\nu}x)$

will be a good approximation of $f(x)={\mathbf{E}}[U^{\hbox{\it\tiny
T}}_{1}\boldsymbol{\sigma}(W_{1}x)]$, for $N$ sufficiently large, by the law
of large numbers. Fig. 1 shows a simple illustration of the ability of such
models to realize neuron-like functions of the form
$f(x)=\sigma(w^{\hbox{\it\tiny T}}x)$.

Figure 1: Left: Sample values of the single-neuron function
$f(x)=\tanh(w^{\hbox{\it\tiny T}}_{0}x)$ on the unit sphere $S^{2}$, where
$w_{0}\in S^{2}$ is a fixed unit vector. Right: The values of $f(x)$ predicted
by a model of the form (7) trained on random samples of $(X,f(X))$ with $X$
drawn from the uniform distribution on $S^{2}$. The weights $W_{t}$ evolve on
$S^{2}$ according to (15) with $m=2$.

The role of the function $h$ in (7b) is to bias or regularize the trajectories
$(W_{t})$ in some manner, as far as their contribution to the value of
$f(x;\theta)$ goes. For example, if some reference trajectory
$(\xi_{t})_{0\leq t\leq 1}$ is given, we could take

$\displaystyle h(w,t)=-|w-\xi_{t}|^{2}.$

In that case, we have

$\displaystyle f(x;\theta)={\mathbf{E}}\Bigg{[}\sigma(W^{\hbox{\it\tiny
T}}_{1}x)\exp\Bigg{(}-\int^{1}_{0}|W_{t}-\xi_{t}|^{2}\operatorname{d\\!}t\Bigg{)}\Bigg{]}$

which has the effect of penalizing those trajectories $(W_{t})$ that differ
too much from the reference.

The class of functions realized by the model (11) will be in general richer
than that realized by (7). By way of illustration, consider the SDE (16) for a
diffusion process on the orthogonal group $O(n)$. Suppose that the skew-
symmetric matrices $A,B_{1},\dots,B_{m}$ are such that the deterministic
system

$\displaystyle\dot{w}(t)=\Bigg{(}A+\sum^{m}_{i=1}v_{i}(t)B_{i}\Bigg{)}w(t)$
(17)

is controllable on $O(n)$ [22], i.e., any two matrices $w_{0},w_{1}\in O(n)$
can be joined by a curve lying along the trajectory $w(t)$ of (17) generated
by some piecewise constant controls $v_{1}(t),\dots,v_{m}(t)$. (For this, it
suffices to consider the case $w_{0}=I_{n}$ [27].) Now, the probability law of
the random path $(W_{t})_{t\geq 0}$ in $O(n)$ starting at $W_{0}=w_{0}$ is
concentrated on the closure of the set of all trajectories of (17) with
$w(0)=w_{0}$ generated by piecewise constant controls [21], which by
controllability covers the entire $O(n)$. This, in turn, will allow us to
sample high-complexity trajectories for controlling (12).

Finally, some comments on the generality of the considered models are in
order. Both (7) and (11) have a cascade structure, where the stochastic
dynamics of the weights $W_{t}$ is autonomous (for each choice of the
parameters $\theta$), and then these weights are used as generalized inputs
[18] to the controlled dynamics of the internal activations $Z_{t}$. Moreover,
the parameters $\theta$ explicitly enter only the dynamics of $W_{t}$.
However, we can consider other models of neural SDEs discussed in the
literature, for example a stochastically excited continuous-time recurrent
neural net

$\displaystyle\operatorname{d\\!}X_{t}$
$\displaystyle=\big{(}-X_{t}+\boldsymbol{\sigma}(X_{t})\big{)}\operatorname{d\\!}t+\sum^{m}_{i=1}B_{i}X_{t}\circ\operatorname{d\\!}V^{i}_{t}$

parametrized by $m$ matrices $B_{1},\dots,B_{m}$. In models of this type,
there is no explicit separation of the internal state into weights and
activations. However, we can use the results of Freedman and Willems [28] and
Krener and Lobry [29] to show that such neural SDEs can be simulated, in a
certain sense, by systems of the cascade type.

To that end, consider the $n$-dimensional stochastic system

$\displaystyle\operatorname{d\\!}X_{t}$
$\displaystyle=f(X_{t})\operatorname{d\\!}t+g(X_{t};\theta)\circ\operatorname{d\\!}V_{t}$
(18a) $\displaystyle Y_{t}$ $\displaystyle=v^{\hbox{\it\tiny T}}X_{t}$ (18b)

with initial condition $X_{0}=x$, where the columns of the $n\times m$ matrix
$g(x;\theta)$ are smooth vector fields on ${\mathbb{R}}^{n}$ with a smooth
dependence on a $k$-dimensional parameter vector $\theta$. Then we have the
following result:

######  Theorem 2.

Suppose that the Lie algebra $\mathcal{A}$ generated by
$g_{1}(\cdot;\theta),\dots,g_{m}(\cdot;\theta)$, $\theta\in{\mathbb{R}}^{k}$,
has finite dimension $d$. Then there exist smooth vector fields
$b_{1},\dots,b_{d}$ on ${\mathbb{R}}^{d}$, smooth functions
$\beta_{ij}\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{k}\to{\mathbb{R}}$, $1\leq i\leq m,1\leq j\leq d$,
$h\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{n}\times{\mathbb{R}}^{d}\to{\mathbb{R}}^{n}$,
$\varphi\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{d}\times{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$, and an
almost surely positive stopping time $\tau$, such that the cascade system

$\displaystyle\operatorname{d\\!}W_{t}$
$\displaystyle=\sum^{m}_{i=1}\Bigg{(}\sum^{d}_{j=1}\beta_{ij}(\theta)b_{j}(W_{t})\Bigg{)}\circ\operatorname{d\\!}V^{i}_{t}$
(19a) $\displaystyle\operatorname{d\\!}Z_{t}$
$\displaystyle=h(Z_{t},W_{t})\operatorname{d\\!}t$ (19b)

with $W_{0}=0$, $Z_{0}=x$ simulates (18) up to time $\tau$, i.e.,

$\displaystyle Y_{t}=v^{\hbox{\it\tiny T}}\varphi(W_{t},Z_{t}),\qquad 0\leq
t<\tau.$

######  Remark 2.

The above theorem shows that we can simulate a system like (18), at least
locally, by a cascade system like (19) by means of a smooth reparametrization
$\theta\mapsto(\beta_{ij}(\theta))$ and a smooth “decoding map”
$\varphi(w,z)=x$ that converts the internal weights $w\in{\mathbb{R}}^{d}$ and
internal activations $z\in{\mathbb{R}}^{n}$ into the overall state
$x\in{\mathbb{R}}^{n}$. Observe that the transformed parameters
$\beta_{ij}(\theta)$ enter into (19a) linearly and that the vector fields
$b_{1},\dots,b_{d}$ and the maps $h$ and $\varphi$ do not have any dependence
$\theta$. Moreover, since we have not imposed any restrictions on the drift
vector field $f$ (apart from the regularity conditions needed to ensure
existence and uniqueness of Stratonovich solutions), the Lie algebra generated
by the vector fields $h(\cdot,w)$ may be infinite-dimensional.

###### Proof.

We closely follow the proof of Theorem 1 in [29] (the special case of Abelian
$\mathcal{A}$ was considered in an earlier paper by Freedman and Willems
[28]). Let $\tilde{g}_{1},\ldots,\tilde{g}_{d}$ be a basis of $\mathcal{A}$,
so that

$\displaystyle
g_{i}(x;\theta)=\sum^{d}_{j=1}\beta_{ij}(\theta)\tilde{g}_{j}(x),\quad
i=1,\dots,m$ (20)

for some smooth maps $\beta_{ij}\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{k}\to{\mathbb{R}}$. Let

$\displaystyle\varphi(w,z)\mathrel{\mathop{\mathchar
58\relax}}=e^{w_{1}\tilde{g}_{1}}\circ e^{w_{2}\tilde{g}_{2}}\circ\dots\circ
e^{w_{d}\tilde{g}_{d}}z,$

where $e^{t\tilde{g}_{i}}$ denotes the flow map of $\tilde{g}_{i}$, i.e.,
$e^{t\tilde{g}_{i}}z$ is the solution, at time $t$, of the ODE

$\displaystyle\dot{z}(t)=\tilde{g}_{i}(z(t)),\qquad z(0)=z.$

We have $\varphi(0,x)=x$. Using the methods of [29], we can show that there
exist smooth vector fields $b_{1},\dots,b_{d}$ on ${\mathbb{R}}^{d}$ and a
neighborhood $\mathcal{U}$ of $(0,x)$, such that

$\displaystyle\frac{\partial\varphi}{\partial
w}(w,z)b_{i}(w)=\tilde{g}_{i}(\varphi(w,z)),\qquad i=1,\dots,m$

for all $(w,z)\in\mathcal{U}$. Consequently,

$\displaystyle\frac{\partial\varphi}{\partial
w}(w,z)\bigg{(}\sum^{d}_{j=1}\beta_{ij}(\theta)b_{i}(w)\bigg{)}=g_{i}(\varphi(w,z);\theta)$

for all $1\leq i\leq d$ and all $\theta\in{\mathbb{R}}^{k}$. Now,
$\varphi(w,z)$ is equal to the solution, at $t=d$, of the time-inhomogeneous
ODE

$\displaystyle\dot{\xi}(t)=v(\xi(t),t),\qquad z(0)=z$

with

$\displaystyle v(z,t)\mathrel{\mathop{\mathchar
58\relax}}=w_{d-i+1}\tilde{g}_{d-i+1}(z),\qquad i-1\leq t<i$

and the Jacobian $\frac{\partial\varphi}{\partial z}(w,z)$ is equal to the
solution, at $t=d$, of the variational equation

$\displaystyle\dot{\Lambda}(t)=\frac{\partial v}{\partial
z}(\xi(t),t)\Lambda(t),\qquad\Lambda(0)=I_{n}.$

Hence, it is invertible, so we take

$\displaystyle h(z,w)\mathrel{\mathop{\mathchar
58\relax}}=\bigg{(}\frac{\partial\varphi}{\partial
z}(w,z)\bigg{)}^{-1}f(\varphi(w,z))$

for all $(w,z)\in\mathcal{U}$. Let $(W_{t},Z_{t})$ be the solution of (19)
with $(W_{0},Z_{0})=(0,x)$, and let $\tau$ be the first exit time from
$\mathcal{U}$:

$\displaystyle\tau\mathrel{\mathop{\mathchar
58\relax}}=\inf\\{t>0\mathrel{\mathop{\mathchar
58\relax}}(W_{t},Z_{t})\not\in\mathcal{U}\\}.$

Then, since Stratonovich differentials obey the rules of ordinary calculus,
for $0\leq t<\tau$ we have

$\displaystyle\operatorname{d\\!}X_{t}$
$\displaystyle=\operatorname{d\\!}\varphi(W_{t},Z_{t})$
$\displaystyle=\frac{\partial\varphi}{\partial
w}(W_{t},Z_{t})\circ\operatorname{d\\!}W_{t}+\frac{\partial\varphi}{\partial
z}(W_{t},Z_{t})\circ\operatorname{d\\!}Z_{t}$
$\displaystyle=f(\varphi(W_{t},Z_{t}))\operatorname{d\\!}t+\sum^{m}_{i=1}g_{i}(\varphi(W_{t},Z_{t});\theta)\circ\operatorname{d\\!}V^{i}_{t}$
$\displaystyle=f(X_{t})\operatorname{d\\!}t+\sum^{m}_{i=1}g_{i}(X_{t};\theta)\circ\operatorname{d\\!}V^{i}_{t}$

and $v^{\hbox{\it\tiny T}}\varphi(W_{t},Z_{t})=v^{\hbox{\it\tiny
T}}X_{t}=Y_{t}$. This gives the desired simulation property. ∎

## Acknowledgments

T. Veeravalli would like to thank A.J. Havens for a PyTorch crash course. M.
Raginsky would like to thank A. Belabbas for useful suggestions in the early
stages of this work.

## References

  * [1] M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken, “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function,” _Neural Networks_ , vol. 6, pp. 861–867, 1993.
  * [2] D. Yarotsky, “Error bounds for approximations with deep ReLU networks,” _Neural Networks_ , vol. 94, pp. 103–114, 2017.
  * [3] M. Telgarsky, “Benefits of depth in neural networks,” in _Conference on Learning Theory_ , 2016.
  * [4] M. W. Hirsch, “Convergent activation dynamics in continuous time networks,” _Neural Networks_ , vol. 2, pp. 331–349, 1989.
  * [5] E. Haber and L. Ruthotto, “Stable architectures for deep neural networks,” _Inverse Problems_ , vol. 34, no. 1, p. 014004, December 2017.
  * [6] W. E, “A proposal on machine learning via dynamical systems,” _Communications in Mathematics and Statistics_ , vol. 5, pp. 1–11, 2017.
  * [7] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural ordinary differential equations,” in _NeurIPS_ , 2018.
  * [8] Q. Li, T. Lin, and Z. Shen, “Deep learning via dynamical systems: An approximation perspective,” 2019. [Online]. Available: http:arxiv.org/abs/1912.10382
  * [9] D. Ruiz-Balet and E. Zuazua, “Neural ODE control for classification, approximation and transport,” 2021. [Online]. Available: https://arxiv.org/abs/2104.05278
  * [10] P. Tabuada and B. Gharesifard, “Universal approximation power of deep neural networks via nonlinear control theory,” in _ICLR_ , 2021.
  * [11] A. Agrachev and A. Sarychev, “Control on the manifolds of mappings with a view to the deep learning,” _Journal of Dynamical and Control Systems_ , vol. 28, no. 4, pp. 989–1008, 2022.
  * [12] A. J. Krener, “A decomposition theory for differentiable systems,” _SIAM Journal on Control and Optimization_ , vol. 15, no. 5, pp. 813–829, 1977.
  * [13] K. M. Choromanski, J. Q. Davis, V. Likhosherstov, X. Song, J.-J. Slotine, J. Varley, H. Lee, A. Weller, and V. Sindhwani, “Ode to an ODE,” in _Advances in Neural Information Processing Systems_ , H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, Eds., vol. 33, 2020, pp. 3338–3350.
  * [14] E. Wong, “Stochastic neural networks,” _Algorithmica_ , vol. 6, pp. 466–478, June 1991.
  * [15] B. Tzen and M. Raginsky, “Theoretical guarantees for sampling and inference in generative models with latent diffusions,” in _COLT_ , 2019.
  * [16] ——, “Neural stochastic differential equations: Deep latent gaussian models in the diffusion limit,” 2019. [Online]. Available: http://arxiv.org/abs/1905.09883
  * [17] X. Li, T.-K. L. Wong, R. T. Q. Chen, and D. Duvenaud, “Scalable gradients for stochastic differential equations,” in _AISTATS_ , 2020.
  * [18] H. Sussmann, “On generalized inputs and white noise,” in _IEEE Conference on Decision and Control_ , 1976, pp. 809–814.
  * [19] A. R. Barron, “Universal approximation bounds for superpositions of a sigmoidal function,” _IEEE Transactions on Information Theory_ , vol. 39, no. 3, pp. 930–945, 1993.
  * [20] T. Veeravalli and M. Raginsky, “Nonlinear controllability and function representation by neural stochastic differential equations,” in _Conference on Learning for Dynamics and Control (L4DC)_ , 2023.
  * [21] D. L. Elliott, “Diffusions on manifolds arising from controllable systems,” in _Geometric Methods in System Theory_ , D. Q. Mayne and R. W. Brockett, Eds. D. Reidel Publishing Co., 1973, pp. 285–294.
  * [22] R. W. Brockett, “Lie theory and control systems defined on spheres,” _SIAM Journal on Applied Mathematics_ , vol. 25, no. 2, pp. 213–225, 1973\.
  * [23] O. Kallenberg, _Foundations of Modern Porbability_. Springer, 1997.
  * [24] A. A. Agrachev and Y. L. Sachkov, _Control Theory from the Geometric Viewpoint_. Springer, 2004.
  * [25] V. S. Varadarajan, _Lie Groups, Lie Algebras, and Their Representations_. Springer, 1984.
  * [26] L. Gurvits and P. Koiran, “Approximation and learning of convex superpositions,” _Journal of Computer and System Sciences_ , vol. 55, pp. 161–170, 1997.
  * [27] V. Jurdjevic and H. J. Sussmann, “Control systems on Lie groups,” _Journal of Differential Equations_ , vol. 12, pp. 313–329, 1972.
  * [28] M. Freedman and J. Willems, “Smooth representation of systems with differentiated inputs,” _IEEE Transactions on Automatic Control_ , vol. 23, no. 1, pp. 16–21, 1978.
  * [29] A. J. Krener and C. Lobry, “The complexity of stochastic differential equations,” _Stochastics: An International Journal of Probability and Stochastic Processes_ , vol. 4, no. 3, pp. 193–203, 1981.