# Geometry of Radial Basis Neural Networks for
Safety Biased Approximation of Unsafe Regions

Ahmad Abuaish†, Mohit Srinivasan‡, Patricio A. Vela† This work was supported
in part by the National Science Foundation under Award S&AS #1849333, by DARPA
PAI, and by KACST Fellowship.† School of Electrical and Computer Engineering,
Georgia Institute of Technology, Atlanta, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>Ford Motor Company, Dearborn, USA<EMAIL_ADDRESS>
###### Abstract

Barrier function-based inequality constraints are a means to enforce safety
specifications for control systems. When used in conjunction with a convex
optimization program, they provide a computationally efficient method to
enforce safety for the general class of control-affine systems. One of the
main assumptions when taking this approach is the a priori knowledge of the
barrier function itself, i.e., knowledge of the safe set. In the context of
navigation through unknown environments where the locally safe set evolves
with time, such knowledge does not exist. This manuscript focuses on the
synthesis of a zeroing barrier function characterizing the safe set based on
safe and unsafe sample measurements, e.g., from perception data in navigation
applications. Prior work formulated a supervised machine learning algorithm
whose solution guaranteed the construction of a zeroing barrier function with
specific level-set properties. However, it did not explore the geometry of the
neural network design used for the synthesis process. This manuscript
describes the specific geometry of the neural network used for zeroing barrier
function synthesis, and shows how the network provides the necessary
representation for splitting the state space into safe and unsafe regions.

##  I Introduction

Invariant set-based safe control synthesis [1] has become a favorable
technique to enforce safety, as it provides theoretical guarantees when used
to augment base controllers [2]. For closed dynamical systems, the invariant
set is represented by the zero and positive level-sets of a continuously
differentiable implicit function, typically termed _zeroing barrier function_
(ZBF). Subsequently, for controlled systems, _control barrier functions_
(CBFs) are introduced to represent the invariant set. In safety-critical
applications, CBF-based controllers are designed to render safe regions in the
state space a positively invariant set. Prevailing use of CBFs is in a point-
wise optimization program solved via quadratic programming (QP) [3]. CBF-based
QPs have been used in a wide range of applications in robotics such as
collision avoidance [4, 5], multi-robot coordination and task satisfaction [6,
7, 8], and automotive applications [3].

Usually, a CBF is handcrafted based on a priori knowledge of the safe regions
in the state space. However, there are applications where the safe regions
evolve with time, with navigation being one. In these applications, it is
critical to synthesize the CBF online using sensor measurements.
Unfortunately, determining the true invariant region defined by a CBF is
generally an NP-hard problem, but some techniques exist that can estimate the
invariant region [9]. Since the existence of a zeroing barrier function is
needed to formally confirm the existence of a CBF, this manuscript solely
focuses on zeroing barrier function synthesis to separate state space into
_safe_ and _unsafe_ regions. The ZBF is constructed from a two-layer kernel
machine network trained from a labeled dataset of safe and unsafe samples.
Further, the geometry of the kernel functions in the network is explored to
efficiently partition the space. The approach is motivated by the
Kolmogorov–Arnold representation theorem, which implies that two-layer neural
networks may be capable of approximating continuous functions [10].

Recently, several data-driven approaches for constructing a CBF were proposed
to account for uncertainties in either the system dynamics, unsafe regions, or
both. One approach category for learning CBF’s involves supervised offline
learning. Instances include imitation learning where training data is
generated by an expert actor or optimal control simulations [11, 12]. The
offline nature lacks the ability to accommodate real-time changes in the
environment. In contrast, _self-supervised_ approaches permit online learning.
Initial work on self-supervised Bayesian learning system of uncertain dynamics
[13] with known barrier functions was merged with [14] to learn the system
dynamics and an implicit function representation of the unsafe region [15]. In
[14], a signed distance function representing obstacles is modeled as a deep
neural network trained from range sensor data via stochastic gradient descent
(SGD) with replay memory. Similarly, [16] presented the construction of a
probabilistic occupancy map from a kernel-based logistic regression model
trained from range sensor data via SGD. However, there were no hard
constraints on misclassifying unsafe data points in the underlying SGD
optimization process, which nullifies any possible theoretical safety
guarantee.

Our previous work focused on creating a ZBF for navigation applications based
on data collected from LiDAR sensors [17]. A two-layer network with Gaussian
radial basis functions (GRBFs) was synthesized from this data. The first layer
used sparsely distributed (over the domain) GRBF centers, while the second
GRBF layer was learnt during the kernel support vector machine (kSVM)
optimization process. The kSVM optimization specifications provided formal
guarantees regarding the partitioning of the domain into safe and unsafe
regions. However, the work did not discuss the geometry and associated
properties of the GRBF network. This work analyzes the geometry of the two-
layer network and the structure of the optimization problem to prove the
existence of a ZBF with known partitioning properties.

The manuscript organization is as follows: Section II discusses the geometry
of Gaussian kernel functions with respect to the first layer. Section III
covers the second layer construction, geometry, and optimization
specifications. Section IV discusses the qualification of the two-layer kernel
machine network to be a zeroing barrier function along with a kernel basis
selection strategy. Sections V and VI present case studies and concluding
remarks, respectively.

##  II Geometry of the Kernel Hilbert Space

GRBF Neural Networks (GRBF-NNs) have several properties ideal for zeroing
barrier function creation. First they are universal approximators [18], second
they exhibit locality [19], and third they partition space. Since the last
property is less frequently mentioned, but often used, this section devotes
attention to space partitioning, as it is essential for data-driven safe set
generation.

The focus of this section is on GRBF-NNs as kernel machines whose kernel
functions are radial basis functions (RBFs), which generally take the
following form,

${k}(x_{i},x_{j})=\varphi(\left|\left|x_{i}-x_{j}\right|\right|;\sigma),\quad\text{for}\quad\varphi:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+},$
(1)

where $x_{i},x_{j}\in\mathbb{R}^{n_{d}}$, $\left|\left|\cdot\right|\right|$ is
a norm and $\sigma$ is the scalar bandwidth, which influences the sensitivity
of the basis function $\varphi$. A kernel machine based on radial basis
functions generates a mapping through the use of multiple kernel function
mappings with fixed argument elements $c_{j}$ in a center set $\mathcal{C}$.
For $c_{j}\in\mathcal{C}$ the kernel mapping is

${\vec{k}}:\mathbb{R}^{n_{d}}\rightarrow\mathcal{H}^{n_{c}},\ \
x\mapsto[{k}(x,c_{1})\;\cdots\;{k}(x,c_{n_{c}})]^{T}.$ (2)

where $n_{c}=\left|{\mathcal{C}}\right|$ and $\mathcal{H}^{\cdot}$ indicates
that the output space is a Hilbert space. For this kernel machine to define a
scalar function requires specifying $\alpha\in\mathbb{R}^{n_{c}}$ such that

$f_{\text{KM}}(x)={\left<\,\alpha\;,\;{\vec{k}}(x)\,\right>}=\alpha^{T}{\vec{k}}(x).$
(3)

Kernel machines permit more general function classes than RBFs in
$\mathbb{R}^{n_{d}}$ (subsequent sections will use polynomial kernel
functions). That said, RBFs have geometric properties that are implicitly
exploited in kernel machine learning applications. Specializing to the case of
the Gaussian radial basis function (or Gaussian kernel), let the kernel
function be

$k_{G}(x,c)=\exp\left(-\frac{\left|\left|x-c\right|\right|^{2}}{\sigma^{2}}\right).$
(4)

######  Theorem 1.

The kernel mapping, with $n_{c}$ Gaussian kernels, maps the input domain
$\mathcal{D}\subset\mathbb{R}^{n_{d}}$ into a surface in the Hilbert space
$\mathcal{H}^{n_{c}}\subset\mathbb{R}^{n_{c}}$ when $n_{c}\geq n_{d}+1$ and
there are $n_{d}+1$ centers capable of defining a coordinate system in
$n_{d}$.

###### Proof.

Showing that the kernel mapping is 1-1 establishes this property. The pre-
image, ${k}^{-1}(\cdot,c_{i})$, of each coordinate’s kernel mapping is a
sphere in the input space. The intersection of all spheres for all coordinate
mappings establishes the pre-image point, which is unique only if the
intersection is unique. Finite solutions for the intersections has been proven
in the context of rigid body geometry for $n_{c}=n_{d}$ [20], with $n_{c}\geq
n_{d}+1$ necessary for a single valid solution. The requirement for the
centers is that they lead to a basis of $n_{d}$ vectors when using one of the
points as the origin and using $n_{d}$ other points to obtain the basis
vectors relative to that origin. Each pre-image imposes a constraint on the
degrees of freedom of the input point $x\in\mathbb{R}^{n_{d}}$, such that the
$n_{d}+1$ intersecting pre-image spheres do so at a single point. When
$n_{c}>n_{d}+1$, the additional pre-image constraints are redundant and
effectively impose no constraints. The rank of the mapping is $n_{d}<n_{c}$,
hence it maps to a surface of dimension $n_{d}$. ∎

Theorem 1 applies to any RBF with infinite support; an RBF with finite support
will have similar properties but will include an $\epsilon$-covering
constraint. Basis functions generated from other norms also have similar
properties but may require more centers to induce a 1-1 mapping. If the kernel
is differentiable, then the 1-1 mapping is an embedding. Safe set generation
with an appropriate kernel mapping will involve defining the concept of a
kernel embedding111Not to be confused with the kernel mean embedding which is
for probability distributions [21]. and its associated kernel embedding
inducing data set.

######  Definition 1.

The kernel mapping
${\vec{k}}:\mathbb{R}^{n_{d}}\rightarrow\mathcal{H}^{n_{c}}$ is a kernel
embedding if the center set $\mathcal{C}\subset\mathcal{D}$ generating the
kernel mapping are such that it is 1-1 and the kernel function
${k}:\mathbb{R}^{n_{d}}\times\mathbb{R}^{n_{d}}\rightarrow\mathbb{R}$ is
differentiable. The set of centers is called the kernel embedding inducing set
(KEI set).

######  Corollary 1.

Consider a finite set of points $\mathcal{X}_{p}\subset\mathbb{R}^{2}$ with an
associated triangulation. Under a kernel embedding, each triangular region
maps to a surface in $\mathcal{H}^{n_{c}}$ homeomorphic to a 2-simplex.

######  Corollary 2.

Consider a finite set of points $\mathcal{X}_{p}\subset\mathbb{R}^{3}$ with an
associated tetrahedralization. Under a kernel embedding, each tetrahedron maps
to a surface in $\mathcal{H}^{n_{c}}$ homeomorphic to a 3-simplex.

If the domain $\mathcal{D}\subset\mathbb{R}^{n_{d}}$ is a set of disconnected
regions excluding points at infinity, then a kernel embedding will map to a
set of disconnected surfaces. Similarly, a collection of non-intersecting
triangulations/tetrahedralizations maps to a set of disconnected surfaces
under a kernel embedding.

####  II-1 Geometry of a Kernel Embedding

The Gaussian kernel mapping outputs lie in the unit cube of
$\mathcal{H}^{n_{c}}\subset\mathbb{R}^{n_{c}}$. Each center $c_{i}$ in the set
$\mathcal{C}$ maximizes its associated coordinate (evaluates to 1), which
means that there is a neighborhood of $c_{i}$ in $\mathcal{D}$ for which this
same coordinate is also maximal for all points in the neighborhood. Points
tending to infinity map to the origin in $\mathcal{H}^{n_{c}}$ since the
Gaussian radial basis function tends to zero as the input radius tends to
infinity.

######  Corollary 3.

For a kernel embedding defined using the Gaussian kernel, a compact domain
$\mathcal{D}$ maps to a compact surface whose points lie outside of a ball
centered at the origin. Furthermore,
${\vec{k}}(\mathbb{R}^{n_{d}})\cup\\{\vec{0}\\}$ is a compact surface.

With the surface geometry of the kernel embedding ${\vec{k}}$ for a given KEI
set $\mathcal{C}$ established, the next step is to consider partitions of the
surface in $\mathcal{H}^{n_{c}}$. Given that the unbounded input space maps to
a compact surface, defining a partition of the space is equivalent to defining
a cutting surface transverse to the compact surface that divides space into
positive and negative regions. The intersection of the compact surface with
the cutting surface in the Hilbert space maps back to separating boundaries in
the original space. The objective will be to define cutting surfaces by
specific level-sets of an implicit function in the Hilbert space. In short,
safe region generation in the input space involves construction of an implicit
function in the Hilbert space based on sampled points with safe / unsafe
labels.

Figure 1: Depiction of the input domain and the associated Hilbert space along
with the GRBF centers and the separating boundaries and surfaces.

The observation leading to Corollary 3 hints that the geometry induced by the
KEI set plays a role in establishing the surface image in the Hilbert space.
Here, we connect the KEI set to coverings of the input domain. For the sets
and spaces defined here, define an $\epsilon$-covering as follows:

######  Definition 1.

For $\mathbb{R}^{n_{d}}$ with the Euclidean norm (2-norm), let
$\mathcal{D}\subset\mathbb{R}^{n_{d}}$. For $\epsilon>0$, the set
$\mathcal{C}\subset\mathcal{D}$ is an $\epsilon$-cover for $\mathcal{D}$ if
for every $x\in\mathcal{D}$ there exists a $c_{i}\in\mathcal{C}$ such that
$\left|\left|x-c_{i}\right|\right|\leq\epsilon$. Equivalently,
$\mathcal{D}\subset\bigcup_{i}\mathcal{B}_{\epsilon}(c_{i})$.

In the context of a Gaussian radial basis neural network or Gaussian kernel
machines, $\epsilon$-covers partition the original space in a predictable
manner in the kernel mapped Hilbert space. An equivalent $\epsilon$-cover
specification is the following $\infty$-norm inequality applied to the kernel
mapping output:

${\rm{C}}=\left\\{\,x\in\mathbb{R}^{n_{d}}\;|\
\left|\left|{\vec{k}}(x)\right|\right|_{\infty}\geq\varphi(\epsilon;\sigma)=e^{-\epsilon^{2}/\sigma^{2}}\,\right\\}.$
(5)

If the KEI set is defined based on an $\epsilon$-cover of
$\mathcal{D}\subset\mathbb{R}^{n_{d}}$, then the kernel map will map points in
$\mathcal{D}$ to points in $\mathcal{H}^{n_{c}}$ on a surface some minimal
distance from the origin. Points in $\mathbb{R}^{n_{d}}$ receding from
$\mathcal{D}$ will tend towards the origin in $\mathcal{H}^{n_{c}}$.

The top row in Figure 1 depicts a 2D example consisting of a collection of
safe (blue +) and unsafe (orange $\circ$) points generated from a nonlinear
separating boundary (dashed black curve), as well as three GRBFs (black *)
centered in the unsafe region, $\mathcal{X}^{u}\subset\mathcal{D}$, that are
an $\epsilon$-cover. These three centers create a mapping of the 2D domain
into a 3D Hilbert space, depicted in the bottom row where the orange surface
is the unsafe region surface ${\vec{k}}(\mathcal{U}^{*})$ and the blue surface
is the safe region surface ${\vec{k}}(\mathcal{S}^{*})$. The cutting surface
sought is one that separates points “near” to the origin from those in
${\vec{k}}(\mathcal{U}^{*})$. Two such surfaces are depicted in the second and
third columns (bottom row) along with their level-set boundaries in the input
space (top row). The next section describes how to derive cutting surfaces
from solutions to constrained optimization problems.

##  III Cutting Surfaces and Partitions

This section exploits the geometry of the kernel mapping to define cutting
hyperplanes in its output Hilbert space. These hyperplanes translate to
partitions of the input space. Initially, the cutting hyperplanes will be in
the original Hilbert space and are defined by a single-layer kernel machine
network. Next, a second layer is added to the kernel machine network to create
cutting hyperplanes that translate to elliptical volumes in the original
Hilbert space. More complex cutting surfaces are then explored based on the
inclusion of positive (unsafe) and negative (safe) samples.

###  III-A Single-Layer Kernel Machine Network Partitions

Constructing an implicit partition function with barrier function properties
with a single-layer kernel machine involve identifying a cutting hyperplane
for the collection of safe and unsafe data. Consider a kernel embedding built
using the unsafe data where $\mathcal{C}$ is a KEI set and covering of the
unsafe set for some $\epsilon$. One cutting hyperplane results from the
following linear program:

$\displaystyle\min_{\alpha\in\mathbb{R}^{n_{c}}}$
$\displaystyle\vec{1}\,^{T}\alpha$ (6) s.t.
$\displaystyle\alpha^{T}\vec{k}(x_{i})\geq 1,\;\forall i\in\left\\{\,l\;|\
x_{l}\in\mathcal{X}^{u}\,\right\\}$ $\displaystyle\alpha_{j}\geq 0,\;\forall
j=1,\ldots,n_{c}$

where $\vec{1}=[1\,\cdots\,1]^{T}\in\mathbb{R}^{n_{c}}$, and
$\alpha\in\mathbb{R}^{n_{c}}$ are the coefficients of the separating
hyperplane. The one level-set of $f_{\text{KM}}(x)=\alpha^{T}\vec{k}(x)$
defines the hyperplane in the positive hyperoctant that is furthest from the
origin and places all unsafe points in $\mathcal{X}$ in the positive half-
plane. The linear program is guaranteed to have a solution:

######  Theorem 2.

Given a kernel embedding ${\vec{k}}$ defined from an $\epsilon$-covering of
$\mathcal{X}^{\text{u}}$, there is a hyperplanar splitting of
$\mathcal{H}^{n_{c}}$ described by

$\mathcal{S}=\left\\{\,x\in\mathcal{D}\;|\
{\left<\,\vec{b}\;,\;{\vec{{k}}}(x)\,\right>}<1\,\right\\}$ (7)

where $\mathcal{S}\subset\mathcal{S}^{*}\subset\mathcal{D}$ and
$\mathcal{X}^{\text{u}}\subset\bar{\mathcal{S}}$.

###### Proof.

Consider the covering set of centers $\mathcal{C}$. Generate a clustering of
the unsafe data points $\mathcal{X}^{\text{u}}$ based on the cluster
assignment function
$x^{u}\in\mathcal{X}^{\text{u}}\mapsto\operatorname*{arg\,max}_{i}k(x^{u},c_{i})$
for $c_{i}\in\mathcal{C}$. For each cluster set $\mathcal{X}^{\text{u}}_{i}$
find the minimum value $y^{i}=k(x,c_{i})$ for
$x\in\mathcal{X}^{\text{u}}_{i}$. Each cluster minimum defines the $i$-th
coordinate intercept for a cutting plane in the Hilbert space, which defines
$\vec{b}$. All unsafe points $\mathcal{X}^{\text{u}}$ lie on the non-negative
side, possibly with some safe points from $\mathcal{X}^{\text{s}}$. Only safe
points in $\mathcal{X}^{\text{s}}$ lie on the negative side (it is half-plane
containing origin). ∎

######  Corollary 4.

By virtue of defining a splitting with a kernel embedding, this same operation
generates a splitting of the original domain $\mathcal{D}$, and likewise the
full space $\mathbb{R}^{n_{d}}$.

The second column of Figure 1 depicts the separating hyperplane, generated by
the linear program in (6). The bottom plot shows the hyperplane (yellow
surface) in the Hilbert space, and the top plot shows the separating boundary
(yellow closed contour) in the input domain.

###  III-B Two-Layer Neural Network Partitions

The Hilbert space splitting based on a separating hyperplane may not capture
the true classification boundary in the Hilbert space. Consequently, a richer
decision boundary based on a nonlinear separating surface should create a
splitting with fewer misclassified points. Consider the problem of defining a
quadratic safety volume in the Hilber space specified by a quadratic boundary

$\zeta^{T}A_{2}\zeta+b_{2}^{T}\zeta+c_{2}=0,$ (8)

where elements below the surface (evaluate to less than 0) are guaranteed to
be safe, and unsafe elements are guaranteed to be above the surface (evaluate
to greater than 0). Now, the optimization problem becomes a quadratically
constrained linear program (QCLP), which is an NP-hard problem in theory [22,
23]. Although relaxation techniques for solving QCLPs exist, such optimization
problems should be avoided when possible. Adding a quadratic polynomial kernel
as a second layer to the network will preserve the LP formulation in (6) for
the cutting surface optimization problem.

####  III-B1 A Quadratic Polynomial Layer

Using a quadratic polynomial kernel of the form
$p_{2}(x,y)=(x^{T}y+\lambda)^{2}$ leads to quadratic cutting surfaces in
$\mathcal{H}^{n_{q}}$ where $n_{q}$ is the number of kernels. The mapping is
now ${\vec{p}}_{2}\circ\vec{k}:\mathcal{D}\rightarrow\mathcal{H}^{n_{q}}$,
where

$\vec{p}_{2}(x)=\left[p_{2}(x,y_{1})\;\cdots\;p_{2}\left(x,y_{n_{q}}\right)\right]^{T}.$
(9)

For the quadratic polynomial layer, the parametric degrees of freedom are the
vectors $y_{i}$ such that each kernel coordinate is $p_{2}(x,y_{i};\lambda)$.
For simplicity, let $n_{q}\geq n_{c}$ such that $y_{i}=e_{i}$ for
$i\in\left\\{\,1,\dots,n_{c}\,\right\\}$ where $e_{i}$ is the $i^{\text{th}}$
unit coordinate vector. The remaining vectors $y_{i}$, if any, for
$n_{c}<i\leq n_{q}$ can be any basis elements that support the task at hand
(typically selected based on the training data). Creating a separating
boundary is equivalent to establishing the coefficients $\alpha_{i}$ for the
following constraint equation,

$z^{T}\left(\sum_{i}\alpha_{i}y_{i}y_{i}^{T}\right)z+2\lambda\left(\sum_{i}\alpha_{i}y_{i}^{T}\right)z+\lambda^{2}\left(\sum\alpha_{i}\right)=0,$
(10)

where $z\in\mathcal{H}^{n_{c}}$. The structure of (10) matches that of (8).
Note that due to the limited number of kernels, and thus limited number of
basis vectors $y_{i}$, (10) can only represent a subset of solutions generated
by (8). However, the trade-off is that the optimization problem for solving
the coefficients of kernel machines remains a linear program:

$\displaystyle\min_{\alpha\in\mathbb{R}^{n_{q}}}$
$\displaystyle\vec{1}\,^{T}\alpha$ (11) s.t.
$\displaystyle\alpha^{T}\vec{p}_{2}\circ\vec{k}(x_{i})\geq 1,\;\forall
i\in\left\\{\,l\;|\ x_{l}\in\mathcal{X}^{u}\,\right\\}$
$\displaystyle\alpha_{j}\geq 0,\;\forall j=1,\ldots,n_{q}$

######  Theorem 3.

Given a kernel embedding ${\vec{k}}$ defined from an $\epsilon$-covering of
$\mathcal{X}^{\text{u}}$, there is a hyperspherical splitting of
$\mathcal{H}^{n_{c}}$ described by

$\mathcal{S}=\left\\{\,x\in\mathcal{D}\;|\
\left|\left|{\vec{{k}}}(x)+\lambda\vec{1}\right|\right|_{2}<\rho_{P}\
\text{for}\ x\in\mathcal{X}^{\text{u}}\,\right\\}$ (12)

where $\mathcal{S}\subset\mathcal{S}^{*}\subset\mathcal{D}$, $\lambda\geq 0$,
$\rho_{P}=\min\left(\left|\left|{\vec{{k}}}(x)+\lambda\vec{1}\right|\right|_{2}\right)$
for $x\in\mathcal{X}^{\text{u}}$, and
$\mathcal{X}^{\text{u}}\in\bar{\mathcal{S}}$.

###### Proof.

The theorem asserts the existence of a feasible quadratic surface with hard
unsafe constraints for splitting the Hilbert space using the 2-norm. The
2-norm operation is equivalent to setting $\lambda=0$,
$\alpha_{i}=\rho_{P}^{-2}$ for $i\leq n_{c}$, and $\alpha_{i}=0$ for
$n_{c}<i\leq n_{q}$ (should such $i$ exist) for the polynomial quadratic
kernel problem specified by (11). Thus, the set of feasible solutions to (11)
for $\lambda=0$ has at least one element in it. For $\lambda>0$,
$\rho_{P}=\min\left(\left|\left|{\vec{{k}}}(x)+\lambda\vec{1}\right|\right|_{2}\right)$
and gives a hyper-sphere centered at $-\lambda\vec{1}$. ∎

The theorem establishes the existence of a solution to the linear program
defined in (11), thereby showing that the solution space is non-empty when
$\lambda\geq 0$. The case of $\lambda<0$ is possible but trickier to solve
for, and may lead to poor solutions when $-1<\lambda<0$. Referring again to
the example in Figure 1, the third column depicts the separating quadratic
surface, solved by the linear program in (11). The bottom plot shows the
quadratic surface (yellow surface) in the Hilbert space, and the top plot
shows the separating boundary (yellow closed contour) in the input domain.
There are slightly fewer misclassified safe data points compared to the linear
hyperplane results in the second column.

Higher order polynomial kernels may be used to partition the space. Doing so
should provide more space carving degrees of freedom but will require a policy
for selecting the kernel basis elements (beyond the unit coordinate vectors of
the Hilbert space). Solving for the second layer as a kernel SVM using a
sequential minimal optimization (SMO) algorithm [24] will provide a solution
that identifies the basis elements and solves for their coefficients, much
like in [17]. This section focused on a constructive approach to guarantee a
non-empty solution space; it can be augmented with an SMO-like basis expansion
solver to obtain $n_{q}>n_{c}$.

###  III-C Two-Layer Partitions with Safe and Unsafe Samples

The linear programs for the splitting hyperplane and hyperellipsoid in the
lifted Hilbert space $\mathcal{H}^{n_{c}}$ attempt to find the furthest
hyperplane or largest hyperellipsoid. They may be conservative relative to
other options for synthesizing such a boundary due to the curvature of the
Hilbert space or the shape of the separating boundary. As a result,
misclassification of safe points as unsafe will occur, as seen in Figure 1.
Richer surfaces exist that better capture the regions misclassified by the
linear and ellipsoidal cutting surfaces. However, the linear program
establishing the model coefficients $\alpha$ will require knowledge of what
parts of the mapped space should be considered safe versus unsafe, and require
having sufficient samples of both classes. These samples permit recovery of
more complex cutting surfaces.

The coordinate vectors, $y_{i}$, for the polynomial kernels in the multi-order
polynomial layer should at least include $n_{c}$ basis unit vectors for each
polynomial order (i.e., the minimum LP problem size grows linearly with the
maximum polynomial order). These can be complemented with a specialized or
task-specific basis expansion criteria. The polynomial kernel mapping is the
following,

$\displaystyle\vec{p}(z)=[$ $\displaystyle(y_{1}^{T}z+\lambda_{1})\ \cdots\
(y_{n_{c}}^{T}z+\lambda_{1})\ \cdots\ $ (13)
$\displaystyle(y_{1}^{T}z+\lambda_{p_{i}})^{p_{i}}\ \cdots\
(y_{n_{c}}^{T}z+\lambda_{p_{i}})^{p_{i}}\ \cdots]^{T},$

where $p_{i}\in\left\\{\,1,\dots,n_{p}\,\right\\}$, $n_{p}$ is the maximal
polynomial order used, and $N_{p}$ is the total number of polynomial kernels
used across all orders. The minimum basis elements requirement means that
$N_{p}\geq n_{p}n_{c}$. The linear-in-$\alpha$ constraint equation for a
multi-order polynomial surface cut is

$\alpha^{T}p(z)=\sum_{i}\alpha_{i}\left(y_{i}^{T}z+\lambda_{p_{i}}\right)^{p_{i}}=0.$
(14)

Per [17], the inclusion of positive and negative samples should involve hard
constraints on the unsafe points and soft constraints for the safe points

$\displaystyle\min_{\alpha\in\mathbb{R}^{N_{p}},\xi\in\mathbb{R}^{n_{s}}}\quad$
$\displaystyle\vec{1}\,^{T}\xi$ (15) s.t.
$\displaystyle\alpha^{T}\vec{p}\circ\vec{k}(x_{i})\geq 1,\ \forall
i\in\left\\{\,l\;|\ x_{l}\in\mathcal{X}^{u}\,\right\\}$
$\displaystyle\alpha^{T}\vec{p}\circ\vec{k}(x_{j})\leq-1+\xi_{j},$
$\displaystyle\xi_{j}\geq 0\,\quad\qquad\qquad\forall j\in\left\\{\,l\;|\
x_{l}\in\mathcal{X}^{s}\,\right\\}$

######  Theorem 4.

Defining the two layer neural network $f=\vec{p}\circ\vec{k}$ with $n_{p}\geq
2$ and performing a hyper-planar partitioning of the space in the output space
provides an equivalent or better splitting of the space, based on safe
misclassification counts.

###### Proof.

The solution space for (15) contains the linear and quadratic cases by setting
the appropriate coefficients $\beta_{i}$ to zero. The optimization problem
will either return one of these options or it will find a better one. By
virtue of having a slack minimizing cost, better solutions will involve either
the same slack or smaller. If more slack variables evaluate to zero, there
will be less safe point classification errors. ∎

The original two-layer formulation in [17] uses a 2-layer GRBF network created
using a kernel-SVM process for the second layer. The Gaussian kernel in the
second layer defines an alternative nonlinear space for generating the cutting
surface, which may provide a richer cutting surface solution space. The
theorems and corollaries here provide a more formal analysis of the
partitioning properties of the 2-layer kernel machine, which is an instance of
a shallow layer kernel machine (SKM), and the specification of learning
problems on the output layer that are linear program formulations with a non-
empty solution space. Similar results apply for a Gaussian kernel second
layer.

(a)(b)(c)$n_{p}=2$$n_{p}=3$ Figure 2: Depiction of cutting surfaces and
separating curves for two-layer, multi-order polynomial partitions using safe
and unsafe data points.

##### Discussion

The provision of negative (safe) samples changes the structure of the problem
by adding a point set complementary to the unsafe set, deemed to be safe.
Incorporating this data into the optimization problem adds information about
what regions of the Hilbert space $\mathcal{H}^{n_{c}}$ should lie on the
negative side of the boundary, thereby pushing the boundary outwards when
there are negative samples on the other side of the hyperplanar or
hyperellipsoidal splitting (based on available degrees of freedom in the
polynomial kernel layer). Solutions should improve the accuracy of the
estimated bound and reduce the quantity of safe misclassification errors, per
Theorem 4.

Figure 2 depicts the cutting surfaces for the case of $n_{p}=2$ (left plot)
and $n_{p}=3$ (middle plot), and their input domain boundaries (right plot)
for the same example problem data of Figure 1. Both surfaces were generated
using positive and negative samples per (15). There are less misclassified
safe data points, with the $n_{p}=3$ case achieving zero slack cost.

##### Maximally Flexible Shallow Network Design

The constructions to date are based on a covering of the positive set (unsafe
set). For potentially more accurate results, additional basis centers can be
chosen from the negative (safe) set. Extremizing coordinates associated to
these centers indicate movement towards safe regions and away from unsafe.
Their inclusion changes the structure of the optimization since these are
points in the Hilbert space surface to avoid. The optimization problem
specification should seek to generate a cutting surface that carves out such
regions. Because of the flipped nature, a bias term will be needed. The bias
term can be easily appended to the multi-order polynomial layer mapping,
updating $\vec{p}$ from (13) to

$\begin{split}\vec{p}(z)=[&(y_{1}^{T}z+\lambda_{1})\ \cdots\
(y_{n_{c}}^{T}z+\lambda_{1})\ \cdots\ \\\
&(y_{1}^{T}z+\lambda_{p_{i}})^{p_{i}}\ \cdots\
(y_{n_{c}}^{T}z+\lambda_{p_{i}})^{p_{i}}\ \cdots\ 1]^{T},\end{split}$ (16)

where, $p_{i}\in\\{1,\dots,n_{p}\\}$. The optimization formulation (15) still
applies. Repeating the earlier analysis will show that (15) has at least one
point in the solution space for separating hyperplanes (and one for separating
hyperellipsoids). Consequently, (15) has a non-empty solution space.

LinearQuadraticMulti-Poly-2Multi-Poly-3Multi-Scale-Poly-2(a)(b)(c)(d)(e)
Figure 3:  Depiction of SKM-ZBF values as a color map along with its zero
level-set (yellow contour) for different network constructions. (a) and (b)
are replica of Figure 1(b,c) (top row); (c) and (d) are replica of Figure
2(c); (e) is the same as (c) but with multi-scale Gaussian kernels.

##  IV Construction of the Barrier Function

The previous sections described methods to synthesize level-set functions for
approximating the boundary of the safe set from labeled, sampled data. The
level-set functions have the necessary properties to serve as zeroing barrier
functions after a constant shift and rescaling. We call the barrier function
designed from concatenated kernel machines a shallow kernel machine-zeroing
barrier function (SKM-ZBF). An additional extended class $\mathcal{K}$ outer
function $\psi\in\mathcal{K}_{e}$ may be needed to adjust the sensitivity
based on the application, e.g., $h(x)=\psi_{0}\circ\bar{h}(x)$, where
$\bar{h}(\cdot)$ is the SKM-ZBF. This section connects the cutting surface
layer with ZBF properties and concludes with some practical considerations for
implementing the SKM-ZBF.

###  IV-A Suitability as Zeroing Barrier Functions

To serve for safety-critical control applications, there are conditions on the
level-set function to be a valid ZBF.

##### Monotonicity

ZBFs must be monotonic when evaluated across the boundary. The specification
of the boundary as a cutting surface constraint means that the domain surface
and the constraint surface are transverse to each other. The intersection of
the two surfaces occur at the zero level-set. Movement on the (Hilbert space)
domain surface away from the boundary necessarily has monotonic behavior,
locally. The hard constraint in (15) for the unsafe samples ensures local
monotonicity, at least from the -1 to +1 level-sets (0 to +2 for the
equivalent SKM-ZBF) and for some non-trivial band away from their boundaries.
Figure 3 depicts SKM-ZBF as a color map for various architectures, all of
which exhibit monotonicity across the boundary.

##### Continuously Differentiable

The canonical control barrier function constraint used for safety depends on
the gradient of the barrier function [2]. The function needs to be
continuously differentiable to avoid undesirable behavior in the constraint.
The Gaussian and polynomial layers consist of smooth basis functions, thus
their composition and linear combination is also smooth, thereby resulting in
smooth gradients of the ZBF. The gradient has a closed form solution based on
equations (2) and (16) which can be used in the CBF constraint [2] to enforce
safety.

##### Dead Gradients

Based on the CBF constraint discussed in [2], vanishing gradients for $h(x)$
for $x\in\mathcal{D}$ may lead to loss of control [25], which is undesirable
in safety critical applications. Monotonicity in the vicinity of the zero
level-set guarantees no vanishing gradients in that region. Due to the nature
of the first layer and the mapping of $\mathbb{R}^{n_{d}}$ to a compact
surface in $\mathcal{H}^{n_{c}}$, the function has known limits, which give
local extrema. One such will be the output of $z=0$ (points at infinity in the
original input space). These points will tend to the same ZBF value (visible
as nearly constant regions in the level-set plots of Fig. 3). At the other
extreme will be the center locations; they locally maximize distance from
$z=0$ in $\mathcal{H}^{n_{c}}$, and may have similar extremizing properties
relative to the transverse cutting surface. These points may be local extrema
of the ZBF, in which case the gradients will vanish there (visible as local
minima in the level-set plots of Fig. 3). Since the focus of this SKM-ZBF
construction is on spatially meaningful or location-based ZBFs related to
navigation, the existence of the aforementioned extrema on sets of measure
zero do not have a strong impact as for other use cases of ZBFs. Due to the
possible topology of unsafe space, attempts to define one signed distance
function for all disconnected regions will result in an extrema set of zero
measure (related to junctions in the Voronoi partition).

##### Guarantees on the Safe Set

The ability to capture the interface between safe and unsafe regions is a
function of the measurement density. There is an operating assumption of an
$\epsilon$-covering for the unsafe regions where the $\epsilon$ value is
related to the bandwidth of the Gaussian kernels, $\sigma$. If the sensor
resolution is too coarse to capture boundary variation, then the unsafe set
may not lie within the resulting covering and the asserted guarantees cannot
hold. It is impossible to guarantee safety when the sensors cannot measure or
sample the space as needed (unless it is known how much the local structure
can change for small changes in sample location).

The optimization problem guarantees that
$h(\mathcal{X}^{\text{u}})\subset\mathbb{R}^{-}$, which implies the existence
of a point-dependent closed covering ${\bar{\rm{C}}}(\mathcal{X}^{\text{u}})$
of the unsafe points for which
$h({\bar{\rm{C}}}(\mathcal{X}^{\text{u}}))\subset\mathbb{R}^{-}$. Assuming
that the sensor can indeed measure and capture the local structure, and the
GRBF layer reflects it, then
$\mathcal{U}^{*}\subset{\bar{\rm{C}}}(\mathcal{X}^{\text{u}})\subset\mathcal{D}$
and $\mathcal{S}\subset\mathcal{S}^{*}\subset\mathcal{D}$.

###  IV-B Practical Method for Barrier Function Synthesis

The standard kernel machine mapping consists of a single spatial scale, which
places a limit on the boundary that can be captured for a given
$\epsilon$-cover. This limit is related to the Fourier spectral properties of
the functions generated by the kernel machine [26]. Multi-scale
implementations improve this limitation [27], with each Gaussian center of the
first layer having a bandwidth parameter chosen from a finite set.

In effect, each $c_{i}$ has a bandwidth
$\sigma_{i}\in\left\\{\,\beta_{1},\,\dots\,,\beta_{N_{b}}\,\right\\}$ where
$N_{b}$ is the number of unique bandwidths used. For multiple Gaussian
bandwidths the $\epsilon$-cover concept changes to be a Gaussian kernel cover
extending the single-bandwidth version. The same equation (5) is used but the
equivalent balls $\mathcal{B}_{\epsilon_{i}}(x)$ in the original input space
will have differing radii $\epsilon_{i}$ based on the kernel function
bandwidth $\beta_{i}$ for the $i^{\text{th}}$ kernel coordinate mapping,
similar to how $p_{i}$ varies for a multi-order polynomial layer. These
(radius variable) balls should cover the space of interest. The cover
${\rm{C}}$ from (5) satisfies

$\mathcal{D}\subset{\rm{C}}=\bigcup_{i}\mathcal{B}_{\epsilon_{i}}(c_{i})\quad\text{for}\
\epsilon_{i}\in\left\\{\,\varepsilon_{1},\,\dots\,,\varepsilon_{N_{b}}\,\right\\},$
(17)

and each $\varepsilon_{j}$ depends on $\beta_{j}$. The same linear programs
for hyperplanar, hyperellipsoidal, and multi-order polynomial cutting surfaces
apply, as do the existence of at least one solution in the feasible space for
each LP specification. Earlier work [28] established a greedy selection policy
for recovering a domain covering center set $\mathcal{C}$ for optimized
function approximation. A similar policy applies to the case of multi-scale,
first layer synthesis from data, modified with finer bandwidth values based on
residual error minimization of the coarser single layer basis functions.
Likewise, regression or function recovery from Gaussian kernel machines
operate best when capturing the deviation from some known parametric model
[29], leading to a semi-parametric first layer (fixed + data-adaptive
elements).

##  V Examples and Implementation

(a)(b)(c)(d)

Figure 4: Curved lane modeling using SKM-ZBF. (a) curved road; (b) road with
labels for synthesized safe (blue) and unsafe (red) points along with GRBF
center locations (black *); (c) synthesized zero level-set of SKM-ZBF (dashed
yellow lines); (d) synthesized SKM-ZBF as a 3D surface plus zero level-set.

(c)(d)(e)(b)(a)

Figure 5: Mobile robot planar navigation problem. Synthesis of a global SKM-
ZBF using (a) multi-polynomial second layer and (b) Gaussian kernel second
layer, to define unsafe region boundaries (yellow contours). Local navigation
SKM-ZBF generated at three marked locations (gray dots) in the global maps,
are depicted in (c), (d), and (e). Orange dots are LiDAR points, the purple
circle is the robot’s location, and the yellow contour is the SKM-ZBF
boundary.

###  V-A Curved Lane Modeling

This case study shows how the proposed method generates a suitable SKM-ZBF
that models curved lanes, as shown in the low-speed driving scenario of Figure
5(a). Closed-form expressions for curved roads can be challenging to
determine, and thus, using control barrier function policy synthesis
frameworks for such applications can be difficult. To that end, we show SKM-
ZBF synthesis alleviates this issue. The safe and unsafe data points used for
network training are shown in Figure 5(b) and are generated synthetically. In
practice, such points would be generated from visual sensing and processing,
such as from semantic segmentation of data from a forward-facing camera on the
car.

A simple strategy for choosing the GRBFs centers is to uniformly place them
along the center line of the lane. Seven GRBFs with a bandwidth of
$\sigma=2w_{r}/\sqrt{\log(2)}$ are used in this case study, as shown in Figure
5(b) by the black (*); where $w_{r}$ is the width of the road. The synthesized
SKM-ZBF boundary is shown in Figure 5(c) as the yellow dashed line and as a
surface in Figure 5(d). The SKM-ZBF compute time was 10 msec (MATLAB, Ubuntu
20.04, Intel i7-8750H CPU), which supports a 100 Hz update rate.

###  V-B Planar Mobile Robot Navigation in ROS STDR Simulator

This study involves a differential drive robot equipped with a LiDAR sensor
navigating a 2D environment with obstacles, see Figure 5. The robot tracked a
predetermined, dynamically feasible, collision-free trajectory. Ideal
localization was assumed to be available for the robot. The procedure to
generate the dataset of safe and unsafe LiDAR samples was done similar to
[17]. Global and local SKM-ZBF synthesis was performed. Global synthesis
aggregated all of the sensor data and applied either a multi-polynomial second
layer, Figure 5(a), or a Gaussian second layer via kernel SVM optimization
[17], Figure 5(b). The latter generates a richer solution space and will be
more accurate. The first layer used a fixed $5\times 8$ uniform grid. Local
SKM-ZBF synthesis used a fixed robot-centered $4\times 4$ grid for three
locations, shown as gray dots in the global maps. The zero level-set of the
locally and globally synthesized SKM-ZBFs are shown as yellow curves. Notice
that the local SKM-ZBF zero level-set cuts into some sides of obstacles that
do not have LiDAR measurements. However, once measurements are taken on those
sides the SKM-ZBF will include them.

To further improve the accuracy of the SKM-ZBF in capturing the domain, multi-
scale kernels can be added to the first layer in addition to the fixed grid
kernels. The centers and bandwidths of the multi-scale kernels are turned
according to the adaptive strategy described in section IV.B.

The global SKM-ZBF compute time was 3.9 sec for the multi-polynomial SKM-ZBF
and 4.3 sec for the Gaussian SKM-ZBF when solved using the SMO solver in the
LIBSVM package (dual QP formation takes 18 sec). The local SKM-ZBF compute
time averaged 18 msec, which supports a 50 Hz update rate. The training time
of the local SKM-ZBF using the proposed LP formulation is significantly
shorter than the training time of neural network-based ZBFs. For example, the
training time of a Softplus-based deep neural network via gradient descent for
local ZBF synthesis was 98 msec on GPU and 407 msec on CPU [14].

##  VI Conclusion

This paper presented a zeroing barrier function synthesis method based on a
two-layer shallow kernel machine network architecture. The first layer is
constructed from Gaussian kernels whose geometry in the associated Hilbert
space was analyzed while considering positive samples (unsafe points) only.
The analysis provided theoretical guarantees on the existence of a separating
hyperplane, synthesized via a linear program (LP), that splits the Hilbert
space into safe and unsafe regions via its level-sets. A second layer was
added to the network that uses polynomial kernels and negative samples (safe
points) to reduce misclassifications, while also maintaining the LP
formulation. The cutting surface geometry of this second layer also applies to
the previous work [17]. Two case studies demonstrate the efficacy of the SKM-
ZBF architecture in generating valid zeroing barrier functions for motion
planning applications. Future work will explore optimized construction of SKM-
ZBFs for boundary estimation accuracy in navigation contexts.

## References

  * [1] F. Blanchini, “Set invariance in control,” Automatica, vol. 35, no. 11, pp. 1747–1767, 1999.
  * [2] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in European Control Conference, pp. 3420–3431, 2019.
  * [3] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” IEEE Trans. on Automatic Control, vol. 62, no. 8, pp. 3861–3876, 2017.
  * [4] L. Wang, A. D. Ames, and M. Egerstedt, “Safety barrier certificates for collisions-free multirobot systems,” IEEE Trans. on Robotics, vol. 33, no. 3, pp. 661–674, 2017.
  * [5] L. Wang, A. D. Ames, and M. Egerstedt, “Safe certificate-based maneuvers for teams of quadrotors using differential flatness,” in IEEE International Conference on Robotics and Automation, pp. 3293–3298, 2017.
  * [6] P. Pierpaoli, A. Li, M. Srinivasan, X. Cai, S. Coogan, and M. Egerstedt, “A sequential composition framework for coordinating multirobot behaviors,” IEEE Trans. on Robotics, vol. 37, no. 3, pp. 864–876, 2021.
  * [7] M. Srinivasan and S. Coogan, “Control of mobile robots using barrier functions under temporal logic specifications,” IEEE Trans. on Robotics, vol. 37, no. 2, pp. 363–374, 2021.
  * [8] L. Lindemann and D. V. Dimarogonas, “Decentralized control barrier functions for coupled multi-agent systems under signal temporal logic tasks,” in European Control Conference, pp. 89–94, June 2019.
  * [9] Y. Chen, M. Jankovic, M. Santillo, and A. D. Ames, “Backup control barrier functions: Formulation and comparative study,” in IEEE Conference on Decision and Control, pp. 6835–6841, 2021.
  * [10] J. Schmidt-Hieber, “The kolmogorov–arnold representation theorem revisited,” Neural Networks, vol. 137, pp. 119–126, 2021.
  * [11] A. Robey, H. Hu, L. Lindemann, H. Zhang, D. V. Dimarogonas, S. Tu, and N. Matni, “Learning control barrier functions from expert demonstrations,” in Conference on Decision and Control, pp. 3717–3724, 2020.
  * [12] Y. Meng, Z. Qin, and C. Fan, “Reactive and safe road user simulations using neural barrier certificates,” in IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 6299–6306, 2021.
  * [13] V. Dhiman, M. J. Khojasteh, M. Franceschetti, and N. Atanasov, “Control barriers in bayesian learning of system dynamics,” IEEE Trans. on Automatic Control, pp. 1–1, 2021.
  * [14] K. Long, C. Qian, J. Cortés, and N. Atanasov, “Learning barrier functions with memory for robust safe navigation,” IEEE Robotics and Automation Letters, vol. 6, no. 3, pp. 4931–4938, 2021.
  * [15] K. Long, V. Dhiman, M. Leok, J. Cortés, and N. Atanasov, “Safe control synthesis with uncertain dynamics and constraints,” IEEE Robotics and Automation Letters, vol. 7, no. 3, pp. 7295–7302, 2022.
  * [16] F. Ramos and L. Ott, “Hilbert maps: Scalable continuous occupancy mapping with stochastic gradient descent,” The International Journal of Robotics Research, vol. 35, no. 14, pp. 1717–1730, 2016.
  * [17] M. Srinivasan, A. Dabholkar, S. Coogan, and P. A. Vela, “Synthesis of control barrier functions using a supervised machine learning approach,” in IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 7139–7145, 2020.
  * [18] J. Park and I. W. Sandberg, “Universal approximation using radial-basis-function networks,” Neural Computation, vol. 3, no. 2, pp. 246–257, 1991.
  * [19] D. B. McDonald, W. J. Grantham, W. L. Tabor, and M. J. Murphy, “Global and local optimization using radial basis function response surface models,” Applied Mathematical Modelling, vol. 31, no. 10, pp. 2095–2110, 2007.
  * [20] K. M. Lynch and F. C. Park, Modern Robotics: Mechanics, Planning, and Control.  USA: Cambridge University Press, 1st ed., 2017.
  * [21] K. Muandet, K. Fukumizu, B. Sriperumbudur, and B. Schölkopf, “Kernel mean embedding of distributions: A review and beyond,” Foundations and Trends in Machine Learning, vol. 10, no. 1-2, pp. 1–141, 2017.
  * [22] J. Linderoth, “A simplicial branch-and-bound algorithm for solving quadratically constrained quadratic programs,” Mathematical Programming, vol. 103, pp. 251–282, Jun 2005.
  * [23] K. Chatterjee, H. Fu, A. K. Goharshady, and E. K. Goharshady, “Polynomial invariant generation for non-deterministic recursive programs,” in Proceedings of the ACM SIGPLAN Conference on Programming Language Design and Implementation, p. 672–687, 2020.
  * [24] J. Platt, “Sequential minimal optimization: A fast algorithm for training support vector machines,” Tech. Rep. MSR-TR-98-14, Microsoft, April 1998.
  * [25] M. Srinivasan, N.-s. P. Hyun, and S. Coogan, “Weighted polar finite time control barrier functions with applications to multi-robot systems,” in IEEE Conference on Decision and Control, pp. 7031–7036, 2019.
  * [26] F. Girosi, M. Jones, and T. Poggio, “Regularization theory and neural networks architectures,” Neural Computation, vol. 7, no. 2, pp. 219–269, 1995.
  * [27] H. Wendland, Scattered Data Approximation.  Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 2004.
  * [28] H. Kingravi, P. Vela, and A. Gray, “Reduced set KPCA for improving the training and execution speed of kernel machines,” in SIAM International Conference on Data Mining, pp. 441–449, 2013.
  * [29] A. Chang, C. Hubicki, J. Aguilar, D. Goldman, A. Ames, and P. Vela, “Learning terrain dynamics: A Gaussian process modeling and optimal control adaptation framework applied to robot jumping,” IEEE Trans. on Control Systems Technology, vol. 29, no. 4, pp. 1581–1596, 2020.


# Leveraging arbitrary mobile sensor trajectories with shallow recurrent
decoder networks for full-state reconstruction

Megan R. Ebers†, Jan P. Williams†, Katherine M. Steele†, J. Nathan Kutz∗,∗∗
†Department of Mechanical Engineering, University of Washington, Seattle, WA
98195  ∗Department of Applied Mathematics, University of Washington, Seattle,
WA 98195  ∗∗Department of Electrical and Computer Engineering, University of
Washington, Seattle, WA 98195

###### Abstract

Sensing is one of the most fundamental tasks for the monitoring, forecasting
and control of complex, spatio-temporal systems. In many applications, a
limited number of sensors are mobile and move with the dynamics, with examples
including wearable technology, ocean monitoring buoys, and weather balloons.
In these dynamic systems (without regions of statistical-independence), the
measurement time history encodes a significant amount of information that can
be extracted for critical tasks. Most model-free sensing paradigms aim to map
current sparse sensor measurements to the high-dimensional state space,
ignoring the time-history all together. Using modern deep learning
architectures, we show that a sequence-to-vector model, such as an LSTM (long,
short-term memory) network, with a decoder network, dynamic trajectory
information can be mapped to full state-space estimates. Indeed, we
demonstrate that by leveraging mobile sensor trajectories with shallow
recurrent decoder networks, we can train the network (i) to accurately
reconstruct the full state space using arbitrary dynamical trajectories of the
sensors, (ii) the architecture reduces the variance of the mean-square error
of the reconstruction error in comparison with immobile sensors, and (iii) the
architecture also allows for rapid generalization (parameterization of
dynamics) for data outside the training set. Moreover, the path of the sensor
can be chosen arbitrarily, provided training data for the spatial trajectory
of the sensor is available. The exceptional performance of the network
architecture is demonstrated on three applications: turbulent flows, global
sea-surface temperature data, and human movement biomechanics.

##  I Introduction

Figure 1: Summary figure of a shallow recurrent decoder network (SHRED)
leveraging mobile sensors to reconstruct full state-space estimates from
sparse dynamical trajectories. (Left) Sensor trajectory history encodes global
information of the spatio-temporal dynamics of the sparsely measured system.
In this work, we evaluate three challenging datasets, including forced
isotropic turbulence, global sea-surface temperature, and human biomechanics.
(Middle) The mobile SHRED architecture can (i) embed the multiscale physics of
a system into a compact and low-dimensional latent space, and (ii) provide a
mapping from the sparse mobile sensors to a full state estimate. (Right) The
high-dimensional and complex system states can be reconstructed, provided
training data for the dynamical trajectory of the sensor(s) is available.

Sensing is a ubiquitous task, being critically important to every scientific
and engineering discipline Manohar  _et al._ (2018). Sensor technologies
accelerate scientific advancement by providing improved observational data
that is capable of advancing new scientific theories and improving
technological development. From the James Webb telescope Gardner  _et al._
(2006) to single neuron recordings Stevenson and Kording (2011), sensors
provide an interface with the physical world that allows us to interrogate
diverse and complex systems across scientific domains. In many application
fields, sensors are mobile by design, moving with dynamical trajectories as
they record quantities of interest. Application areas such as motion tracking
of human biomechanics or robotic systems, buoys that measure ocean dynamics,
and weather balloons are all examples of technologies where sensors are not
(or cannot be) immobile. From these mobile sensors, some or all of the
canonical tasks of sensing are enacted, including reconstruction, forecasting,
model discovery, control and uncertainty quantification Manohar  _et al._
(2018). Emerging sensor technologies are capable of producing exceptional
quantities of data which can now be leveraged by machine learning algorithms.
Specifically, we utilize a recurrent neural network Medsker and Jain (2001);
Salehinejad  _et al._ (2017); Zaremba  _et al._ (2014) to learn the temporal
sequences in a latent space that maps to full state-space estimates via a
shallow decoder network Erichson  _et al._ (2020); Williams  _et al._ (2022).
This shallow recurrent decoder network (SHRED) architecture has been
successful for immobile sensors Williams  _et al._ (2023), but this work shows
how it can be used in a mobile sensing framework. Indeed, mobile sensing
affords distinct advantages over immobile sensors, including the reduction of
the variance of the mean-square reconstruction error in comparison with
immobile sensors and improved generalization (parameterization of dynamics)
for data outside the training set.

Traditional system identification techniques such as the Kalman filter and its
variants have long achieved state-of-the-art performance Law  _et al._ (2015),
but are limited to applications where a model of the underlying dynamics is
readily available. With the recent explosion in quantity and quality of data
from high-dimensional systems, we are increasingly required to perform state
estimates of systems with unknown or computationally infeasible dynamics.
State estimation in these cases is typically achieved by exploiting low-rank
structures in the data learned by a singular value decomposition (SVD),
otherwise known as proper orthogonal decomposition (POD) Kutz (2013); Brunton
and Kutz (2019). Gappy POD was one of the first such techniques to reconstruct
high-dimensional fields from limited, sparse sensor measurements Everson and
Sirovich (1995). Since then, much work has been aimed at optimizing sensor
placement to well-condition the linear inverse problem upon which state
estimates rely Manohar  _et al._ (2018); Drmač and Gugercin (2016); Barrault
_et al._ (2004); Chaturantabut and Sorensen (2010). There exist numerous
techniques for determining optimal sensor locations in small search spaces
Boyd and Vandenberghe (2004); Joshi and Boyd (2009); Caselton and Zidek
(1984); Krause  _et al._ (2008); Lindley (1956); Sebastiani and Wynn (2000);
Paninski (2005), but in high-dimensional systems a greedy approximation, such
as QR decomposition with column pivoting, is needed. These techniques for
state-estimation in the absence of a governing model also differ from Kalman
filtering in that they only utilize static measurements, rather than condition
upon a history of sensor data. While this lends the methods flexibility to
reconstruct data from a dataset of faces or a dynamical system, it ignores the
rich information that can be gleaned from incorporating past sensor
measurements in a method for reconstruction.

Leveraging the temporal histories of sensor measurements offers an opportunity
to capitalize on both the rich information from measurement trajectories and
enable greater inclusion of measurement modalities – beyond immobile sensors.
To date, the inclusion and design of dynamic trajectories are addressed with
sensor path planning. Path planning is a long-standing challenge in the field
of engineering and robotics related to navigation and estimation in dynamical
environments Gunnarson  _et al._ (2021); Krishna  _et al._ (2022); Buzzicotti
_et al._ (2020); Biferale  _et al._ (2019); Madridano  _et al._ (2021). The
majority of path planning has focused on the modeling and control of the
mobile sensor positions Shriwastav  _et al._ (2022); Lynch  _et al._ (2008);
Leonard  _et al._ (2007); DeVries  _et al._ (2013); Ogren  _et al._ (2004);
Zhang and Leonard (2010); Paley and Wolek (2020); Peng  _et al._ (2014). For
example, recent work proposed the use of mobile sensors with Kalman filtering
to estimate spatio-temporal data Mei  _et al._ (2022). The authors’ approach
to sensor path planning for dynamic estimation includes two main steps: (i)
dynamic mode decomposition finds a low-rank representation of the data for
Kalman filtering, and (ii) greedy path finding optimizes the observability
matrix along the path and improves Kalman filter estimation. However, path
planning for mobile sensing and state estimation may be unnecessary. The SHRED
architecture has demonstrated the successful reconstruction of complex spatio-
temporal data with arbitrarily placed – albeit immobile – sensors Williams
_et al._ (2023). Therefore, this work aims to employ mobile sensors with SHRED
to understand whether the use of arbitrary dynamic trajectories can feasibly
and reliably reconstruct complex spatio-temporal data, thereby removing the
challenge of path planning altogether.

As a mathematical architecture, SHRED leverages emerging deep learning
paradigms LeCun  _et al._ (2015) and their universal approximation
capabilities Hornik  _et al._ (1989) to map from sparse and minimal sensor
measurements to high-dimensional spatio-temporal data. In contrast with common
sensing modalities where principled sensor placement is critical Manohar  _et
al._ (2018), the LSTM embedded in the SHRED architecture extracts a latent
representation of the low-dimensional dynamics using sensor trajectory (time
history) information which is agnostic to the sensor location. Thus the sensor
trajectory history encodes global information of the spatio-temporal dynamics
of the measured system. And while previous research has considered mobile
sensor features like the timescale of the spatio-temporal dynamics, velocity
of the sensors, and rate of sampling as part of the path planning optimization
Mei  _et al._ (2022), in this work, such features are treated hyperparameters
to be tuned in the SHRED architecture. We demonstrate the performance of SHRED
on three challenge data sets related to turbulence, sea-surface temperature,
and human biomechanics. In all three cases, the mobile SHRED architecture
provides a high-quality and accurate algorithm for estimating the original
complex, multiscale, and high-dimensional system. In all the examples, it is
demonstrated that by leveraging mobile sensor trajectories with SHRED, we can
train the network to accurately reconstruct the full state space while also
reducing the variance of the mean-square error of the reconstruction error in
comparison with immobile sensors. Moreover, the path of the sensor can be
chosen arbitrarily, provided training data for the spatial trajectory of the
sensor is available. SHRED works with high-probability given that the data do
not have regions of statistically independent data. In addition, when full
state data is not available, proxy computational data with empirically similar
data can be used for training.

##  II Mathematical Formulation

The shallow recurrent decoder network (SHRED) can be understood as the
amalgamation of an LSTM for processing a time-series of sensor measurements
followed by a feedforward neural network, or decoder, for reconstructing a
high-dimensional state from the learned latent representation of the LSTM
Williams  _et al._ (2023). Let the high-dimensional state to be reconstructed
be denoted as $x_{T}\in\mathbb{R}^{n}$ and assume access to a set of sensor
measurements $y_{t}=Cx_{t}\in\mathbb{R}^{m}$ for
$t\in\\{T-K,T-K+1,\dots,T-1,T\\}.$ $K$ can be determined by empirical analyses
of the system at hand. We assume the measurements are sparse point
measurements, that is $m<<n$ and $C$ consists of rows of the $n\times n$
identity matrix, although there is some evidence to suggest that neural
network based reconstructions can be performed with nonlinear measurements
Erichson  _et al._ (2020). The set of sensor measurements serve as inputs to
an LSTM Hochreiter and Schmidhuber (1997) with recursive update equations

$\displaystyle
h_{t}=\sigma\left(W_{o}\begin{bmatrix}h_{t-1},y_{t}\end{bmatrix}+b_{o}\right)\odot\tanh(c_{t})$
(1) $\displaystyle
c_{t}=\sigma\left(W_{f}\begin{bmatrix}h_{t-1},y_{t}\end{bmatrix}+b_{f}\right)\odot
c_{t-1}$ (2)
$\displaystyle+\sigma\left(W_{i}\begin{bmatrix}h_{t-1},y_{t}\end{bmatrix}+b_{f}\right)\odot\tanh\left(W_{g}\begin{bmatrix}h_{t-1},y_{t}\end{bmatrix}+b_{g}\right)$

where $W_{RN}=\\{W_{o},W_{f},W_{i},W_{g},b_{o},b_{f},b_{i},b_{g}\\}$ are the
trainable weights and biases of the LSTM. We denote

$h_{T}=\mathcal{G}(\\{y_{t}\\}_{T-K}^{T};W_{RN}).$ (3)

The latent state $h_{T}$ learned by the LSTM has a variety of interesting
properties that will be discussed in a later section.

The feedforward component of the SHRED architecture is a shallow decoder with
$b$ layers denoted by

$\mathcal{F}(h;W_{SD}):=R(W^{b}R(W^{b-1}\cdots R(W^{1}h))),$ (4)

parameterized by trainable weights $W_{SD}=\\{W^{1},\dots,W^{b}\\}$ and with
nonlinear scalar activation function $R$ (chosen to be ReLU). In total, the
SHRED network is given by

$\mathcal{H}(\\{y_{t}\\}_{T-K}^{T})=\mathcal{F}(\mathcal{G}(\\{y_{t}\\}_{T-K}^{T};W_{RN});W_{SD}).$
(5)

The network is trained to minimize reconstruction loss over a set of training
states $\\{x_{t}\\}_{1}^{N}$,

$\mathcal{H}\in\operatorname*{arg\rm{}min}_{\widetilde{\mathcal{H}}\in\mathscr{H}}\sum_{t=1}^{N}||x_{i}-\widetilde{\mathcal{H}}\left(\\{y_{i}\\}_{i=t-K}^{t}\right)||_{2},$
(6)

using the ADAM optimizer Kingma and Ba (2017). The assumption of access to a
set of high-dimensional states for training in this manner is a strong one;
simultaneous measurement of an entire high-dimensional system is sometimes
simply impossible. In such cases, a high-fidelity simulation can be used to
train the network, provided the simulation accurately approximates the
statistics of the real system. Alternatively, if full-state measurements are
possible, but prohibitively expensive in the long-term, the generation of
training data can be viewed as a one-time up front cost.

Previous work has demonstrated that such networks outperform traditional, POD
based techniques for state estimation while requiring fewer available sensors
Williams  _et al._ (2023). In this work, we consider the case of a time-
dependent measurement matrix $C.$ That is,

$y_{t}=C_{t}x_{t}$ (7)

where $C_{t}$ varies in time. This corresponds to a mobile sensor, in contrast
to existing work which considers only immobile sensors. We emphasize that
while the set of measurement matrices $C_{t}$ can be chosen arbitrarily,
corresponding training data is necessary in order to train the network; the
architecture as it stands cannot extrapolate to unseen sensor location
trajectories. In the cases of periodic or quasi-periodic phenomena, this
requirement amounts to the evolution of $C_{t}$ being periodic over some
characteristic time-scale in order to extrapolate beyond a temporally
partitioned training dataset.

Despite the reputation of neural networks as lacking interpretability, we
believe that there are a variety of potential theoretical connections between
the latent state learned by the LSTM of a SHRED model and other domains of
mathematics. In particular, there are parallels between the function of the
LSTM and Takens’ embedding, which states that with a sufficiently long time-
history, the dynamics of a time-delayed state variable are diffeomorphic to
the dynamics of the high-dimensional state space. While the LSTM does not
explicitly construct a time-delayed state, we conjecture that, once trained,
it does so implicitly in its processing of sensor measurements. The
feedforward network that follows the LSTM layer might then be thought of as
having learned a diffeomorphism between the analog of a time-delayed vector
and the high-dimensional state space.

We also believe that theoretical exploration of the the LSTM’s latent
dimension size is warranted. The plot in Fig. 1 show the MSE of
reconstructions of sea-surface temperature vs the number of nodes per hidden
layer and the SVD singular value spectrum. Both plots exhibit an clear elbow
at $\sim 5$ nodes per hidden layer and the $\sim 5$ singular value. In
traditional dimensionality reduction terms, the data exhibits a low-rank
structure and can be well approximated using the $\sim 5$ dominant modes
learned by the SVD. We conjecture that a similar relation might exist between
the minimum required LSTM state dimension and the underlying structure of the
data. In that case, the values of each entry in the latent state might be
understood in analogy to the entries in the matrix $V$ of the SVD $X=U\Sigma
V^{T}.$ The plot in Fig. 1 illustrates that these entries do appear to evolve
as some quasi-periodic phenomena.

Figure 2: (Top) Example flow fields to be reconstructed from the test set.
(Bottom) Histograms of the difference between ground truth and reconstruction
across all nodes and samples in the test set for immobile sensors (left) and
mobile sensors (right). Both distributions are approximately zero mean, but
the variance for mobile sensors is lower. (Center) Box plot of the MSE
evaluated across all samples in the test set for the 100 SHRED models with
mobile and immobile sensors. While the median performance is similar, there
are notably more outliers with poor performance in the case of immobile
sensors.

##  III Mobile sensing applications

To demonstrate the mobile SHRED architecture, three example problems are
selected: isotropic turbulence, sea-surface temperature and human
biomechanics. The physics of all three systems produce high-dimensional,
multiscale dynamics which is challenging to model in practice. However, we
will show in the figures that follow that SHRED provides exceptional
performance by mapping between arbitrary trajectories of sensor measurements
and the full state space.

###  III.1 Forced isotropic turbulent flow

The first dataset upon which we evaluate the performance of SHRED with mobile
sensors is a forced isotropic turbulent flow from the Johns Hopkins turbulence
database Li  _et al._ (2008). The data was generated via the pseudo-spectral
method with $1024^{3}$ nodes and utilized time steps of 0.0002 seconds. From
the simulation, we select a $50\times 50\times 50$ node volume with stride
length of 4 in all dimensions. We select 1,257 snapshots of this volume with
temporal spacing of 0.024 seconds. Of these 1,257 snapshots, the first 100 are
set aside to account for our use of a time-history of 100 sensor measurements,
900 of the remaining snapshots are then randomly selected for training while
the remainder are divided evenly between validation and test sets. In this
example, we do not consider a strict, temporal partition of training and test
data.

We train 100 SHRED models with one immobile, randomly placed sensors and 100
SHRED models with one mobile sensor following a modified random walk. Each
random walk is initialized from the center of the considered volume and every
three time-steps takes a one-node step in one of the three spatial dimensions,
each with equal probability. With high probability, SHRED models with either
mobile or immobile sensors accurately reconstruct the high-dimensional state
as measured by MSE across the test set. However, the box plot in Fig. 2 shows
that a greater proportion of models using mobile sensors achieve good
performance in comparison to the models using immobile sensors. The histograms
in Fig. 2 plot the deviation of the reconstruction from the ground truth at
each node and across all samples in the test set. The estimates for both
mobile and immobile sensors appear to be unbiased ($\mu=0.0015$ immobile,
$\mu=-0.0003$ mobile), but the variance of the mobile sensing reconstructions
is less ($\sigma^{2}=0.0998^{2}$ immobile, $\sigma^{2}=0.0617^{2}$ mobile).
Neither distribution is well approximated by a Gaussian, although the fit for
mobile sensors is better.

Sea surface temperature reconstruction error
---
Route |  | MSE
---
(random parition)
| MSE
---
(temporal partition)
Immobile sensors (3) | 0.0279 | 0.0460
Immobile sensor (1) | 0.0317 | 0.0481
Atlantic Ocean | 0.0358 | 0.0442
Antarctica | 0.0321 | 0.0516
USA West Coast | 0.0321 | 0.0499
Atlantic Ocean + Antarctica | 0.0395 | 0.1549
Atlantic Ocean + USA West Coast | 0.0440 | 0.1349
Antarctica + USA West Coast | 0.0328 | 0.1278
| Atlantic Ocean + Antarctica +
---
USA West Coast
0.0342 | 0.1345

Table 1: Mean-squared error for SHRED reconstructing sea-surface temperature
for both randomly and temporaly partitioned training/test/validation data.
Dynamic trajectories from mobile sensors in the Atlantic Ocean, Antarctic
Ocean, and along the USA’s West Coast were compared to one and three immobile
sensors for reconstruction the complex spatio-temporal sea-surface temperature
data Reynolds  _et al._ (2002). Figure 3: (Top) Example snapshots of global
sea-surface temperature to be reconstructed from either the (i) randomly or
(ii) temporally partitioned test set. (Bottom left) Histograms of the
difference between ground truth and reconstructed states across all grid
spaces and snapshots in the randomly-partitioned test set for immobile (1 and
3) and mobile sensors; dynamical trajectories of mobile sensors are one or
more combinations of a year-long circuit in (i) the Atlantic Ocean, Antarctica
Southern Ocean, and along the USA West Coast. (Bottom right) Histograms of the
differences between ground truth and reconstructed states across all grid
spaces and snapshots in the temporally-partitioned test set, again for
immobile and mobile sensors; sensor positions and trajectories were the same
as the randomly-partitioned dataset.

Marker-based Dataset results reconstructing kinematic states
---
Type | Input Sensors | MSE ($\pm$ s.d.) of Rotational Kinematics [∘]
SHRED | SDN | Linear
Individual |  | 3 random mobile sensors
---
(transverse-plane pelvis rotation angle,
medio-lateral pelvis position,
right hip adduction angle)
0.056 $\pm$ 0.018 | 0.220 $\pm$ 0.097 | 0.419 $\pm$ 0.243
| 3 non-random mobile sensors
---
(right hip flexion angle,
right knee flexion angle,
right ankle dorsiflexion angle)
0.064 $\pm$ 0.027 | 0.153 $\pm$ 0.097 | 0.514 $\pm$ 0.328
| 1 random mobile sensor
---
(medio-lateral pelvis position)
0.080 $\pm$ 0.031 | 1.033 $\pm$ 0.440 | 1.096 $\pm$ 0.470
| 1 non-random mobile sensor
---
(right ankle dorsiflexion angle)
0.087 $\pm$ 0.042 | 0.607 $\pm$ 0.244 | 1.083 $\pm$ 0.705
Population |  | 3 non-random mobile sensors
---
(right hip flexion angle,
right knee flexion angle,
right ankle dorsiflexion angle)
0.098 $\pm$ 0.062 | 0.311 $\pm$ 0.168 | 0.944 $\pm$ 0.847
| 1 non-random mobile sensor
---
(right ankle dorsiflexion angle)
0.121 $\pm$ 0.073 | 1.208 $\pm$ 0.822 | 1.843 $\pm$ 2.110

Table 2: Results from the human biomechanics dataset Rosenberg  _et al._
(2020). SHRED with random and non-random measurement trajectory inputs is
compared to a SDN and linear model for reconstructing kinematic states for
individual-specific models and population-based models. Overall, SHRED
successfully reconstructed kinematic states (rotational variables showed here)
with low mean-squared error, far outperforming the other architectures.
Translational kinematic variables (x, y, and z-direction pelvis position) were
excluded from these MSE calculations, as rotational and translational
variables have different units ([∘] and [$m$], respectively; see Figs. 4 and 5
for all kinematic results.

###  III.2 Sea-surface temperature

For our second dataset, we evaluate sea-surface temperature (SST) as reported
by NOAA Reynolds  _et al._ (2002). This dataset includes the weekly mean sea-
surface temperature from the years 1992 to 2019. Notably, NOAA produced the
SST data using optimum interpolation with in situ and satellite data. The data
comprises 1400 snapshots of a 180 by 360 grid; the exclusion of landmass
leaves 44,219 spatial locations corresponding to the sea surface. As in
Williams  _et al._ (2023), we denote SHRED’s input trajectory length to
include 52 temporal measurements, resulting in one year of measurement time
histories. The remaining 1348 snapshots are partitioned into training, test,
and validation sets, comprised of 1000, 174, and 174 snapshots respectively.
Regarding the dynamic trajectories of the mobile sensors, we consider three
distinct moving sensors that complete a circuit over a year span: (i) an
Atlantic Ocean route that travels from the east coast of the USA towards
Europe and back, (ii) an out-and-back route near Antarctica Southern Ocean,
and (iii) a route that follows the West Coast of the USA.

In Table 1, we compare how well each route – as well as their combinations –
reconstruct SST, relative to one and three immobile sensors, for both random
and temporal data partitions. For the random partition of
training/test/validation data – an interpolation problem – we observe
comparable mean-squared reconstruction errors across all mobile sensor routes
relative to the immobile sensor reconstruction errors. For the temporal
partition – an extrapolation problem – we see comparable errors for the single
dynamic trajectory routes relative to the immobile sensors. The 2+ route
combinations demonstrate much higher reconstruction errors; such results could
be a potential indicator of overfitting, or even demonstrative of the
challenges of extrapolating complex spatio-temporal data.

We additionally evaluated the mean and variance in the mean-squared
reconstruction errors across routes and data-partitions, as seen in Fig. 3.
For the randomly-partitioned data of global sea-surface temperature, immobile
sensors are more biased than mobile sensors ($\mu=-0.0306$ 3 immobile sensors,
$\mu=-0.0577$ 1 immobile sensor; across all mobile sensor routes, minimum
$\mu=-0.0004$, maximum $\mu=0.0446$, absolute mean $\mu=0.0241$). However, the
variance of the mobile versus immobile sensors are more variable
($\sigma^{2}=0.5080^{2}$ 3 immobile sensors, $\sigma^{2}=0.5958^{2}$ 1
immobile sensor; across all mobile sensor routes, minimum $\sigma^{2}=0.5410$,
maximum $\sigma^{2}=1.4131$, mean $\sigma^{2}=0.7756$). For the temporally-
partitioned data, immobile sensors are less biased than mobile sensors
($\mu=-0.0291$ 3 immobile sensors, $\mu=-0.0132$ 1 immobile sensor; across all
mobile sensor routes, minimum $\mu=-0.0896$, maximum $\mu=0.1907$, absolute
mean $\mu=0.1406$). We see comparable variance for the single dynamic
trajectory routes relative to the immobile sensors; the 2+ route combinations
have much higher error variance; ($\sigma^{2}=0.7881^{2}$ 3 immobile sensors,
$\sigma^{2}=1.0108^{2}$ 1 immobile sensor; across all mobile sensor routes,
minimum $\sigma^{2}=0.8281$, maximum $\sigma^{2}=2.8122$, mean
$\sigma^{2}=1.7959$).

Figure 4: Visualization of individual-specific results with non-random sensor
inputs for reconstructing human biomechanics Rosenberg  _et al._ (2020); see
Table 2 for additional individual-specific results using random sensor inputs.
A unique mapping was trained for each individual to transform their sparse set
of sensor measurements to their full measurement data; each individual’s
mapping was evaluated by reconstructing their held-out test data. Performance
was calculated using mean-square error (MSE); box plots show the aggregate MSE
for all individuals and across kinematic states. Rotational and translational
kinematic states are displayed separately, as rotational states use units of
degrees while translational states use units of meters. Three modeling
paradigms were trained to learn a mapping: linear regression (orange), a
shallow decoder network (SDN) (green), and a shallow recurrent demcoder
network (SHRED) (blue). Note the independent y-axes and respective scales for
each graph. (Top) Three dynamic trajectories were purposefully chosen as
sensor inputs to the reconstruction mapping: right hip flexion angle
(degrees), right knee flexion angle (degrees), and right ankle dorsiflexion
angle (degrees). (Bottom) One dynamic trajectory was purposefully chosen as
sensor inputs to the reconstruction mapping: right ankle dorsiflexion angle
(degrees). Figure 5: Visualization of population-based results for
reconstructing human biomechanics Rosenberg  _et al._ (2020). The population
models trained using 11 of the 12 individuals’ data and were evaluated by
reconstructing the held-out subject’s full measurement data from their sparse
input sensor(s); performance was calculated by using mean-square error (MSE).
Rotational and translational kinematic states are displayed separately, as
rotational states use units of degrees while translational states use units of
meters. Three modeling paradigms were trained to learn the mapping: linear
regression (orange), a shallow decoder neural network (green), and a shallow
recurrent neural network (blue). Note the independent y-axes and respective
scales for each graph. (Top) Three dynamic trajectories were purposefully
chosen as sensor inputs to the reconstruction mapping: right hip flexion angle
(degrees), right knee flexion angle (degrees), and right ankle dorsiflexion
angle (degrees). (Bottom) One dynamic trajectory was purposefully chosen as
sensor inputs to the reconstruction mapping: right ankle dorsiflexion angle
(degrees).

###  III.3 Human biomechanics

The final dataset for which we leverage mobile sensor trajectories with SHRED
is for capturing human biomechanics. Human motion tracking and analysis is
essential for monitoring disease progression, guiding rehabilitation
treatment, evaluating sports performance, and informing assistive device
design. We use an open-source dataset that captures nondisabled human
biomechanics during steady, rhythmic walking; such an approach can be
generalized to other motion tracking such as robotic manipulation or computer
animation. We use kinematics from 12 nondisabled adults (six female/six male;
age 23.9 $\pm$ 1.8 years; height = 1.69 $\pm$ 0.10m; mass = 66.5 $\pm$ 11.7
kg) Rosenberg  _et al._ (2020). Participants walked at a self-selected speed
(speed = 1.36 $\pm$ 0.11 m/s) on a split-belt instrumented treadmill (Bertec
Corp., Columbus, USA) for six minutes. While kinematics can be captured with
multiple modalities, gold-standard marker data was recorded using a 10-camera
optical motion capture system (Qualisys AB, Gothenburg, SE). Kinematics were
estimated from marker data using the Inverse Kinematics algorithm in OpenSim
3.3 with a dynamically-constrained 19 degree-of-freedom skeletal model scaled
to each participant Delp  _et al._ (2007); Hicks  _et al._ (2015). We
evaluated the 18 kinematic states previously calculated in the open-source
dataset which included three translational degrees of freedom at the pelvis,
pelvis tilt/list/rotation, lumbar extension/bending, bilateral hip
flexion/adduction/rotation, bilateral knee flexion, and bilateral ankle
dorsiflexion. Kinematics were low-pass filtered at 6Hz using a fourth-order
Butterworth filter (Rosenberg  _et al._ , 2020). An important distinction with
this human biomechanics data, as opposed to the isotropic turbulence or sea-
surface temperature, is that marker data recorded during motion capture is in
a global coordinate frame, but through the standard procedures for marker data
processing are transformed into kinematics in a relative frame of reference to
the human subject.

Because this kinematic dataset cannot be compared to SHRED with immobile
sensors, as the markers must move with the human, we compare the performance
of mobile sensor trajectories and SHRED with that of a shallow decoder network
(SDN) (i.e., SHRED without leveraging time histories) and a linear model.
Therefore, the aim of this example was to evaluate each modeling
architecture’s ability to learn a mapping between a sparse set of mobile
sensors to reconstruct the full set of 18 kinematic states. We evaluate two
conditions: (1) individual-specific models, in which a model is developed for
each participant to reconstruct their respective set of full kinematics, and
(2) population-based models, which uses data from all participants to
reconstruct the full set of kinematics for an unseen individual, thus enabling
rapid generalization to data outside the training set.

Overall, SHRED successfully maps between a sparse set of measurements to
reconstruct the full biomechanical states for both the individual and
population models. Additionally, SHRED far outperforms the linear and SDN
architectures. While more input sensors generally increases reconstruction
accuracy, performance with as little as 0.068 and 0.048 degrees (rotational
variables) mean-squared error can be achieved with one and three sensors with
SHRED, well below the accuracy and repeatability recommended even for clinical
gait analyses ($2\deg$ in sagittal-plane and $5\deg$ in frontal-plane) Schmitz
_et al._ (2014); Hicks  _et al._ (2015). Table 2 provides a summary of the
mean-squared errors for the rotational kinematic variables across modeling
architectures and conditions. We take a closer look at each condition’s
results in the following sections:

####  III.3.1 Individual-specific biomechanics models

For the individual-specific biomechanics models, a mapping is trained for each
of the 12 participants using a linear, SDN and SHRED architecture. Four
combinations of dynamic trajectory inputs are tested: (1) three randomly
chosen kinematic states (transverse-plane pelvis rotation angle, medio-lateral
pelvis position, right hip adduction angle), (2) three non-randomly chosen
kinematic states (right hip flexion angle, right knee flexion angle, right
ankle dorsiflexion angle), (3) one randomly chosen kinematic state (medio-
lateral pelvis position), and (4) one non-randomly chosen kinematic state
(right ankle dorsiflexion angle). The non-randomly chosen states represent the
kinematics most commonly used for biomechanical assessments, and which
generally have the greatest accuracy and inter-session reliability (i.e,
sagittal-plane kinematics) Schmitz  _et al._ (2014).

Fig. 4 shows the performance of SHRED relative to the SDN and linear model for
reconstructing human motion across individuals. SHRED reconstructs the full
set of 18 biomechanical variables more accurately that the shallow decoder and
linear models, regardless of the number or choice of input sensors. Using
three mobile sensors, SHRED has errors less than $0.064\pm 0.027$ degrees
(rotational variables), as much as $4.0$x more accurate on average than the
SDN and $8.0$x more accurate than the linear model. Using a single mobile
sensor, SHRED has errors less than $0.087\pm 0.042$ degrees, as much as
$12.8$x more accurate on average than the SDN and $13.6$x more accurate than
the linear model.

####  III.3.2 Population-based biomechanics models

For the population-based biomechanics models, a reconstruction mapping was
trained using 11 of the 12 participants’ kinematic data, for a one-subject
hold-out test using each architecture. Two combinations of dynamic trajectory
inputs were tested: (1) three non-randomly chosen kinematic states (right hip
flexion angle, right knee flexion angle, right ankle dorsiflexion angle), and
(2) one non-randomly chosen kinematic state (right ankle dorsiflexion angle).

Fig. 5 shows the performance of SHRED relative to the SDN and linear model for
reconstructing human motion within a population. SHRED far outperforms the SDN
and linear model regardless of number or choice of input sensors. Using three
mobile sensors, SHRED has errors less than $0.098\pm 0.062$ degrees, which was
on average $3.3$x more accurate than the SDN and $9.6$x more accurate than the
linear model. Using a single mobile sensor, SHRED has errors less than
$0.121\pm 0.073$ degrees, which was on average $9.8$x more accurate than the
SDN and $14.8$x more accurate than the linear model. These results demonstrate
that SHRED can enable rapid generalization (parameterization of dynamics) for
data outside the training set, which may be especially useful in human
biomechanics where population-level models can be used to estimate parameters
for a new individual (e.g., a new patient visits a clinic for the first time).

##  IV Discussion & Conclusion

Sensing continues to be one of the most important tasks for science and
engineering across disciplines. Machine learning methods are offering
exceptional new paradigms for maximally exploiting sensor information for the
diversity of tasks required of sensors, including reconstruction, forecasting,
model discovery, control, and uncertainty quantification. Although there is a
well established theory for sensing and sensor placement using immobile
sensors, mobile sensing has remained exceptionally challenging with only
recent and limited advancements made for leveraging the dynamical trajectories
of the sensor. The mobile SHRED architecture advocated here provides a new
paradigm for sensing whereby the time-history of the sensor is used to encode
global information of the measured high-dimensional state space. SHRED not
only allows for the embedding of the multiscale physics into a compact and
low-dimensional latent space, it also provides through its decoder network a
mapping from the minimal mobile sensors to full state estimates. Thus the
SHRED architecture provides a nonlinear generalization of classical low-rank
and linear embeddings such as the singular value decomposition.

The performance of SHRED is demonstrated on three challenging data sets which
are all characterized by complex, nonlinear, and multiscale interactions. Thus
the data considered represents some of the most challenging environments for
sensing. Additionally, such applications often require the mobile sensing.
Given the limited work on mobile sensing to date, along with the performance
demonstrated by the method, the mobile SHRED architecture represents a
significant innovation in mobile sensing technologies. Not only does the
architecture reduce the variance in state estimation in comparison to immobile
sensors, the decoder requires less data to train than typical networks due to
its shallow structure. Empirically it is also observed that the mobile SHRED
architecture is robust to hyper-parameter tuning, requiring little effort to
train high performing models. The advantageous of such an architecture provide
an application agnostic scheme which has significant potential for broad usage
in science and engineering.

The mobile SHRED architecture works because the time history of the sensors
embed global information of the system measured. For spatio-temporal data, for
instance, the dynamics typically are governed by partial differential
equations whose spatial derivative estimates rely on neighboring data points.
These neighbors in turn rely on their neighbors and so forth. Thus a local
spatial position is coupled to neighboring positions and ultimately globally
by spatial derivatives. The sensor trajectories learn a representation of the
spatio-temporal field via this global coupling. Of course, such an
architecture would break down if there are regions in the dynamics that are
statistically independent from each other, i.e. there are isolated dynamics
that are not coupled in any way. There can also be “dead spots” in the
measurement space, but the moving sensors can then typically move to new
locations where information is captured. Thus with high-probability, a random
sensor trajectory can be used to encode global information and estimate the
full state space.

Finally, for the examples shown, there is clearly a requirement in the
training data to have both the sensor trajectory data as well as the
corresponding high-dimensional state space data. In many applications, there
is no access to the high-dimensional state space. Thus proxy data, in the form
of simulation data, can be used as a surrogate for the real data. As long as
the proxy data is representative of the true data, it can be used to train the
SHRED architecture. This gives a hybrid approach to training the neural
network whereby simulation data can be used in partnership with real data to
build a sensing model for producing otherwise inaccessible full state
estimates.

## Acknowledgements

JNK acknowledges support from the Air Force Office of Scientific Research
(FA9550-19-1-0386). The authors also acknowledge support in part by the US
National Science Foundation (NSF) AI Institute for Dynamical Systems
(dynamicsai.org), grant 2112085.

## References

  * Manohar  _et al._ (2018) K. Manohar, B. W. Brunton, J. N. Kutz, and S. L. Brunton, IEEE Control Systems Magazine 38, 63 (2018).
  * Gardner  _et al._ (2006) J. P. Gardner, J. C. Mather, M. Clampin, R. Doyon, M. A. Greenhouse, H. B. Hammel, J. B. Hutchings, P. Jakobsen, S. J. Lilly, K. S. Long, _et al._ , Space Science Reviews 123, 485 (2006).
  * Stevenson and Kording (2011) I. H. Stevenson and K. P. Kording, Nature neuroscience 14, 139 (2011).
  * Medsker and Jain (2001) L. R. Medsker and L. Jain, Design and Applications 5, 2 (2001).
  * Salehinejad  _et al._ (2017) H. Salehinejad, S. Sankar, J. Barfett, E. Colak, and S. Valaee, arXiv preprint arXiv:1801.01078 (2017).
  * Zaremba  _et al._ (2014) W. Zaremba, I. Sutskever, and O. Vinyals, arXiv preprint arXiv:1409.2329 (2014).
  * Erichson  _et al._ (2020) N. B. Erichson, L. Mathelin, Z. Yao, S. L. Brunton, M. W. Mahoney, and J. N. Kutz, Proceedings of the Royal Society A 476, 20200097 (2020).
  * Williams  _et al._ (2022) J. Williams, O. Zahn, and J. N. Kutz, arXiv preprint arXiv:2202.05330 (2022).
  * Williams  _et al._ (2023) J. P. Williams, O. Zahn, and J. N. Kutz, arXiv preprint arXiv:2301.12011 (2023).
  * Law  _et al._ (2015) K. Law, A. Stuart, and K. Zygalakis,  _Data Assimilation: A Mathematical Introduction_ , 1st ed. (Springer Cham, 2015).
  * Kutz (2013) J. N. Kutz,  _Data-driven modeling & scientific computation: methods for complex systems & big data_ (OUP Oxford, 2013).
  * Brunton and Kutz (2019) S. L. Brunton and J. N. Kutz,  _Data-driven science and engineering: Machine learning, dynamical systems, and control_ (Cambridge University Press, 2019).
  * Everson and Sirovich (1995) R. Everson and L. Sirovich, JOSA A 12, 1657 (1995), publisher: Optica Publishing Group.
  * Drmač and Gugercin (2016) Z. Drmač and S. Gugercin, SIAM Journal on Scientific Computing 38, A631 (2016), publisher: Society for Industrial and Applied Mathematics.
  * Barrault  _et al._ (2004) M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera, Comptes Rendus Mathematique 339, 667 (2004).
  * Chaturantabut and Sorensen (2010) S. Chaturantabut and D. C. Sorensen, SIAM Journal on Scientific Computing 32, 2737 (2010), publisher: Society for Industrial and Applied Mathematics.
  * Boyd and Vandenberghe (2004) S. P. Boyd and L. Vandenberghe,  _Convex optimization_ (Cambridge university press, 2004).
  * Joshi and Boyd (2009) S. Joshi and S. Boyd, IEEE Transactions on Signal Processing 57, 451 (2009), conference Name: IEEE Transactions on Signal Processing.
  * Caselton and Zidek (1984) W. F. Caselton and J. V. Zidek, Statistics & Probability Letters 2, 223 (1984).
  * Krause  _et al._ (2008) A. Krause, A. Singh, and C. Guestrin, The Journal of Machine Learning Research 9, 235 (2008).
  * Lindley (1956) D. V. Lindley, The Annals of Mathematical Statistics 27, 986 (1956), publisher: Institute of Mathematical Statistics.
  * Sebastiani and Wynn (2000) P. Sebastiani and H. P. Wynn, Journal of the Royal Statistical Society. Series B (Statistical Methodology) 62, 145 (2000), publisher: [Royal Statistical Society, Wiley].
  * Paninski (2005) L. Paninski, Neural Computation 17, 1480 (2005).
  * Gunnarson  _et al._ (2021) P. Gunnarson, I. Mandralis, G. Novati, P. Koumoutsakos, and J. O. Dabiri, Nature communications 12, 7143 (2021).
  * Krishna  _et al._ (2022) K. Krishna, Z. Song, and S. L. Brunton, Proceedings of the Royal Society A 478, 20210255 (2022).
  * Buzzicotti  _et al._ (2020) M. Buzzicotti, L. Biferale, F. Bonaccorso, P. Clark di Leoni, and K. Gustavsson, in  _International Conference of the Italian Association for Artificial Intelligence_ (Springer, 2020) pp. 223–234.
  * Biferale  _et al._ (2019) L. Biferale, F. Bonaccorso, M. Buzzicotti, P. Clark Di Leoni, and K. Gustavsson, Chaos: An Interdisciplinary Journal of Nonlinear Science 29 (2019).
  * Madridano  _et al._ (2021) A. Madridano, A. Al-Kaff, D. Martín, and A. De La Escalera, Expert Systems with Applications 173, 114660 (2021).
  * Shriwastav  _et al._ (2022) S. Shriwastav, G. Snyder, and Z. Song, in  _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ (IEEE, 2022) pp. 11910–11915.
  * Lynch  _et al._ (2008) K. M. Lynch, I. B. Schwartz, P. Yang, and R. A. Freeman, IEEE transactions on robotics 24, 710 (2008).
  * Leonard  _et al._ (2007) N. E. Leonard, D. A. Paley, F. Lekien, R. Sepulchre, D. M. Fratantoni, and R. E. Davis, Proceedings of the IEEE 95, 48 (2007).
  * DeVries  _et al._ (2013) L. DeVries, S. J. Majumdar, and D. A. Paley, Journal of intelligent & robotic systems 70, 527 (2013).
  * Ogren  _et al._ (2004) P. Ogren, E. Fiorelli, and N. E. Leonard, IEEE Transactions on Automatic control 49, 1292 (2004).
  * Zhang and Leonard (2010) F. Zhang and N. E. Leonard, IEEE Transactions on Automatic Control 55, 650 (2010).
  * Paley and Wolek (2020) D. A. Paley and A. Wolek, Annual Review of Control, Robotics, and Autonomous Systems 3, 91 (2020).
  * Peng  _et al._ (2014) L. Peng, D. Lipinski, and K. Mohseni, Journal of Intelligent & Robotic Systems 74, 421 (2014).
  * Mei  _et al._ (2022) J. Mei, S. L. Brunton, and J. N. Kutz, arXiv preprint arXiv:2212.08280 (2022).
  * LeCun  _et al._ (2015) Y. LeCun, Y. Bengio, and G. Hinton, nature 521, 436 (2015).
  * Hornik  _et al._ (1989) K. Hornik, M. Stinchcombe, and H. White, Neural networks 2, 359 (1989).
  * Hochreiter and Schmidhuber (1997) S. Hochreiter and J. Schmidhuber, Neural Computation 9, 1735 (1997).
  * Kingma and Ba (2017) D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,”  (2017), arXiv:1412.6980 [cs].
  * Li  _et al._ (2008) Y. Li, E. Perlman, M. Wan, Y. Yang, C. Meneveau, R. Burns, S. Chen, A. Szalay, and G. Eyink, Journal of Turbulence 9, N31 (2008), publisher: Taylor & Francis _eprint: https://doi.org/10.1080/14685240802376389.
  * Reynolds  _et al._ (2002) R. W. Reynolds, N. A. Rayner, T. M. Smith, D. C. Stokes, and W. Wang, Journal of climate 15, 1609 (2002).
  * Rosenberg  _et al._ (2020) M. C. Rosenberg, B. S. Banjanin, S. A. Burden, and K. M. Steele, Journal of the Royal Society Interface 17, 20200487 (2020).
  * Delp  _et al._ (2007) S. L. Delp, F. C. Anderson, A. S. Arnold, P. Loan, A. Habib, C. T. John, E. Guendelman, and D. G. Thelen, IEEE transactions on biomedical engineering 54, 1940 (2007).
  * Hicks  _et al._ (2015) J. L. Hicks, T. K. Uchida, A. Seth, A. Rajagopal, and S. L. Delp, Journal of biomechanical engineering 137, 020905 (2015).
  * Schmitz  _et al._ (2014) A. Schmitz, M. Ye, R. Shapiro, R. Yang, and B. Noehren, Journal of biomechanics 47, 587 (2014).