# Robust Machine Learning Modeling for Predictive Control Using Lipschitz-
Constrained Neural Networks

Wallace Tan Gian Yion  Zhe Wu  Department of Chemical and Biomolecular
Engineering, National University of Singapore, 117585, Singapore.

###### Abstract

Neural networks (NNs) have emerged as a state-of-the-art method for modeling
nonlinear systems in model predictive control (MPC). However, the robustness
of NNs, in terms of sensitivity to small input perturbations, remains a
critical challenge for practical applications. To address this, we develop
Lipschitz-Constrained Neural Networks (LCNNs) for modeling nonlinear systems
and derive rigorous theoretical results to demonstrate their effectiveness in
approximating Lipschitz functions, reducing input sensitivity, and preventing
over-fitting. Specifically, we first prove a universal approximation theorem
to show that LCNNs using SpectralDense layers can approximate any 1-Lipschitz
target function. Then, we prove a probabilistic generalization error bound for
LCNNs using SpectralDense layers by using their empirical Rademacher
complexity. Finally, the LCNNs are incorporated into the MPC scheme, and a
chemical process example is utilized to show that LCNN-based MPC outperforms
MPC using conventional feedforward NNs in the presence of training data noise.

###### keywords:

Lipschitz-Constrained Neural Networks; Robust Machine Learning Model;
Generalization Error; Model Predictive Control; Neural Network Sensitivity;
Over-fitting

††journal: Elsevier

##  1 Introduction

Model Predictive Control (MPC) is an advanced optimization-based control
strategy for various chemical engineering processes, such as batch
crystallization processes (Kwon et al. (2013, 2014)) and continuous stirred
tank reactors (CSTRs) (Chen et al. (1995); Wu (2001)). Machine learning (ML)
techniques such as artificial neural networks (ANN) have been utilized to
develop prediction models that will be incorporated into the design of MPC.
Among various ML-based MPC schemes, an accurate and robust process model for
prediction has always been one of the key components that ensures desired
closed-loop performance.

Despite the success of ANNs in modeling complex nonlinear systems, one
prominent issue is that they could potentially be sensitive to small
perturbations in input features. Such sensitivity issues arise naturally when
a small change in the input (e.g., due to data noise, perturbation, or
artificially generated adversarial inputs) could result in a drastic change in
the output. For example, in classification problems, adversarial input
perturbations have been shown to lead to a large variation in NN output,
thereby leading to misclassified results (Szegedy et al. (2013)).
Additionally, Balda et al. (2019) proposed a novel approach for constructing
adversarial examples for regression problems using perturbation analysis for
the underlying learning algorithm for the neural network. Since the lack of
robustness of NNs could adversely affect the prediction accuracy in ML-based
MPC, it is important to address the sensitivity issue for the implementation
of NNs in performance-critical applications.

To mitigate this issue, adversarial training has been adopted as one of the
most effective approaches to train NNs against adversarial examples (Szegedy
et al. (2013)). For example, Shaham et al. (2015) proposed a robust
optimization framework with a different loss function in the optimization
process that aims to search for adversarial examples near the vicinity of each
training data point. This approach has been empirically shown to be effective
against adversarial attacks by both Shaham et al. (2015) and Madry et al.
(2017). However, there is a lack of a provably performance guarantee on the
input sensitivity of neural networks for this approach.

Therefore, a type of neural network with a fixed Lipschitz constant termed
Lipschitz-Constrained Neural Networks (LCNNs) has received an emerging amount
of interest in recent years (Baldi and Sadowski (2013); Anil et al. (2019)).
One immediate way to control the Lipschitz constant of a neural network is to
bound the norms of the weight matrices, and use activation functions with
bounded derivatives (e.g., ReLU activation function). However, it is
demonstrated in Anil et al. (2019) that this method substantially limits the
model capacity of the NNs if component-wise activation functions are used. A
recent breakthrough using a special activation function termed GroupSort (Anil
et al. (2019)) significantly increases the expressive power of the NNs.
Specifically, the GroupSort activation function enables NNs with norm-
constrained GroupSort architectures to serve as universal approximators of
1-Lipschitz functions. Using the restricted Stone-Weierstrass theorem, Anil et
al. (2019) proved the universal approximation theorem for GroupSort
feedforward neural networks with appropriate assumptions on the norms of the
weight matrices. However, since the proof in Anil et al. (2019) is based on
the $\infty$-norm of matrices, the way to demonstrate the universal
approximation property of LCNNs using the spectral norm (i.e., $\ell_{2}$
matrix norm) remains an open question.

The second prominent issue in the development of ANNs is over-fitting, where
the networks perform very well on the training data but fail to predict
accurately on the test data, which results in a high generalization error. One
possible reason for over-fitting is that the training data has noise that
negatively impacts learning performance (Ying (2019)). Additionally, over-
fitting occurs when there is insufficient training data, or when there is a
high hypothesis complexity, in terms of large weights, a large number of
neurons, and extremely deep architectures. Therefore, designing neural network
architectures that are less prone to over-fitting is a pertinent issue in
supervised machine learning. Sabiri et al. (2022) provides an overview of
popular solutions to prevent over-fitting. For example, one of the most common
solutions to over-fitting is regularization, such as $\ell_{1}$ or $\ell_{2}$
regularization (e.g. Moore and DeNero (2011); Cortes et al. (2012)), which
implements the size of the weights as a soft constraint. Other popular
solutions include dropout (e.g. Baldi and Sadowski (2014); Srivastava et al.
(2014)), where certain neurons are dropped out during training time with a
certain specified probability, and early stopping, where the training is
stopped using a predefined predicate, usually when the validation error
reaches a minimum (e.g. Baldi and Sadowski (2013)). For example, in our
previous work Wu et al. (2021a), the Monte Carlo dropout technique was
utilized in the development of NNs to mitigate the impact of data noise and
reduce over-fitting. In addition to the above solutions, LCNNs have been
demonstrated to be able to efficiently avoid over-fitting by constraining the
Lipschitz constant of a network. However, at this stage, a fundamental
understanding of the capability of LCNNs in reducing over-fitting in terms of
the generalization ability of LCNNs over the underlying data distribution is
missing.

Motivated by the above considerations, in this work, we incorporate LCNNs
using SpectralDense layers in MPC and demonstrate that the LCNNs can
effectively resolve the two aforementioned issues: sensitivity to input
perturbations and over-fitting in the presence of noise. Rigorous theoretical
results are developed to demonstrate that LCNNs are provably robust against
input perturbations because of their low Lipschitz constant and provably
robust against over-fitting due to their lowered hypothesis complexity, i.e.,
low Rademacher complexity. The rest of this article is organized as follows.
In Section 2, the nonlinear systems that are considered and the application
method of FNNs in MPC are first presented. In Section 3, we present the
formulation of LCNNs using SpectralDense layers, followed by a discussion on
their improved robustness against input perturbations. In Section 4, we prove
the universal approximation theorem for 1-Lipschitz continuous functions for
LCNNs using SpectralDense layers. In Section 5, we develop a probabilistic
generalization error bound for LCNNs to show that LCNNs using SpectralDense
layers can effectively prevent over-fitting. This is done by computing an
upper bound on the empirical Rademacher complexity method (ERC) of the
function class represented by LCNNs using SpectralDense layers. Finally, in
Section 6, we carry out a simulation study of a benchmark chemical reactor
example, where we will exhibit the superiority of LCNNs over conventional FNNs
with dense layers in the presence of noisy training data.

##  2 Preliminaries

###  2.1 Notations

$\lVert W\rVert_{F}$ and $\lVert W\rVert_{2}$ denote the Frobenius norm and
the spectral norm of a matrix $W\in\mathbb{R}^{n\times m}$ respectively.
$\mathbb{R}^{\geq 0}$ denotes the set of all nonnegative real numbers. A
function $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is continuously differentiable if
and only if it is differentiable and the Jacobian of $f$, denoted by $J_{f}$,
is continuous. Given a vector $x\in\mathbb{R}^{n}$, let $\lVert x\rVert$
denote the Euclidean norm of $x$. A metric space is a set $X$ equipped with a
metric function $d_{X}:X\times X\to\mathbb{R}^{\geq 0}$ such that 1)
$d_{X}(x,y)=0$ if and only if $x=y$, and 2) for all $x,y,z\in X$, the triangle
inequality $d_{X}(x,z)\leq d_{X}(x,y)+d_{X}(y,z)$ holds. We denote a metric
space as an ordered pair $(X,d_{X})$. Given an event $A$, we denote
$\mathbb{P}(A)$ to be its probability. Given a random variable $X$, we denote
$\mathbb{E}[X]$ to be its expectation.

###  2.2 Class of Systems

The nonlinear systems that are considered in this article can be represented
by the following ordinary differential equation (ODE):

$\dot{x}=F(x,u):=f(x)+g(x)u$ (1)

Here $x\in\mathbb{R}^{n}$ is the current state vector, $u\in\mathbb{R}^{m}$ is
the control vector, and $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$,
$g:\mathbb{R}^{n}\to\mathbb{R}^{n\times m}$ are continuously differentiable
matrix-valued functions. We also assume that $f(0)=0$, so that the origin
$(x,u)=(0,0)$ is an equilibrium point.

We also assume that there is a Lyapunov function $V:D\to\mathbb{R}^{\geq 0}$
that is continuously differentiable and is equipped with a controller
$\Phi:D\to U$, such that the origin $(x,u)=(0,0)$ is an equilibrium point that
is exponentially (closed-loop) stable. Here $D$ and $U$ are compact subsets of
$\mathbb{R}^{n}$ and $\mathbb{R}^{m}$, respectively, that contain an open set
surrounding the origin. In addition, the stability region for this controller
$\Phi(x)$ is then taken to be a sublevel set of $V$, i.e.,
$\Omega_{\rho}:=\\{x\leavevmode\nobreak\ |\,V(x)\leq\rho\\}$ with $\rho$
positive and $\Omega_{\rho}\subset D$. Additionally, since $f,g$ are
continuously differentiable, the following inequalities can be readily derived
for any $x,x^{\prime}\in D$, $u\in U$ and some constants $K_{F}$, $L_{x}$:

$\lVert F(x,u)\rVert\leq K_{F}$ (2a) $\lVert F(x^{\prime},u)-F(x,u)\rVert\leq
L_{x}\lVert x^{\prime}-x\rVert$ (2b)

###  2.3 Feedforward Neural Network (FNN)

This subsection gives a short summary of the development of FNNs for the
nonlinear systems represented by Eq. 1. Specifically, since the FNN is
developed as the prediction model for model predictive controllers, we
consider the FNNs that are built to capture the nonlinear dynamics of Eq. 1 ,
whereby the control input actions are applied using the sample-and-hold
method. This means that given an initial state $x_{0}$, the control action $u$
applied is constant throughout the entire period of the sampling time
$\Delta>0$. Suppose that the system state is currently at $x(0)=x_{0}$. From
Eq. 2b and the Picard-Lindelöf Theorem for ODEs, there exists a unique state
trajectory $x(t)$ such that

$x(t)=x_{0}+\int_{0}^{t}F(x(s),u)ds$ (3)

We can then define
$\tilde{F}:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}$ by
$\tilde{F}(x_{0},u)=x(\Delta)$, which is the function that takes the current
state and the control action as inputs and predicts the state at $\Delta$ time
later. This function can then be approximated by an FNN, such as an LCNN, and
the approximation is denoted by $\tilde{F}_{nn}$. In order to develop the
neural networks, open-loop simulations of the ODE system represented by Eq. 1
using varying control actions will be conducted to generate the required
dataset. Specifically, we perform a sweep of all possible initial states
$x_{0}\in\Omega_{\rho}$ and control actions $u\in U$, and use the forward
Euler method with integration time step $h_{c}\ll\Delta$, to deduce the value
of $\tilde{F}(x_{0},u)$, that is, to deduce the state at $t=\Delta$. The
training dataset consists of all such possible pairs $(x_{0},u)$ (FNN inputs)
$\tilde{F}(x_{0},u)$ (FNN outputs), where $\tilde{F}(x_{0},u)$ is the actual
future state. Therefore, the two functions $x_{t+1}:=\tilde{F}(x_{t},u_{t})$
and $x_{t+1}:=\tilde{F}_{nn}(x_{t},u_{t})$ represent define two distinct
nonlinear discrete-time systems respectively, where $x_{t+1}$ denotes the
state at $t+\Delta$.

Since $\tilde{F}_{nn}$ will be the function used in the MPC optimization
algorithm, ensuring that $\tilde{F}_{nn}$ is an accurate approximation of
$\tilde{F}$ is necessary so that the neural network captures the nonlinear
dynamics well. In general, the FNN model $\tilde{F}_{nn}$ should be developed
with sufficient training data and an appropriate architecture in terms of the
number of neurons and layers in order to achieve the desired prediction
accuracy on both training and test sets. However, in the presence of
insufficient training data or noisy data, over-fitting might occur, leading to
a model that performs poorly on the test data set and generalizes poorly.
Also, the developed FNN model for $\tilde{F}_{nn}$ should not be overtly
sensitive to input perturbations (that is, it should not have overtly large
gradients) in order to ensure that $\tilde{F}_{nn}$ can generalize well to
other data points outside the training dataset but within the desired domain
$D\times U$. Therefore, to address the issues of over-fitting and sensitivity,
we will develop LCNNs for the nonlinear system of Eq. 1 in this work, and show
that LCNNs can overcome the limitations of conventional FNNs with dense layers
by lowering sensitivity and preventing over-fitting.

##  3 Lipschitz-Constrained Neural Network Models Using SpectralDense Layers

In this section, the architecture of LCNNs using SpectralDense layers will
first be introduced, followed by a discussion on the reduced sensitivity of
LCNN to input perturbations as compared to conventional FNNs. First, we begin
with an important definition.

######  Definition 1.

A function $f:X\to Y$ where $X\subset\mathbb{R}^{n}$ and
$Y\subset\mathbb{R}^{m}$ is Lipschitz continuous with Lipschitz constant $L$
(or $L$-Lipschitz) if $\forall x,y\in X$, one has

$\lVert f(x)-f(y)\rVert\leq L\cdot\lVert x-y\rVert$

It is readily shown that if $f$ is $L$-Lipschitz continuous, given a small
perturbation to the input, the output $f(x)$ changes by at most $L$ times the
magnitude of that perturbation. As a result, as long as the Lipschitz constant
of a neural network is constrained to be a small value, it is less sensitive
with respect to input perturbations. In the next subsection, we demonstrate
that LCNNs using SpectralDense layers have a constrained small Lipschitz
constant, where each of the SpectralDense layers has a Lipschitz constant of
1.

###  3.1 SpectralDense layers

The mathematical definition of the SpectralDense layers used to construct an
LCNN is first presented. Recall that by the singular value decomposition
(SVD), for any $W\in\mathbb{R}^{m\times n}$, there exist orthogonal matrices
$U\in\mathbb{R}^{m\times m}$, $V\in\mathbb{R}^{n\times n}$, and a rectangular
diagonal matrix $D\in\mathbb{R}^{m\times n}$ with positive entries such that
$W=UDV^{T}$. First, we recall the definition of a conventional dense layer:

######  Definition 2.

A dense layer is a function $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ of the form

$f:x\to\sigma(Wx+b)$

where $\sigma\in\mathbb{R}^{m}\to\mathbb{R}^{m}$ is an activation function,
$W\in\mathbb{R}^{m\times n}$ is a weight matrix, and $b\in\mathbb{R}^{n}$ is a
bias term. A dense layer is a layer that is deeply connected with its
preceding layer (i.e., the neurons of the layer are connected to every neuron
of its preceding layer).

It should be noted that in conventional dense layers, the activation function
$\sigma$ is applied component-wise, that is,
$\sigma(x_{1},x_{2},\dots,x_{m})=(\sigma^{\prime}(x_{1}),\sigma^{\prime}(x_{2}),\dots,\sigma^{\prime}(x_{m}))$
where $\sigma^{\prime}:\mathbb{R}\to\mathbb{R}$ is a real-valued function,
such as ReLU or $\tanh$. However, in SpectralDense layers, the following
GroupSort function is used as the activation function $\sigma$:

######  Definition 3.

(Anil et al. (2019)) The GroupSort function (of group size 2) is a function
$\sigma:\mathbb{R}^{m}\to\mathbb{R}^{m}$ defined as follows:

If $m$ is even, then
$\sigma([x_{1},x_{2},\cdots,x_{m-1},x_{m}]^{T})=[\max(x_{1},x_{2}),\min(x_{1},x_{2}),\cdots,\max(x_{m-1},x_{m}),\min(x_{m-1},x_{m})]^{T}$
(4a) else, if $m$ is odd,
$\sigma([x_{1},x_{2},\cdots,x_{m-2},x_{m-1},x_{m}]^{T})=[\max(x_{1},x_{2}),\min(x_{1},x_{2}),\cdots,\min(x_{m-2},x_{m-1}),x_{m}]^{T}$
(4b)

For example, in the case where the output layer has dimension $m=4$, we have

$\sigma([0,3,4,2]^{T})=[3,0,4,2]^{T}\;\;\;\sigma([5,3,2,4]^{T})=[5,3,4,2]^{T}$
(4c)

SpectralDense layers can now be defined as follows:

######  Definition 4.

(Serrurier et al. (2021)) SpectralDense layers are dense layers such that 1)
the largest singular value of $W$ is 1, and 2) the activation function
$\sigma$ is GroupSort function.

Therefore, SpectralDense layers are similar to dense layers in terms of their
structure, except that the activation function does not act component-wise and
the weight matrices have a spectral norm of 1. The spectral norm of a matrix
is also equal to the largest singular value in its SVD. Since the largest
singular value of the weight matrix $W$ is 1, the spectral norm $\lVert
W\rVert_{2}$ is also 1. The function $\sigma:\mathbb{R}^{m}\to\mathbb{R}^{m}$
is also 1-Lipschitz continuous (with respect to the Euclidean norm), since it
has Jacobian of spectral norm 1 almost everywhere (everywhere except a set of
measure 0); thus it is also 1-Lipschitz continuous (see Theorem 3.1.6 in
Federer (2014)). We therefore conclude that every SpectralDense layer is also
1-Lipschitz continuous. Next, the definition of the class of LCNNs is given as
follows.

######  Definition 5.

Let $\mathcal{LN}^{m}_{n}$ be the class of Lipschitz-constrained neural
networks (LCNNs) as follows:

$\begin{array}[]{l}\mathcal{LN}_{n}^{m}:=\\{\>f\leavevmode\nobreak\
|\leavevmode\nobreak\ f:\mathbb{R}^{n}\to\mathbb{R}^{m}\>,\>\exists
j\in\mathbb{N}\text{ such that }f=W_{j+1}f_{j}\circ f_{j-1}\circ...\circ
f_{2}\circ f_{1},\\\ \text{ where }f_{i}=\sigma(W_{i}x+b),\text{ and }\,\lVert
W_{i}\rVert_{2}=1,i=1,...,j\,\\}\end{array}$ (5)

where $\sigma$ is the GroupSort activation with group size 2. Thus, each LCNN
in $\mathcal{LN}_{n}^{m}$ consists of many SpectralDense layers (i.e.,
$W_{i}$, $i=1,...,j$) composed together, with a final weight matrix $W_{j+1}$
at the end. The spectral norm constraint is imposed for all weight matrices
except the final weight matrix $W_{j+1}$.

Since the Lipschitz constant for each of the functions $f_{i}$, $i=1,...,j,$
is bounded by 1, it is readily shown that for each neural network in
$\mathcal{LN}_{n}^{m}$, the Lipschitz constant is bounded by the spectral norm
of the final weight matrix $W_{j+1}$. This implies that we can control the
Lipschitz constant of an LCNN by manipulating the spectral norm of the final
weight matrix $W_{j+1}$, and this can be done in a variety of ways, such as
imposing the constraint for the final weight matrix as a whole, or restricting
the absolute value of each entry in the final weight matrix during the
training process. In the simulation study in Section 5, we impose the absolute
value constraint on each entry of the final weight matrix to control the
Lipschitz constant of the LCNN.

######  Remark 1.

The SpectralDense layers adopted in this work differ from those used in Anil
et al. (2019), since the matrix norm used in this work is the spectral norm,
whereas the norm used by Anil et al. is the $\infty$-norm. Although all matrix
norms give rise to equivalent topologies (the same open sets), we use the
spectral norm as it is directly related to the Jacobian of the function.
Specifically, a well-known theorem by Rademacher (see Lemma 3.1.7 of Federer
(2014) for a proof) states that if $X\subset\mathbb{R}^{n}$ is open and
$f:X\to\mathbb{R}^{m}$ is $L$-Lipschitz continuous, then $f$ is almost
everywhere differentiable, and we have

$L=\sup_{x\in X}\,\lVert J_{f}(x)\rVert_{2}$ (6)

Therefore, it is observed that the Lipschitz constant using the Euclidean norm
on the input space $\mathbb{R}^{n}$ and output space $\mathbb{R}^{m}$ is
actually the supremum of the spectral norm of the Jacobian matrix of $f$. If
the $\infty$-norm were to be used, to the best of our knowledge, no such
essential relationship that involves the Lipschitz constant has been proven to
this date.

###  3.2 Robustness of LCNNs

We now discuss how LCNNs using SpectralDense layers can resolve the issue of
sensitivity to input perturbations. Let $f:X\to\mathbb{R}^{m}$ be a neural
network that has been trained using a training algorithm and
$X\subset\mathbb{R}^{n}$ be an open subset such that $f$ is almost everywhere
differentiable with Jacobian $J_{f}$. Given a set of training data, at each
point $x\in X$, one plausible way to maximize the impact of the input
perturbation is to traverse along the direction corresponding to the largest
eigenvalue (the spectral norm) of $J_{f}(x)$ in its SVD, since this is the
direction that leads to the largest variation in output (Szegedy et al.
(2013); Goodfellow et al. (2014)).

For any $f\in\mathcal{LN}^{m}_{n}$, Eq. 5 shows that the Lipschitz constant of
$f$ is bounded by the spectral norm of the final weight matrix. If the
spectral norm of the final matrix is small, the corresponding LCNN will have a
small Lipschitz constant, making it difficult to perturb, even if we travel
along the direction corresponding to the largest eigenvalue of $J_{f}$.
Specifically, if the input perturbation is of size $\delta$, the output change
is at most $L\times\delta$, where $L$ is constrained to be a small value.
Therefore, one plausible method to reduce the sensitivity of neural networks
to input perturbations is to constrain the Lipschitz constant of the networks.

However, since the Lipschitz constant affects the network capacity,
controlling the upper bound of the Lipschitz constant in LCNNs could result in
a reduced network capacity. To address this issue, we will demonstrate in the
next section that the function class $\mathcal{LN}^{m}_{n}$ is a universal
approximator for any Lipschitz continuous target function. Additionally, a
pertinent question that arises is whether, in practice, the Lipschitz
constants of conventional FNNs (e.g., FNNs using conventional dense layers and
ReLU activation functions) are indeed much larger than those in LCNNs. In the
special case of FNNs with ReLU activation functions, Bhowmick et al. (2021)
have designed a provably correct approximation algorithm known as Lipschitz
Branch and Bound (LipBaB), which obtains the Lipschitz constant of such
networks on a compact rectangular domain. In Section 6.7, we will demonstrate
empirically that with noisy training data, FNNs with dense layers could have a
Lipschitz constant several orders of magnitude higher than that of LCNNs.

##  4 Universal Approximation Theorem for LCNNs

This section develops the universal approximation theorem for LCNNs using
SpectralDense layers, which demonstrates that the LCNNs with a bounded
Lipschitz constant can approximate any nonlinear function as long as the
target function is Lipschitz continuous. Before we present the proof for
vector-valued LCNNs that are developed for nonlinear systems with vector-
valued outputs such that of Eq. 1, we first develop a theorem that considers
the approximation of real-valued functions. Then, the results for real-valued
functions can be generalized to the multi-dimensional output case. We first
define the real-valued function class as follows.

$\mathcal{LN}_{n}:=\\{\>f\leavevmode\nobreak\ |\leavevmode\nobreak\
f:\mathbb{R}^{n}\to\mathbb{R}\>,\>\exists
j\in\mathbb{N},s.t.\leavevmode\nobreak\ f=f_{j}\circ...f_{2}\circ
f_{1},f_{i}=\sigma(W_{i}x+b),\lVert W_{i}\rVert_{2}=1,i=1,...,j\\}$ (7)

where $\sigma$ is the GroupSort activation with group size 2. The definition
of Eq. 7 is similar to that of Eq. 5, except that the functions are real-
valued, and the spectral norm for each weight matrix is 1, including the final
weight matrix. Note that the final map $f_{j}$ is an affine map, without any
sorting, since the output of $f_{j}$ is a single real number. It readily
follows that any function in $\mathcal{LN}_{n}$ is also 1-Lipschitz continuous
since the spectral norm of the final weight matrix is one.

Given a target function $F:D\to\mathbb{R}$ where $D$ is a compact and
connected domain and $F$ is Lipschitz continuous, we prove that real-valued
functions from $\mathcal{LN}_{n}$ are universal approximators of 1-Lipschitz
functions, provided that we allow for an amplification of at most $\sqrt{2}$
at the end. In principle, this implies that LCNNs can approximate any
Lipschitz function (i.e., they are dense with respect to the uniform norm on a
compact set). We will follow the notation in Anil et al. (2019), but with
slight modifications. We first present the following definitions, which will
be used in the proof of the universal approximation theorem.

######  Definition 6.

Let $(X,d_{X})$ be a metric space. We use $C_{L}(X,\mathbb{R})$ to denote the
set of all 1-Lipschitz real-valued functions on $X$.

######  Definition 7.

Let $A$ be a set of functions from $\mathbb{R}^{n}$ to $\mathbb{R}$, and let
$k$ be a real number. We define $kA$ as follows.

$kA:=\\{\;cf\;|\;|c|\leq k,\;f\in A\\}$ (8)

######  Definition 8.

A lattice $\mathcal{L}$ in $C_{L}(X,\mathbb{R})$ is a set of functions that is
closed under point-wise maximums and minimums, that is, $\forall
f,g\in\mathcal{L}$, $\min(f,g),\max(f,g)\in\mathcal{L}$.

The following restricted Stone-Weierstrass theorem allows us to approximate
1-Lipschitz continuous functions using lattices.

######  Theorem 1 (Restricted Stone-Weierstrass Anil et al. (2019)).

Let $(X,d_{X})$ be a compact metric space and $\mathcal{L}$ be a lattice in
$C_{L}(X,\mathbb{R})$. Suppose that for all $a,b\in\mathbb{R}$ and $x,y\in X$
such that $|a-b|\leq d_{X}(x,y)$, there exists an $f\in\mathcal{L}$ such that
$f(x)=a$ and $f(y)=b$. Then $\mathcal{L}$ is dense in $C_{L}(X,\mathbb{R})$
with respect to the uniform topology, that is, for every $\epsilon>0$ and for
every $f\in C_{L}(X,\mathbb{R})$, there exists an $\tilde{f}\in\mathcal{L}$
such that

$\sup_{x\in X}\,|f(x)-\tilde{f}(x)|<\epsilon$ (9)

The proof of Theorem 1 is given in Anil et al. (2019) and is omitted here.
Based on Theorem 1, we develop the following theorem to prove that the LCNN
networks constructed using SpectralDense layers can also serve as universal
approximators, if we allow for an amplification of the output at the end. The
proof uses Theorem 1 to show that $\sqrt{2}\mathcal{LN}_{n}\cap
C_{L}(D,\mathbb{R})$ is a lattice. The proof techniques and structure are
similar to those in Anil et al. (2019), while the key difference is the use of
the SVDs of the weight matrices since we are using the spectral norm instead.

######  Theorem 2.

Let $D\subset\mathbb{R}^{n}$ be a compact subset, and $\mathcal{LN}_{n}$ be
the set of LCNNs defined in Eq. 5. $\sqrt{2}\mathcal{LN}_{n}\cap
C_{L}(D,\mathbb{R})$ is dense in $C_{L}(D,\mathbb{R})$ with respect to the
uniform topology, i.e., for every $\epsilon>0$ and for every $f\in
C_{L}(D,\mathbb{R})$, there exists an
$\tilde{f}\in\sqrt{2}\mathcal{LN}_{n}\cap C_{L}(D,\mathbb{R})$ such that

$\sup_{x\in D}\,\lvert f(x)-\tilde{f}(x)\rvert<\epsilon$ (10)

###### Proof.

Theorem 2 states that if we allow for an amplification of $\sqrt{2}$ at the
end, then the LCNNs using SpectralDense layers defined in Eq. 5 can
approximate 1-Lipschitz continuous functions arbitrarily accurately. To prove
Eq. 10, we first show that $\mathcal{L}:=\sqrt{2}\mathcal{LN}_{n}\cap
C_{L}(D,\mathbb{R})$ satisfies the assumptions needed for the restricted
Stone-Weierstrass theorem.

Note that for all $a,b\in\mathbb{R}$ and $x,y\in D$ such that $|a-b|\leq\lVert
x-y\rVert$, there exists an $f\in\mathcal{L}$ such that $f(x)=a$ and $f(y)=b$.
To construct such an $f$, let $f(v)=w^{T}(v-x)+a$ and choose
$w\in\mathbb{R}^{n}$ with $\lVert w\rVert=1$ carefully so that this holds. We
need to ensure that $w^{T}(y-x)=b-a$ such that $f(y)=b$. Therefore, we need to
choose $w$ such that the following equation holds.

$b-a=w^{T}(y-x)\leq\lVert w\rVert\cdot\lVert y-x\rVert=\lVert y-x\rVert$ (11)

We can first set $w$ to be in the same direction as $y-x$, and then gradually
rotate $w$ away from $y-x$ (for example, using a suitable orthogonal matrix)
so that we eventually have $b-a=w^{T}(y-x)$.

Next, we need to show that $\mathcal{L}$ is a lattice. This is equivalent to
showing that $\mathcal{L}$ is closed under pointwise maximums and minimums by
Definition 8. Following the same proof as in Anil et al. (2019), we assume
that $f,g\in\mathcal{L}$ are defined by their weights and biases:

$[W_{1}^{f},b_{1}^{f},W_{2}^{f},b_{2}^{f},\dots,W_{d_{f}}^{f},b_{d_{f}}^{f}]\;\;\;[W_{1}^{g},b_{1}^{g},W_{2}^{g},b_{2}^{g},\dots,W_{d_{g}}^{g},b_{d_{g}}^{g}]$
(12)

where $d_{f}$ and $d_{g}$ represent the depths of the networks $f$ and $g$,
respectively. We assume without loss of generality that the two neural
networks $f$ and $g$ have the same depth, i.e., $d_{f}=d_{g}$. This is
possible since if $d_{f}>d_{g}$, one can pad the neural network with identity
matrix weights and zero biases until they are of the same depth, and likewise
for $d_{f}<d_{g}$. We also assume without loss of generality that each of the
weight matrices (except the final weight matrix) has an even number of rows.
In the case where the neural network has a weight matrix with an odd number of
rows, the weight matrix can be padded with a zero row under the last row and
with a bias $-M$ where $M>0$ is sufficiently large in that row (this is
possible since $D$ is compact and by the extreme value theorem). This is to
prevent a different sorting configuration of the output vector of that layer.
Then, in the next matrix, we add a column of zeros to remove the $-M$ entry.

We now construct a neural network $h$ for $\max(f,g)$ and $\min(f,g)$ with new
weights such that each of the weights satisfies $\lVert W_{i}^{h}\rVert_{2}=1$
for $i=1,\cdots,d_{f}$, and the scaling factor of $\sqrt{2}$ will be applied
later. We first construct a suitable neural network with the layers of $f$ and
$g$ side by side, but with some modifications. Specifically, the first matrix
$W_{1}^{h}$ and the first bias $b_{1}^{h}$ are designed as follows:

$W_{1}^{h}=c[W_{1}^{f},W_{1}^{g}]^{T}\,\,\,\,b_{1}^{h}=[b_{1}^{f},b_{1}^{g}]$
(13)

where $c$ is a positive constant chosen based on $W_{1}^{f}$ and $W_{1}^{g}$
to ensure that $\lVert W_{1}^{h}\rVert_{2}=1$. From Eq. 13, it is shown that
the output of the first layer is obtained by concatenating the output of the
first layer of $f$ and $g$ together, and multiplied by a positive constant
$c$. Note that $1\leq c\leq\frac{1}{\sqrt{2}}$ since each of $W_{1}^{f}$ and
$W_{1}^{g}$ have spectral norms 1. Then, for the rest of the layers ($i\geq
2$), we define

$W_{i}^{h}=\begin{bmatrix}W_{i}^{f}&0\\\ 0&W_{i}^{g}\end{bmatrix}$ (14)

To prove that $\lVert W_{i}^{h}\rVert_{2}=1$, we use the singular value
decomposition.

$\displaystyle W_{i}^{h}$ $\displaystyle=\begin{bmatrix}W_{i}^{f}&0\\\
0&W_{i}^{g}\end{bmatrix}=\begin{bmatrix}U_{i}^{f}D_{i}^{f}{V_{i}^{f}}^{*}&0\\\
0&U_{i}^{g}D_{i}^{g}{V_{i}^{g}}^{*}\end{bmatrix}$ (15a)
$\displaystyle=\begin{bmatrix}U_{i}^{f}&0\\\
0&U_{i}^{g}\end{bmatrix}\begin{bmatrix}D_{i}^{f}&0\\\
0&D_{i}^{g}\end{bmatrix}\begin{bmatrix}{V_{i}^{f}}^{*}&0\\\
0&{V_{i}^{g}}^{*}\end{bmatrix}$ (15b) $\displaystyle=UDV^{T}$ (15c)

It is noted that $\lVert W_{i}^{h}\rVert_{2}$ is simply the largest singular
value of its singular value decomposition. On the right-hand side (RHS) of Eq.
15c, even if the matrix $D$ is not rectangular block diagonal, we can permute
the columns of $D$ and the rows of $V^{T}$ simultaneously, and then the rows
of $D$ and columns of $U$ simultaneously, to obtain a new rectangular block
diagonal matrix. Permuting the columns of $U$ and the rows of $V^{T}$ does not
change the unitary property of these matrices. Therefore, the largest singular
value of $W_{i}^{h}$ is 1 since both $D_{i}^{f}$ and $D_{i}^{g}$ have the
largest singular values of 1, and therefore $\lVert W_{i}^{h}\rVert_{2}=1$. We
also take the following.

$b_{i}^{h}=c[b_{i}^{f},b_{i}^{g}]$ (16)

for each of the biases. Finally, after passing through the GroupSort
activation function, the output of the last layer is

$c[\max(f(x),g(x)),\min(f(x),g(x))]^{T}$ (17)

By passing Eq. 17 through the weight matrix $[1,0]$ or $[0,1]$, we obtain
$h(x)=c\max(f(x),g(x))$ or $h(x)=c\min(f(x),g(x))$ respectively. Since we
consider the set $\mathcal{L}=\sqrt{2}\mathcal{LN}_{n}\cap
C_{L}(D,\mathbb{R})$, and $1\leq c^{-1}\leq\sqrt{2}$, we observe that
$c^{-1}\,h(x)=\max(f(x),g(x))$ or $c^{-1}\,h(x)=\min(f(x),g(x))$. This proves
that $\max(f(x),g(x))$ and $\min(f(x),g(x))$ are in
$\sqrt{2}\mathcal{LN}_{n}\cap C_{L}(D,\mathbb{R})$, which implies that
$\mathcal{L}$ is a lattice. Therefore, $\sqrt{2}\mathcal{LN}_{n}\cap
C_{L}(D,\mathbb{R})$ satisfies the assumptions needed for the restricted
Stone-Weierstrass theorem in Theorem 1, and this completes the proof. ∎

Theorem 2 implies that for 1-Lipschitz target functions, if we allow for an
amplification of $\sqrt{2}$ at the last layer, then LCNNs using SpectralDense
layers are universal approximators for the target function. Similarly, LCNNs
can approximate an $L$-Lipschitz continuous function so long as an
amplification of $\sqrt{2}L$ is allowed. The final amplification can be easily
implemented by using a suitable weight matrix such as a constant multiple of
the identity matrix as the final layer.

######  Remark 2.

The above theorem can be generalized to regression problems from
$\mathbb{R}^{n}$ to $\mathbb{R}^{m}$. Specifically, a similar result can be
derived for LCNNs developed for $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ by
approximating each component function $f_{i}$ for $i=1,...,m$ with a
1-Lipschitz neural network, and then padding the neural networks together to
form a large neural network with $m$ outputs. However, note that, as a result,
the resulting approximating function might be $\sqrt{m}$-Lipschitz continuous.

##  5 Preventing Over-fitting and Improving Generalization with LCNNs

In this section, we provide a theoretical argument using the empirical
Rademacher complexity to show that LCNNs can prevent over-fitting and
therefore generalize better than conventional FNNs by computing the empirical
Rademacher complexity of the LCNNs, and comparing that with usual FNNs made
using the same architecture (i.e., same number of neurons in each layer and
same depth). Specifically, we develop a bound of the empirical Rademacher
complexity (ERC) for any FNN that utilizes the GroupSort function (hereby
called GroupSort Neural Networks), and subsequently use this to obtain a bound
for LCNNs using SpectralDense layers as an immediate corollary. Then, we
compare this bound with the bound for the empirical Rademacher complexity of
FNNs using 1-Lipschitz component-wise activation functions (e.g. ReLU) and
show that the GroupSort NNs achieve a tighter generalization error bound.

###  5.1 Assumptions and Preliminaries

We first assume that the input domain is a bounded subset
$\mathcal{X}\subset\mathbb{R}^{d_{x}}$, where for each $x\in\mathcal{X}$, one
has $\lVert x\rVert\leq B$, and the output space is a subset of the vector
space $\mathcal{Y}\subset\mathbb{R}^{d_{y}}$. We also assume that each input-
output pair in the dataset $(x,y)$ is drawn from some distribution
$\mathcal{D}\subset\mathcal{X}\,\times\,\mathcal{Y}$ with some probability
distribution $\mathbb{P}$. We assume that
$L:\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}^{\geq 0}$ is a loss function that
satisfies the following two properties: 1) there exists a positive $M>0$ such
that for all $y,y^{\prime}\in\mathbb{R}$ $L(y,y^{\prime})\leq M$, and 2) for
all $y^{\prime}\in\mathbb{R}^{d_{y}}$, the function $y\to L(y,y^{\prime})$ is
$L_{r}$-Lipschitz for some $L_{r}>0$.

Let $h:\mathcal{X}\to\mathcal{Y}$ be a function that represents a neural
network model (termed hypothesis) in a hypothesis class.

######  Definition 9.

We define generalization error as

$\underset{(x,y)\sim\mathcal{D}}{\mathbb{E}}\,[L(h(x),y)]=\int_{\mathcal{X}\times\mathcal{Y}}L(h(x),y)\,\mathbb{P}(dx\times
dy)$ (18)

######  Definition 10.

For any dataset
$S=\\{s_{1},s_{2},\cdots,s_{m}\\}\subset\mathcal{X}\times\mathcal{Y}$ with
$s_{i}=(x_{i},y_{i})$, we define the empirical error as

$\hat{\mathbb{E}}_{S}[L(h(x),y)]=\frac{1}{m}\sum_{i=1}^{m}L(h(x_{i}),y_{i})$
(19)

The generalization error is a measure of how well a hypothesis
$h:\mathcal{X}\to\mathcal{Y}$ generalizes from the training dataset to the
entire domain being considered. On the other hand, the empirical error is a
measure of how well the hypothesis $h$ performs on the available data points
$S$.

######  Remark 3.

In this work, we use the $\ell_{2}$ error for the loss function $L$ in the NN
training process, i.e., $L(y,y^{\prime})=\lVert y-y^{\prime}\rVert^{2}$, which
is locally Lipschitz continuous but not globally Lipschitz continuous.
However, since we consider an input domain of $D\times U$ which is a compact
set, and the function $\tilde{F}$ is also continuous, the range of $\tilde{F}$
is also compact and bounded. As a result, the output domain
$\mathcal{Y}\subset\mathbb{R}^{d_{y}}$ is bounded, and the assumptions for the
loss function $L(\cdot,\cdot)$ are satisfied using the $\ell_{2}$ error.

###  5.2 Empirical Rademacher Complexity Bound for GroupSort Neural Networks

We first define the empirical Rademacher complexity of a real-valued function
hypothesis class.

######  Definition 11.

Given an input domain space $\mathcal{X}\subset\mathbb{R}^{d_{x}}$, suppose
$H$ is a class of real-valued functions from $\mathcal{X}$ to $\mathbb{R}$.
Let $S=\\{x_{1},x_{2},\cdots,x_{m}\\}\subset\mathcal{X}$, which is a set of
samples from $\mathcal{X}$. The empircial Rademacher complexity (ERC) of $S$
with respect to $H$, denoted by $\mathcal{R}_{S}(H)$, is defined as:

$\mathcal{R}_{S}(H):=\underset{\epsilon}{\mathbb{E}}\,\sup_{h\in
H}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})$ (20)

where each of $\epsilon_{i}$ are i.i.d Rademacher variables, i.e., with
$\mathbb{P}(\epsilon_{i}=1)=\frac{1}{2}$ and
$\mathbb{P}(\epsilon_{i}=-1)=\frac{1}{2}$.

The ERC measures the richness of a real-valued function hypothesis class with
respect to a probability distribution. A more complex hypothesis class with a
larger ERC is likely to represent a richer variety of functions, but may also
lead to over-fitting, especially in the presence of noise or insufficient
data. The ERC is often used to obtain a probabilistic bound on the
generalization error in statistical machine learning. Specifically, we first
present the following theorem in Mohri et al. (2018) that obtains a
probabilistic upper bound on the generalization error for a hypothesis class.

######  Theorem 3 (Theorem 3.3 in Mohri et al. (2018)).

Let $\mathcal{H}$ be a hypothesis class of functions
$h:\mathcal{X}\subset\mathbb{R}^{d_{x}}\to\mathcal{Y}\subset\mathbb{R}^{d_{y}}$
and $L:\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}^{\geq 0}$ be the loss
function that satisfies the properties in Section 5.1. Let $\mathcal{G}$ be a
hypothesis class of loss functions associated with the hypotheses
$h\in\mathcal{H}$:

$\mathcal{G}:=\\{g:(x,y)\to L(h(x),y)\,,\,h\in\mathcal{H}\\}$ (21)

For any set of $m$ i.i.d. training samples $S=\\{s_{1},s_{2},\cdots,s_{m}\\}$,
$s_{i}=(x_{i},y_{i})$, $i=1,...,m$, drawn from a probability distribution
$\mathcal{D}\subset\mathcal{X}\times\mathcal{Y}$, for any $\delta\in(0,1)$,
the following upper bound holds with probability $1-\delta$ :

$\underset{(x,y)\sim\mathcal{D}}{E}\,[L(h(x),y)]\leq\frac{1}{m}\sum_{i=1}^{m}L(h(x_{i}),y_{i})+2\mathcal{R}_{S}(\mathcal{G})+3M\sqrt{\frac{\log\frac{1}{\delta}}{2m}}$
(22)

Eq. 21 in Theorem 3 explains why an overtly large hypothesis class with a high
degree of complexity could inevitably increase the generalization error.
Specifically, if the hypothesis class is overtly large, the first term in the
RHS of Eq. 21 (i.e., the empirical error
$\hat{E}_{S}[L(h(x),y)]=\sum_{i}L(h(x_{i}),y_{i})$) is expected to be
sufficiently small, since one could find a hypothesis $h$ that minimizes the
training error using an appropriate training algorithm. However, the trade-off
is that the ERC $\mathcal{R}_{S}(\mathcal{G})$ will increase with the size of
the hypothesis class. On the contrary, if the hypothesis class is limited, the
model complexity represented by the ERC of $\mathcal{R}_{S}(\mathcal{G})$ is
reduced, while the empirical error $\hat{E}_{S}[L(h(x),y)]$ is expected to
increase. Since the ERC of the hypothesis class of loss functions
$\mathcal{R}_{S}(\mathcal{G})$ appears on the RHS of the inequality above, it
implies that a lower ERC leads to a tighter generalization error bound.
Therefore, in this section, we demonstrate that LCNNs using SpectralDense
layers have a smaller ERC $\mathcal{R}_{S}(\mathcal{G})$ than conventional
FNNs using dense layers.

Since $\mathcal{Y}\subset\mathbb{R}^{d_{y}}$ is a subset of a real vector
space, while the definition of ERC in Definition 11 is with respect to a real-
valued class of functions, to simplify the discussion, we can first apply the
following contraction inequality.

######  Lemma 4 (Vector Contraction Inequality Maurer (2016)).

Let $\mathcal{H}$ be a hypothesis class of functions from
$\mathcal{X}\subset\mathbb{R}^{d_{x}}$ to
$\mathcal{Y}\subset\mathbb{R}^{d_{y}}$. For any set of data points
$S=\\{s_{1},s_{2},\cdots,s_{m}\\}\subset\mathcal{X}\times\mathcal{Y}$ and
$L_{r}$-Lipschitz loss function $y\to L(y,y^{\prime})$ for some $L_{r}>0$, we
have

$\mathcal{R}_{S}(\mathcal{G})=\underset{\epsilon}{\mathbb{E}}\,\sup_{h\in\mathcal{H}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}L(h(x_{i}),y)\leq\sqrt{2}L_{r}\underset{\epsilon}{\mathbb{E}}\,\sup_{h\in\mathcal{H}}\frac{1}{m}\sum_{i=1}^{m}\sum_{k=1}^{d_{y}}\epsilon_{ik}h_{k}(x_{i})$
(23)

where $h_{k}$ is the $k^{\text{th}}$ component function of $h\in\mathcal{H}$
and each of $\epsilon_{ik}$ are i.i.d. Rademacher variables.

The proof of the above inequality can be found in Maurer (2016), and is
omitted here. The RHS of the inequality in Eq. 23 can be bounded by taking out
the sum in the supremum:

$\underset{\epsilon}{\mathbb{E}}\sup_{h\in\mathcal{H}}\frac{1}{m}\sum_{i=1}^{m}\sum_{k=1}^{d_{y}}\epsilon_{ik}h_{k}(x_{i})\leq\sum_{k=1}^{d_{y}}\underset{\epsilon}{\mathbb{E}}\sup_{h_{k}\in\mathcal{H}_{k}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{ik}h_{k}(x_{i})$
(24)

where each of the $\mathcal{H}_{k}$ is a real-valued function class that
correspond to the $k^{\text{th}}$ component function of each function in
$h\in\mathcal{H}$. The RHS of Eq. 24 is the sum of the Rademacher complexities
of the hypothesis classes $\mathcal{H}_{k}$. Therefore, to obtain an upper
bound for $\mathcal{R}_{S}(\mathcal{G})$, we can first consider the case of a
real-valued function hypothesis class, and then extend the results to the
multidimensional case by applying Eq. 24. Specifically, we first develop an
ERC bound for GroupSort Neural Networks, which refers to any FNN that utilizes
the GroupSort activation function. Since LCNNs using SpectralDense layers are
FNNs that also utilize the GroupSort activation function, the ERC bound for
LCNNs using SpectralDense layers will follow as an immediate corollary.

The following definitions are first presented to define the classes of real-
valued GroupSort neural networks and of real-valued LCNNs using SpectralDense
layers.

######  Definition 12.

We use $\mathcal{H}_{d}$ to denote the hypothesis class of GroupSort neural
networks with depth $d$ that map the input domain
$\mathcal{X}\subset\mathbb{R}^{d_{x}}$ to $\mathbb{R}$:

$x\to W_{d}\,\sigma(W_{d-1}\,\sigma(\cdots\sigma(W_{1}x)))$ (25)

where each of the weight matrices $W_{i}\in\mathbb{R}^{m_{i}\times n_{i}}$ has
a bounded Frobenius norm, that is, $\lVert W_{i}\rVert_{F}\leq R_{i}$ for some
$R_{i}\geq 0$, $m_{i}$ and $n_{i}$ are even for all $i=1,...,d-1$ except for
$n_{1}$, and $\sigma$ is the GroupSort function with group size 2.

######  Definition 13.

We use $\mathcal{H}^{SD}_{d}$ to denote the hypothesis class of LCNNs using
SpectralDense layers with depth $d$ that map the input domain
$\mathcal{X}\subset\mathbb{R}^{d_{x}}$ to $\mathbb{R}$:

$x\to W_{d}\,\sigma(W_{d-1}\,\sigma(\cdots\sigma(W_{1}x)))$ (26)

where each of the weight matrices satisfies $\lVert W_{i}\rVert_{2}=1$,
$m_{i}$ and $n_{i}$ are even for all $i=1,...,d-1$ except for $n_{1}$, and
$\sigma$ is the GroupSort function with group size 2.

Note that the key difference between $\mathcal{H}_{d}$ and
$\mathcal{H}^{SD}_{d}$ is that $\mathcal{H}_{d}$ is defined for any FNN using
GroupSort activation functions, while $\mathcal{H}^{SD}_{d}$ is defined for
LCNNs that use the GroupSort activation function and satisfy $\lVert
W_{i}\rVert_{2}=1$. Additionally, the Frobenius norm is used in
$\mathcal{H}_{d}$, while the spectral norm is used in $\mathcal{H}^{SD}_{d}$.
Despite the differences between $\mathcal{H}_{d}$ and $\mathcal{H}^{SD}_{d}$,
it will be demonstrated from the following lemmas and theorems that the two
classes $\mathcal{H}_{d}$ and $\mathcal{H}^{SD}_{d}$ are highly related.

######  Remark 4.

The bias term is ignored in Eq. 25 and Eq. 26 to simplify the formulation, as
in principle, we can take into account the bias term by padding the input with
a vector consisting of ones, and introducing another weight matrix into each
of the $W_{i}$. Additionally, we assume that $m_{i}$ and $n_{i}$ are even
without any loss of generality, since this does not affect the expressiveness
of the class of networks using LCNNs using SpectralDense layers, as shown in
the proof of Theorem 2.

Next, we develop a bound on the ERC of $\mathcal{H}_{d}$, i.e.,
$\mathcal{R}_{S}(\mathcal{H}_{d})$, and then the ERC bound for
$\mathcal{H}^{SD}_{d}$ follows as an immediate corollary. The main intuition
to obtain such an upper bound for $\mathcal{R}_{S}(\mathcal{H}_{d})$ is to
recursively “peel off” the weight matrices and activation functions. Such
methods were used in Golowich et al. (2018); Neyshabur et al. (2015); Wu et
al. (2021b, 2022), where the activation function was applied element-wise. The
key difficulty in our setting stems from the fact that the GroupSort
activation function $\sigma$ is not an element-wise function. To address this
issue, we first represent the functions $\max(a,b)$ and $\min(a,b)$ as
follows:

$\max(a,b)=\frac{1}{2}(a+b+|b-a|)\,\,\,\,\min(a,b)=\frac{1}{2}(a+b-|b-a|)$
(27)

Before we present the results for peeling off the weight matrices of FNNs, the
following definition is first given and will be used in the proof of Lemma 5
that peels off one GroupSort activation function layer.

######  Definition 14.

For $d\geq 1$, we define $\tilde{\mathcal{H}}_{d}$ as the class of (vector-
valued) functions on the input domain $\mathcal{X}\subset\mathbb{R}^{d_{x}}$
of the form

$x\to\sigma(W_{d}\,\sigma(\cdots\sigma(W_{1}x)))$ (28)

where each of the weight matrices $W_{i}\in\mathbb{R}^{m_{i}\times n_{i}}$ has
a bounded Frobenius norm, that is, $\lVert W_{i}\rVert_{F}\leq R_{i}$ for some
$R_{i}\geq 0$, $m_{i}$ and $n_{i}$ are even for all $i=1,...,d$ except
$n_{1}$, and $\sigma$ is the GroupSort function with group size 2. If $d=0$,
we define $\tilde{H}_{0}$ as a hypothesis class that contains only the
identity map on $\mathcal{X}$.

Note that Definition 14 is very similar to Definition 12, but the last layer
has been removed, so that the resultant hypothesis class is vector-valued.
Subsequently, the following lemma provides a way to “peel off” layers using
the GroupSort function, which is the main tool used in the derivation of ERC
for GroupSort NNs.

######  Lemma 5.

Let $\tilde{H}_{d}$ be a vector-valued hypothesis class defined in Definition
14, with $d\geq 1$, and suppose that $\lVert W_{d}\rVert_{F}\leq R_{d}$. For
any dataset with $m$ data points, we have the following inequality:

$\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert\leq
2R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\left\lVert\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert$
(29)

###### Proof.

Letting $w_{1}^{T},w_{2}^{T},\cdots,w_{k}^{T}$ represent the rows of $W_{d}$,
we have

$\displaystyle\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert$
$\displaystyle=\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}\sigma(W_{d}h(x_{i}))\right\rVert$
(30a) $\displaystyle=\underset{\epsilon}{\mathbb{E}}\sup_{\lVert
W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}\sigma\begin{bmatrix}w_{1}^{T}h(x_{i})\\\
w_{2}^{T}h(x_{i})\\\ \vdots\\\ w_{k}^{T}h(x_{i})\end{bmatrix}\right\rVert$
(30b)

By expanding the components of $\sigma$ using the identities for the maximum
and minimum in Eqs. 27, Eq. 30 can be written as

$\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{2m}\sum_{i=1}^{m}\epsilon_{i}\begin{bmatrix}w_{1}^{T}h(x_{i})+w_{2}^{T}h(x_{i})+|w_{1}^{T}h(x_{i})-w_{2}^{T}h(x_{i})|\\\
w_{1}^{T}h(x_{i})+w_{2}^{T}h(x_{i})-|w_{1}^{T}h(x_{i})-w_{2}^{T}h(x_{i})|\\\
\vdots\\\
w_{k-1}^{T}h(x_{i})+w_{k}^{T}h(x_{i})+|w_{k-1}^{T}h(x_{i})-w_{k}^{T}h(x_{i})|\\\
w_{k-1}^{T}h(x_{i})+w_{k}^{T}h(x_{i})-|w_{k-1}^{T}h(x_{i})-w_{k}^{T}h(x_{i})|\end{bmatrix}\right\rVert$
(31)

We can bound Eq. 31 using the triangle inequality with $A_{1}$ and $A_{2}$
defined as follows.

$\displaystyle A_{1}=$
$\displaystyle\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{2m}\sum_{i=1}^{m}\epsilon_{i}\begin{bmatrix}w_{1}^{T}h(x_{i})\\\
w_{2}^{T}h(x_{i})\\\ \vdots\\\ w_{k-1}^{T}h(x_{i})\\\
w_{k}^{T}h(x_{i})\end{bmatrix}\right\rVert+\underset{\epsilon}{\mathbb{E}}\sup_{\lVert
W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{2m}\sum_{i=1}^{m}\epsilon_{i}\begin{bmatrix}w_{2}^{T}h(x_{i})\\\
w_{1}^{T}h(x_{i})\\\ \vdots\\\ w_{k}^{T}h(x_{i})\\\
w_{k-1}^{T}h(x_{i})\end{bmatrix}\right\rVert$ (32a)

$A_{2}=\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{2m}\sum_{i=1}^{m}\epsilon_{i}\begin{bmatrix}|w_{1}^{T}h(x_{i})-w_{2}^{T}h(x_{i})|\\\
-|w_{1}^{T}h(x_{i})-w_{2}^{T}h(x_{i})|\\\ \vdots\\\
|w_{k-1}^{T}h(x_{i})-w_{k}^{T}h(x_{i})|\\\
-|w_{k-1}^{T}h(x_{i})-w_{k}^{T}h(x_{i})|\end{bmatrix}\right\rVert$ (32b)

We first bound $A_{1}$ appropriately. By noting that the first term in Eq. 32a
is equal to the second term since we only permuted the components, the
following equation is derived for $A_{1}$.

$A_{1}=\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}\begin{bmatrix}w_{1}^{T}h(x_{i})\\\
w_{2}^{T}h(x_{i})\\\ \vdots\\\ w_{k-1}^{T}h(x_{i})\\\
w_{k-1}^{T}h(x_{i})\end{bmatrix}\right\rVert=\underset{\epsilon}{\mathbb{E}}\sup_{\lVert
W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}W_{d}h(x_{i})\right\rVert$
(33)

The RHS of Eq. 33 can be bounded using similar strategies that can be found in
Lemma 3.1 of Golowich et al. (2018). Subsequently, we rewrite the expression
inside the supremum and expectation of the first term as follows.

$\sqrt{\sum_{j=1}^{k}\lVert
w_{j}\rVert^{2}\biggl{(}\sum_{i=1}^{m}\biggl{(}\epsilon_{i}\frac{w_{j}^{T}}{\lVert
w_{j}\rVert}h(x_{i})\biggl{)}\biggl{)}^{2}}$ (34)

Note that Eq. 34 is maximized when one of the $\lVert w_{j}\rVert=R_{d}$ and
the rest of $\lVert w_{l}\rVert=0$ for $l\neq j$ (this is simply because we
maximize a positive linear function over $\lVert w_{1}\rVert^{2},\lVert
w_{2}\rVert^{2},\cdots,\lVert w_{h}\rVert^{2}$). Therefore, it follows that

$\displaystyle\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}\begin{bmatrix}w_{1}^{T}h(x_{i})\\\
w_{2}^{T}h(x_{i})\\\ \vdots\\\ w_{k-1}^{T}h(x_{i})\\\
w_{k-1}^{T}h(x_{i})\end{bmatrix}\right\rVert$
$\displaystyle=\underset{\epsilon}{\mathbb{E}}\sup_{\lVert w\rVert\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}w^{T}h(x_{i})$
$\displaystyle\leq\underset{\epsilon}{\mathbb{E}}\sup_{\lVert w\rVert\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\lVert
w\rVert\bigg{\|}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\bigg{\|}$
$\displaystyle=R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\bigg{\|}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\bigg{\|}$
(35)

The second term $A_{2}$ can be rewritten in the following way.

$A_{2}=\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{2m}\sum_{i=1}^{m}\epsilon_{i}\begin{bmatrix}|w_{1}^{T}h(x_{i})+w_{2}^{T}h(x_{i})|\\\
|w_{1}^{T}h(x_{i})+w_{2}^{T}h(x_{i})|\\\ \vdots\\\
|w_{k-1}^{T}h(x_{i})+w_{k}^{T}h(x_{i})|\\\
|w_{k-1}^{T}h(x_{i})+w_{k}^{T}h(x_{i})|\end{bmatrix}\right\rVert$ (36)

Note that the negative signs in even coordinates can be removed since we take
the magnitude of the vector, and removal of negative signs in even rows
$w_{2},w_{4},\cdots,w_{k}$ is possible by the symmetry of the set $\lVert
W_{d}\rVert_{F}\leq R_{d}$. Then, we rewrite the inner expression of Eq. 36 as
follows.

$\frac{1}{2m}\sqrt{\sum_{j=1}^{k/2}2\lVert
w_{2j}+w_{2j-1}\rVert^{2}\biggl{(}\sum_{i=1}^{m}\epsilon_{i}\bigg{|}\frac{(w_{2j}+w_{2j-1})^{T}}{\lVert
w_{2j}+w_{2j-1}\rVert}h(x_{i})\bigg{|}\biggl{)}^{2}}$ (37)

Using the triangle inequality property of norms, Eq. 37 can be bounded by

$\frac{1}{2m}\sqrt{\sum_{j=1}^{k/2}4\bigg{(}\lVert w_{2j}\rVert^{2}+\lVert
w_{2j-1}\rVert^{2}\biggl{)}\biggl{(}\sum_{i=1}^{m}\epsilon_{i}\bigg{|}\frac{(w_{2j}+w_{2j-1})^{T}}{\lVert
w_{2j}+w_{2j-1}\rVert}h(x_{i})\bigg{|}\biggl{)}^{2}}$ (38)

It is readily shown that Eq. 38 is maximized when $\lVert
w_{2j}\rVert^{2}+\lVert w_{2j-1}\rVert^{2}=R_{d}^{2}$ for some $j$, by a
similar argument to the one for $A_{1}$. Thus, Eq. 38 can be further bounded
by

$\frac{1}{2m}\sqrt{4R_{d}^{2}\biggl{(}\sum_{i=1}^{m}\epsilon_{i}\bigg{|}\frac{(w_{2j}+w_{2j-1})^{T}}{\lVert
w_{2j}+w_{2j-1}\rVert}h(x_{i})\bigg{|}\biggl{)}^{2}}=\frac{R_{d}}{m}\bigg{|}\sum_{i=1}^{m}\epsilon_{i}\bigg{|}\frac{(w_{1}+w_{2})^{T}}{\lVert
w_{1}+w_{2}\rVert}h(x_{i})\bigg{|}\bigg{|}$ (39)

Therefore, we finally derive the bound for $A_{2}$ defined in Eq. 32b as
follows.

$\displaystyle A_{2}$
$\displaystyle=\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{2m}\sum_{i=1}^{m}\epsilon_{i}\begin{bmatrix}|w_{1}^{T}h(x_{i})+w_{2}^{T}h(x_{i})|\\\
|w_{1}^{T}h(x_{i})+w_{2}^{T}h(x_{i})|\\\ \vdots\\\
|w_{k-1}^{T}h(x_{i})+w_{k}^{T}h(x_{i})|\\\
|w_{k-1}^{T}h(x_{i})+w_{k}^{T}h(x_{i})|\end{bmatrix}\right\rVert$
$\displaystyle\leq\,\underset{\epsilon}{\mathbb{E}}\sup_{\lVert
W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\frac{R_{d}}{m}\bigg{|}\sum_{i=1}^{m}\epsilon_{i}\bigg{|}\frac{(w_{1}+w_{2})^{T}}{\lVert
w_{1}+w_{2}\rVert}h(x_{i})\bigg{|}\bigg{|}$ $\displaystyle\leq
R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\bigg{|}\sum_{i=1}^{m}\epsilon_{i}\frac{(w_{1}+w_{2})^{T}}{\lVert
w_{1}+w_{2}\rVert}h(x_{i})\bigg{|}\qquad\text{by Talagrand's Contraction
Lemma}$ $\displaystyle\leq R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{\lVert
W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\left\lVert\frac{(w_{1}+w_{2})^{T}}{\lVert
w_{1}+w_{2}\rVert}\right\rVert\left\lVert\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert$
$\displaystyle\leq
R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\left\lVert\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert$
(40)

In the second inequality, a slightly different version of Talagrand’s
Contraction Lemma is used, which can be found in Mohri et al. (2018).
Combining $A_{1}$ and $A_{2}$, we derive Eq. 29 as follows.

$\displaystyle\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert$
$\displaystyle\leq A_{1}+A_{2}$ $\displaystyle\leq
R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\left\lVert\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert+R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\left\lVert\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert$
$\displaystyle=2R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\left\lVert\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert$
(41)

This completes the proof. ∎

######  Theorem 6.

Assume that $d\geq 1$ and $\mathcal{H}_{d}$ is the hypothesis class of real-
valued functions defined in Definition 12. Suppose that
$\mathcal{X}\subset\mathbb{R}^{d_{x}}$ is a bounded subset such that for all
$x\in\mathcal{X},\lVert x\rVert\leq B$. For any set of $m$ training samples
$S=\\{x_{1},x_{2},\cdots,x_{m}\\}$, we have

$\mathcal{R}_{S}(\mathcal{H}_{d})\leq\frac{B}{\sqrt{m}}2^{d-1}\Pi_{i=1}^{d}R_{i}$
(42)

where $\lVert W_{i}\rVert_{F}\leq R_{i}$ for each weight matrix in
$\mathcal{H}_{d}$.

###### Proof.

Since the proof is very similar to that in Golowich et al. (2018), we provide
a proof sketch only for clarity. Note that

$\mathcal{R}_{S}(\mathcal{H}_{d})=\underset{\epsilon}{\mathbb{E}}\sup_{h\in\mathcal{H}_{d}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})=\underset{\epsilon}{\mathbb{E}}\sup_{\lVert
W_{d}\rVert_{2}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}W_{d}h(x_{i})$
(43)

Since $W_{d}$ has only 1 row, it follows that $\lVert W_{d}\rVert_{2}=\lVert
W_{d}\rVert_{F}$ by Lemma 7. By the Cauchy-Schwarz inequality, the RHS of Eq.
42 can be bounded by:

$\displaystyle\underset{\epsilon}{\mathbb{E}}\sup_{\lVert W_{d}\rVert_{F}\leq
R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}W_{d}h(x_{i})$
$\displaystyle\leq\underset{\epsilon}{\mathbb{E}}\sup_{\lVert
W_{d}\rVert_{2}\leq R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\lVert
W_{d}\rVert_{2}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert$
$\displaystyle\leq
R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert$
(44)

By recursively applying Lemma 5, we derive the following inequality:

$R_{d}\,\underset{\epsilon}{\mathbb{E}}\sup_{h\in\tilde{\mathcal{H}}_{d-1}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\rVert\leq(2^{d-1}\Pi_{i=1}^{d}R_{i})\,\underset{\epsilon}{\mathbb{E}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}x_{i}\right\rVert$
(45)

The expected value can then be bounded using Jensen’s inequality:

$\displaystyle\underset{\epsilon}{\mathbb{E}}\left\lVert\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}x_{i}\right\rVert$
$\displaystyle\leq\frac{1}{m}\sqrt{\underset{\epsilon}{\mathbb{E}}\left\lVert\sum_{i=1}^{m}\epsilon_{i}x_{i}\right\rVert^{2}}$
$\displaystyle=\frac{1}{m}\sqrt{\underset{\epsilon}{\mathbb{E}}\sum_{i=1}^{m}\sum_{j=1}^{m}\epsilon_{i}\epsilon_{j}x_{j}^{T}x_{i}}$
$\displaystyle=\frac{1}{m}\sqrt{\underset{\epsilon}{\mathbb{E}}\,m\lVert
x\rVert^{2}}=\frac{B}{\sqrt{m}}$ (46)

Therefore, the following inequality is derived, and this completes the proof.

$\mathcal{R}_{S}(\mathcal{H}_{d})\leq\frac{B}{\sqrt{m}}2^{d-1}\Pi_{i=1}^{d}R_{i}$
(47)

∎

Theorem 6 develops the upper bound for the ERC of GroupSort NNs. Based on the
results derived in Theorem 6, we subsequently derive the bound for the ERC of
LCNNs. We begin with a lemma that relates the Frobenius norm to the spectral
norm.

######  Lemma 7.

Given a weight matrix $W\in\mathbb{R}^{m\times n}$, the following inequality
holds:

$\lVert W\rVert_{2}\leq\lVert W\rVert_{F}\leq\min(m,n)\lVert W\rVert_{2}$ (48)

The proof can be readily obtained based on the fact that $\lVert W\rVert_{F}$
is the norm of the vector consisting of all singular values of $W$, and the
spectral norm is simply the largest singular value. Using this lemma, we
develop the following bound for $\mathcal{R}_{S}(\mathcal{H}^{SD}_{d})$, where
$\mathcal{H}^{SD}_{d}$ is the hypothesis class of LCNNs using SpectralDense
layers of depth $d$.

######  Corollary 1.

Let $\mathcal{H}^{SD}_{d}$ be the real-valued function hypothesis class
defined in Definition 13 with $d\geq 1$. Suppose that
$\mathcal{X}\subset\mathbb{R}^{d_{x}}$ is a bounded subset such that for all
$x\in\mathcal{X},\lVert x\rVert\leq B$. For any set of $m$ training samples
$S=\\{x_{1},x_{2},\cdots,x_{m}\\}$, we have the following inequality:

$\mathcal{R}_{S}(\mathcal{H}^{SD}_{d})\leq\frac{B}{\sqrt{m}}2^{d-1}\Pi_{i=1}^{d}\min(m_{i},n_{i})=\frac{B}{\sqrt{m}}2^{d-1}\Pi_{i=1}^{d-1}\min(m_{i},n_{i})$
(49)

where $m_{i}$ and $n_{i}$ are the number of rows and columns of the $i^{th}$
weight matrix $W_{i}$.

###### Proof.

The proof can be readily obtained using the fact that $\lVert
W_{i}\rVert_{2}=1$ for all functions in $\mathcal{H}^{SD}_{d}$ and Lemma 7
(i.e., $\lVert W_{i}\rVert_{F}\leq\min(m_{i},n_{i})\lVert
W_{i}\rVert_{2}=\min(m_{i},n_{i})$). Then, by substituting this inequality
into Theorem 6, we derive the inequality in Eq. 49. The last equality
immediately follows since $m_{d}=1$, and this completes the proof. ∎

######  Remark 5.

Note that the bound derived in Eq. 49 is a completely size-dependent bound for
the ERC of $\mathcal{H}^{SD}_{d}$. Therefore, once the architecture of the
LCNNs using SpectralDense layers has been decided, the ERC of the set of LCNNs
is bounded by a constant, which only depends on the neurons in each layer, and
most importantly, not on the choice of weights in each weight matrix.

Subsequently, using the vector contraction inequality in Lemma 4, we derive
the following corollary that generalizes the results to the multi-dimensional
output case.

######  Corollary 2.

Suppose that $\mathcal{X}\subset\mathbb{R}^{d_{x}}$ is a bounded subset such
that for all $x\in\mathcal{X},\lVert x\rVert\leq B$. Let $\mathcal{G}$ be the
hypothesis class of loss functions defined in Eq. 21, where $\mathcal{H}$ is
the hypothesis class of functions from $\mathcal{X}$ to
$\mathcal{Y}\subset\mathbb{R}^{d_{y}}$ with each component function
$h_{k}\in\mathcal{H}^{SD}_{d}$ for $k=1,...,d_{y}$, and
$L:\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}^{\geq 0}$ is a $L_{r}$-Lipschitz
loss function that satisfies the properties in Section 5.1. Then we have

$\mathcal{R}_{S}(\mathcal{G})\leq\frac{\sqrt{2}d_{y}L_{r}B}{\sqrt{m}}2^{d-1}\Pi_{i=1}^{d-1}\min(m_{i},n_{i})$
(50)

where $m_{i}$ and $n_{i}$ are the number of rows and columns in $W_{i}$.

###### Proof.

The proof of Eq. 50 can be readily obtained by applying Corollary 1 to Eq. 23
in Lemma 4 and using Eq. 24 as follows.

$\displaystyle\mathcal{R}_{S}(\mathcal{G})=\underset{\epsilon}{\mathbb{E}}\,\sup_{h\in\mathcal{H}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}L(h(x_{i}),y)$
$\displaystyle\leq\sqrt{2}L_{r}\underset{\epsilon}{\mathbb{E}}\,\sup_{h\in\mathcal{H}}\frac{1}{m}\sum_{i=1}^{m}\sum_{k=1}^{d_{y}}\epsilon_{ik}h_{k}(x_{i})$
$\displaystyle\leq\sqrt{2}L_{r}\sum_{k=1}^{d_{y}}\underset{\epsilon}{\mathbb{E}}\sup_{h_{k}\in\mathcal{H}^{SD}_{d}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{ik}h_{k}(x_{i})$
$\displaystyle=\sqrt{2}L_{r}d_{y}\,\underset{\epsilon}{\mathbb{E}}\sup_{h_{k}\in\mathcal{H}^{SD}_{d}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{ik}h_{k}(x_{i})$
$\displaystyle=\frac{\sqrt{2}d_{y}L_{r}B}{\sqrt{m}}2^{d-1}\Pi_{i=1}^{d-1}\min(m_{i},n_{i})$
(51)

∎

Finally, we develop the following theorem for the generalization error bound
of LCNNs using SpectralDense layers as an immediate result of Corollary 2.

######  Theorem 8.

Suppose that $\mathcal{X}\subset\mathbb{R}^{d_{x}}$ is a bounded subset such
that for all $x\in\mathcal{X},\lVert x\rVert\leq B$. Let $\mathcal{G}$ be the
hypothesis class of loss functions defined in Eq. 21, where $\mathcal{H}$ is
the hypothesis class of functions from $\mathcal{X}$ to
$\mathcal{Y}\subset\mathbb{R}^{d_{y}}$ with each component function
$h_{k}\in\mathcal{H}^{SD}_{d}$ for $k=1,...,d_{y}$, and
$L:\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}^{\geq 0}$ is a $L_{r}$-Lipschitz
loss function that satisfies the properties in Section 5.1. For any set of $m$
i.i.d. training samples $S=\\{s_{1},s_{2},\cdots,s_{m}\\}$ with
$s_{i}=(x_{i},y_{i})$, $i\in[1,m]$ drawn from a probability distribution
$\mathcal{D}\subset\mathcal{X}\times\mathcal{Y}$, then for any
$\delta\in(0,1)$, the following upper bound holds with probability $1-\delta$
:

$\underset{(x,y)\sim\mathcal{D}}{E}\,[L(h(x),y)]\leq\frac{1}{m}\sum_{i=1}^{m}L(h(x_{i}),y_{i})+\frac{\sqrt{2}d_{y}L_{r}B}{\sqrt{m}}2^{d}\,\Pi_{i=1}^{d-1}\min(m_{i},n_{i})+3M\sqrt{\frac{\log\frac{1}{\delta}}{2m}}$
(52)

where $m_{i}$ and $n_{i}$ are the number of rows and columns in $W_{i}$.

###### Proof.

The proof of Eq. 52 follows immediately from substituting Eq. 50 (i.e., the
upper bound for $\mathcal{R}_{S}(\mathcal{G})$) from Corollary 2 into Eq. 22
in Theorem 3. ∎

###  5.3 Comparison of Empirical Rademacher Complexity Bounds

The bound in Corollary 1 shows that the ERC of LCNNs using SpectralDense
layers is simply bounded by a constant that depends on the size of the network
(i.e., the number of neurons in each layer and the depth of the network).
However, if we consider the ERC of the class of conventional FNNs with
1-Lipschitz activation functions of depth $d$, denoted by
$\mathcal{R}_{S}(\mathcal{H}^{D}_{d})$, we have the following bound by
Golowich et al. (2018) :

$\mathcal{R}_{S}(\mathcal{H}^{D}_{d})\leq\frac{B(\sqrt{2\log(2)d})\,\Pi_{i=1}^{d}R_{i}}{\sqrt{m}}$
(53)

A slight advantage of this bound is that we have an $O(\sqrt{d})$ dependence
on the depth, while the bound in Corollary 1 has an exponential dependence on
$d$. However, the main advantage of the bound in Corollary 1 is that it is
only dependent on the size of the network, i.e., the number of neurons in each
layer. For conventional FNNs with no constraints on the norm of the weight
matrices, even though the training error could be rendered sufficiently small,
the term $\Pi_{i=1}^{d}R_{i}$ is not controlled, which could result in a
sufficiently large ERC for the class $\mathcal{H}^{D}_{d}$, and subsequently a
large generalization error.

Furthermore, given an FNN using conventional dense layers and an LCNN using
SpectralDense layers of the same size parameters, we have demonstrated that
the ERC of LCNNs using SpectralDense layers is bounded by a fixed constant,
while there is no constraint on the ERC bound of the FNNs using dense layers.
This is a tremendous advantage since all the parameters in Eq. 52 in Theorem 8
are fixed except the training sample size $m$ as long as our target function
is Lipschitz continuous and the neural network architecture is fixed. This
implies that there is a provably correct probabilistic guarantee that the
generalization error is at most $O(\frac{1}{\sqrt{m}})$. Therefore, it is
demonstrated that LCNNs using SpectralDense layers can not only approximate a
wide variety of nonlinear functions (i.e., the universal approximation theorem
in Section 4), but also effectively mitigate the issue of over-fitting to data
noise and exhibit better generalization properties, since LCNNs improve the
robustness against data noise as compared to conventional FNNs and are
developed with reduced model complexity.

##  6 Application of LCNNs to Predictive Control of a Chemical Process

In order to show the robustness of LCNNs, we develop LCNNs using SpectralDense
layers and incorporate the LCNNs into a model predictive controller (MPC). We
will demonstrate that LCNNs using SpectralDense layers are able to accurately
model the nonlinear dynamics of a chemical process, and furthermore, they are
able to prevent over-fitting in the presence of data noise.

###  6.1 Chemical Process Description

The chemical process in consideration, which was first developed by Wu et al.
(2019a), is a non-isothermal, well-mixed CSTR that contains an exothermic
second-order irreversible reaction that converts molecule $A$ to molecule $B$.
There is a feed stream of $A$ into the CSTR, as well as a jacket that either
supplies heat or cools the CSTR. We denote $C_{A}$ to be the concentration of
$A$ and $T$ the temperature of the CSTR. The CSTR dynamics can be modeled with
the governing equations as follows:

$\begin{split}\frac{dC_{A}}{dt}&=\frac{F}{V}(C_{A0}-C_{A})-k_{0}e^{-\frac{E}{RT}}C_{A}^{2}\\\
\frac{dT}{dt}&=\frac{F}{V}(T_{0}-T)-\frac{\Delta
H}{\rho_{L}C_{p}}k_{0}e^{-\frac{E}{RT}}C_{A}^{2}+\frac{Q}{\rho_{L}C_{p}V}\end{split}$
(54)

where $F$ is the feed flowrate, $\Delta H$ is the molar enthalpy of reaction,
$k_{0}$ is the rate constant, $R$ is the ideal gas constant, $\rho_{L}$ is the
fluid density, $E$ is the activation energy, and $C_{p}$ is the specific heat
capacity. $C_{A0}$ is the feed flow concentration of $A$, and $Q$ is the rate
at which heat is transferred to the CSTR. The values of the process parameters
used in this work are omitted, as they are exactly the same as those in Wu et
al. (2019b). The CSTR of Eq. 54 has an unstable steady-state
$[C_{As},T_{s},C_{A0_{s}},Q_{s}]=[1.95\leavevmode\nobreak\ \text{mol /
dm}^{3},402\leavevmode\nobreak\ \text{K},4\leavevmode\nobreak\ \text{mol /
dm}^{3},0\leavevmode\nobreak\ \text{kJ / h}]$. We use $x$ to denote the system
states: $x^{T}=[C_{A}-C_{As},T-T_{s}]$ and $u$ to denote the manipulated
inputs: $u^{T}=[C_{A0}-C_{A0_{s}},Q-Q_{s}]$, so that the equilibrium point
$(x_{s},u_{s})$ is located at the origin. Following the formulation presented
in Section 2.2, after extensive numerical simulations, a Lyapunov function
$V(x):=x^{T}Px$ with $P=\left[\begin{array}[]{cc}1060&22\\\
22&0.52\end{array}\right]$ was constructed. The region $\Omega_{\rho}$ with
$\rho=372$ is found via open-loop simulations by sweeping through many
feasible initial conditions within the domain space $D$, such that for any
initial state $x_{0}\in\Omega_{\rho}$, the controller $\Phi$ renders the
origin of the state space exponentially stable.

###  6.2 Data Generation

In order to train the LCNNs, we conducted open-loop simulations of the CSTR
dynamics shown in (54) using many different possible control actions to
generate the required dataset. Specifically, we performed a sweep of all
possible initial states $x_{0}\in\Omega_{\rho}$ and control actions $u\in U$.
The forward Euler method with step size $h_{c}=10^{-5}\leavevmode\nobreak\
\text{hr}$ was used to obtain the value of $\tilde{F}(x_{0},u)$, that is, to
deduce the state after $\Delta=10^{-3}\leavevmode\nobreak\ \text{hr}$ time has
passed. The inputs in the dataset will be all such possible pairs $(x_{0},u)$
and the outputs will be $\tilde{F}(x_{0},u)$. We collected 20000 input-output
pairs in the dataset, and split it into training (52.5 %), validation (17.5
%), and testing (30 %) datasets. Before the training process, the dataset was
pre-processed appropriately by standard scalers to ensure that each variable
has a variance of the same order of magnitude.

###  6.3 Model Training With Noise-Free Data

For the prediction model $\tilde{F}_{nn}$ that approximates the function
$\tilde{F}$, we first demonstrate that in the absence of training data noise,
the LCNNs using SpectralDense layers can capture the dynamics of
$\tilde{F}_{nn}$ well in the operating region $\Omega_{\rho}$. Specifically,
we used an LCNN using two SpectralDense layers with 40 neurons each, followed
by a dense layer with weights that have absolute values bounded by 1.0 and
linear activation functions as the final layer. The weight bound was
implemented using the max_norm constraint function provided by Tensorflow. The
neural network was trained using the Adam Optimizer package provided by
Tensorflow. The loss function that was utilized to train the neural networks
was the mean squared error (MSE), and the testing error for the LCNN reached
$2.83\times 10^{-5}$, which was considered sufficiently small using normalized
training data.

###  6.4 Lyapunov-based MPC

Subsequently, we incorporate the LCNN model into the design of Lyapunov-based
MPC (LMPC). The control actions in LMPC are implemented using the sample-and-
hold method, where $\Delta$ is the sampling period. Let $N$ be a positive
integer that represents the prediction horizon of the MPC. The state at time
$t=k\Delta$ is denoted by $x_{k}$, and the LMPC using LCNN model corresponds
to the optimization problem described below, following the notation in Wu et
al. (2019a):

$\displaystyle\mathcal{J}$
$\displaystyle=\min_{\tilde{u}}\sum_{i=1}^{N}L(\tilde{x}_{k+i},u_{k+i-1})$
(55a) s.t.
$\displaystyle\tilde{x}_{t+1}=\tilde{F}_{nn}(\tilde{x}_{t},u_{t}),\leavevmode\nobreak\
\leavevmode\nobreak\ \forall t\in[k,k+N)$ (55b) $\displaystyle u_{t}\in
U,\leavevmode\nobreak\ \leavevmode\nobreak\ \forall t\in[k,k+N)$ (55c)
$\displaystyle\tilde{x}_{k}=x_{k}$ (55d) $\displaystyle
V(\tilde{F}_{nn}(x_{k},u_{k}))\leq
V(\tilde{F}_{nn}(x_{k},\Phi(x_{k}))),\leavevmode\nobreak\ \leavevmode\nobreak\
\text{ if }\>x_{k}\in\Omega_{\rho}\setminus\Omega_{\rho_{nn}}$ (55e)
$\displaystyle V(x_{t})\leq\rho_{nn}\>\>\forall
t\in[k,k+N],\leavevmode\nobreak\ \leavevmode\nobreak\ \text{if
}\>x_{k}\in\Omega_{\rho_{nn}}$ (55f)

where $\tilde{u}=[u_{k},u_{k+1},u_{k+2},....,u_{k+N-1}]$, and
$\Omega_{\rho_{nn}}$ is a much smaller sublevel set than $\Omega_{\rho}$. The
state predicted by the LCNN model $\tilde{F}_{nn}(\tilde{x}_{t},u_{t})$ is
represented by $\tilde{x}$. Eq. 55c describes the input constraints imposed,
since we assume that $U$ is the set of all possible input constraints. The
initial condition of the prediction model is obtained from the feedback state
measurement at $t=k\Delta$, as shown in Eq. 55d. The constraints of Eqs.
55e-55f guarantee closed-loop stability, i.e., it ensures that the state which
is originally in the set $\Omega_{\rho}$ will eventually converge to the much
smaller sublevel set $\Omega_{\rho_{nn}}$, provided that a sufficiently small
$\Delta>0$ is used such that the controller $\Phi$ when applied with the
sample-and-hold method still guarantees convergence to the origin. Note that
the key difference between the LMPC of Eq. 55 in this work and the one in our
previous work Wu et al. (2019a) is that the LCNN model, which is a type of
feedforward neural network, is used as the underlying model for prediction of
future states in Eq. 55b to predict the state one sampling time forward, while
in Wu et al. (2019a), a recurrent neural network was developed to predict the
trajectory of future states within one sampling period. Therefore, the
formulated objective function shown in Eq. 55a and the constraints of Eqs.
55e-55f only account for the predicted states in the sampling instance. The
above optimization problem is solved using IPOPT, which is a package for
solving large-scale nonlinear and non-convex optimization problems. The LMPC
of Eq. 55 for the CSTR example is designed with the following parameters:
$\rho=372$, $\rho_{nn}=2$, $N=2$, and $L(x,u)=x^{T}Q_{1}x+u^{T}Q_{2}u$, where
$Q_{1}=\left[\begin{array}[]{cc}6.25\times 10^{-4}&0\\\ 0&1\end{array}\right]$
and $Q_{2}=\left[\begin{array}[]{cc}0.01&0\\\ 0&4.0\times
10^{-12}\end{array}\right]$.

###  6.5 Closed-Loop Simulation Results

The closed-loop simulation results under LMPC using LCNN models (termed LCNN-
LMPC) are shown in Figure 1 and Figure 2, where the initial condition is
$x=(72\leavevmode\nobreak\ \text{K}\>,\>-1.65\leavevmode\nobreak\ \text{ kmol
/ m}^{3})$. The closed-loop simulation under LMPC using the first-principles
model of Eq. 54 is also carried out as the reference for comparison purposes.
As shown in Figure 1, the state trajectory when then LCNN is used overlaps
with the trajectory when the first-principles model is used in MPC, suggesting
that the neural network modeling error is sufficiently small such that the
state can be driven to the equilibrium point under MPC. The control actions
taken by the two models are similar as shown in Fig. 2, where slight
deviations and oscillations occur under the MPC using LCNN model. The
oscillations in the manipulated input profile could be due to the following
factors: 1) slight model discrepancies between the model and the first-
principles model, and 2) the IPOPT software being trapped within local minima
of the objective function, or a combination of the two factors (in fact, the
first factor might inadvertently cause the second because of irregularities in
the LMPC loss function $L$ when LCNN is used). Nevertheless, the LCNN is able
to drive the state to the origin effectively, which shows that the LCNN can
effectively model the nonlinear dynamics of the CSTR.

######  Remark 6.

The weight bound of 1.0 in the last layer was chosen to ensure that the class
of functions that can be represented using this neural network architecture is
sufficiently large to contain the target function. If a smaller weight bound
is chosen, the target function might not be approximated well, since the
Lipschitz constant of the network does not meet the one for the target
function.

###  6.6 Robustness against Data Noise

As discussed in Section 5, one of the most crucial advantages of LCNNs using
SpectralDense layers is their robustness against data noise and potential
over-fitting during the training process. For example, when the number of
neurons per hidden layer is exceptionally large, the neural network tends to
overestimate the complexity of the problem and ultimately learns the data
noise (see Ke and Liu (2008) and Sheela et al. (2013) for more details on this
phenomenon). Therefore, in this subsection, we will demonstrate that when
Gaussian data noise is introduced into the training datasets, the LCNNs
outperform the FNNs using conventional dense layers (termed “Dense FNNs”).

We followed the data generation process in Section 6.2, and added Gaussian
noise with a standard deviation of 0.1 or 0.2 to the training dataset. To show
the difference between SpectralDense LCNNs and conventional Dense FNNs in
robustness against over-fitting, we trained SpectralDense LCNNs and Dense FNNs
with the same set of hidden layer architectures. Specifically, we used two
hidden layers in each type of neural network with 640 or 1280 neurons each.
The LCNNs were developed with SpectralDense hidden layers, while the
conventional Dense FNNs were developed using the dense layers from Tensorflow
with ReLU activation functions. Throughout the training process for all the
networks, the same training hyperparameters were used, such as the number of
epochs, the early stopping callback, and the batch size used in the Adam
Optimizer.

Table 1 shows the testing errors for the various neural networks trained. As
seen in Table 1, the testing errors of the conventional Dense FNNs have an
order of magnitude of $10^{-3}$ to $10^{-2}$, which are significantly larger
than those of SpectralDense LCNNs (i.e., the testing errors of SpectralDense
LCNNs have an order of magnitude of $10^{-5}$ to $10^{-4}$). The increase in
testing error in the Dense FNNs is due to over-fitting of the noise since the
testing error has the same order of magnitude as the variance of the Gaussian
noise.

Additionally, we integrated the LCNN and Dense FNN with 640 neurons per layer
and 0.1 standard deviation Gaussian Noise into MPC, similar to the process
described in Section 6.5. The results when both the LCNNs and Dense FNNs are
integrated into MPC are shown in Figure 3 and 4. From the plot of the Lyapunov
function value $V(x)$ in Figure 3, it is demonstrated that the Dense FNN is
unable to effectively drive the state to the origin compared to the LCNN. This
is readily observed because not only is the Lyapunov function value for the
Dense FNN much higher, but there are also considerable oscillations in the
function value, especially in the time frame between $0.15$ hr and $0.3$ hr.
In addition, in Figure 4, it is observed that, while the predicted control
actions of the LCNN are very similar to those of the first-principles model,
the predicted control actions under the Dense FNN show large oscillations and
differ largely from those of the first-principles model. This large disparity
between the Dense FNN and the first-principles model shows that the Dense FNN
has become incapable of accurately modeling the process dynamics when embedded
into MPC.

###  6.7 Comparison of Lipschitz Constants between LCNNs and Dense FNNs

Additionally, we compare the Lipschitz constants of the LCNNs and the Dense
FNNs developed for the CSTR of Eq. 54 and demonstrate that the conventional
Dense FNNs have a much larger Lipschitz constant than the SpectralDense LCNNs
as a result of noise over-fitting and the lack of constraints on the weight
matrices. For the Dense FNNs, we used the LipBaB algorithm to obtain the
Lipschitz constant for the FNNs using dense layers, and for the SpectralDense
LCNNs, we took the SVD of the last weight matrix and obtained the spectral
norm of the last weight matrix as the upper bound of the Lipschitz constant.
The results are shown in Table 2. The Dense FNNs have Lipschitz constants in
the order of magnitude $10^{2}$ compared to SpectralDense LCNNs with Lipschitz
constants in the order of magnitude $10^{0}$. The comparison of Lipschitz
constants demonstrates that the LCNNs are also provably and certainly less
sensitive to input perturbations as compared to the Dense FNNs since they have
a Lipschitz constant several orders of magnitude lower. Furthermore, the
calculation of Lipschitz constants demonstrates that LCNNs are able to prevent
over-fitting data noise by maintaining a small Lipschitz constant, while
conventional Dense FNNs with a large number of neurons could over-fit data
noise.

##  7 Conclusions

In this work, we developed LCNNs for the general class of nonlinear systems,
and discussed how LCNNs using SpectralDense layers can mitigate sensitivity
issues and prevent over-fitting to noisy data from the perspectives of
Lipschitz constants and generalization error. Specifically, we first proved
the universal approximation theorem for LCNNs using SpectralDense layers to
demonstrate that LCNNs are capable of retaining expressive power for Lipschitz
target functions despite having a small hypothesis class. Then, we derived the
generalization error bound for SpectralDense LCNNs using the Rademacher
complexity method. The above results provided the theoretical foundations to
demonstrate that LCNNs can improve input sensitivity due to their constrained
Lipschitz constant and generalize better to prevent over-fitting. Finally,
LCNNs using SpectralDense layers were integrated into MPC and applied to a
chemical reactor example. The simulations show that the LCNNs effectively
captured the process dynamics and outperformed conventional Dense FNNs in
terms of smaller testing errors, higher prediction accuracy in MPC, and
smaller Lipschitz constants in the presence of noisy training data.

##  8 Acknowledgments

Financial support from the NUS Start-up grant R-279-000-656-731 is gratefully
acknowledged.

## References

  * Anil et al. (2019) Anil, C., Lucas, J., Grosse, R., 2019.  Sorting out lipschitz function approximation, in: International Conference on Machine Learning, PMLR. pp. 291–301.
  * Balda et al. (2019) Balda, E.R., Behboodi, A., Mathar, R., 2019.  Perturbation analysis of learning algorithms: Generation of adversarial examples from classification to regression.  IEEE Transactions on Signal Processing 67, 6078–6091.
  * Baldi and Sadowski (2014) Baldi, P., Sadowski, P., 2014\.  The dropout learning algorithm.  Artificial intelligence 210, 78–122.
  * Baldi and Sadowski (2013) Baldi, P., Sadowski, P.J., 2013\.  Understanding dropout.  Advances in neural information processing systems 26\.
  * Bhowmick et al. (2021) Bhowmick, A., D’Souza, M., Raghavan, G.S., 2021.  LipBaB: Computing exact lipschitz constant of relu networks, in: Artificial Neural Networks and Machine Learning–ICANN 2021: 30th International Conference on Artificial Neural Networks, Bratislava, Slovakia, September 14–17, 2021, Proceedings, Part IV 30, Springer. pp. 151–162.
  * Chen et al. (1995) Chen, H., Kremling, A., Allgöwer, F., 1995.  Nonlinear predictive control of a benchmark CSTR, in: Proceedings of 3rd European control conference, sn. pp. 3247–3252.
  * Cortes et al. (2012) Cortes, C., Mohri, M., Rostamizadeh, A., 2012.  L2 regularization for learning kernels.  arXiv preprint arXiv:1205.2653 .
  * Federer (2014) Federer, H., 2014.  Geometric measure theory.  Springer.
  * Golowich et al. (2018) Golowich, N., Rakhlin, A., Shamir, O., 2018.  Size-independent sample complexity of neural networks, in: Conference On Learning Theory, PMLR. pp. 297–299.
  * Goodfellow et al. (2014) Goodfellow, I.J., Shlens, J., Szegedy, C., 2014.  Explaining and harnessing adversarial examples.  arXiv preprint arXiv:1412.6572 .
  * Ke and Liu (2008) Ke, J., Liu, X., 2008.  Empirical analysis of optimal hidden neurons in neural network modeling for stock prediction, in: Proceedings of 2008 IEEE Pacific-Asia Workshop on Computational Intelligence and Industrial Application, IEEE. pp. 828–832.
  * Kwon et al. (2013) Kwon, J.S.I., Nayhouse, M., Christofides, P.D., Orkoulas, G., 2013\.  Modeling and control of protein crystal shape and size in batch crystallization.  AIChE Journal 59, 2317–2327.
  * Kwon et al. (2014) Kwon, J.S.I., Nayhouse, M., Christofides, P.D., Orkoulas, G., 2014\.  Protein crystal shape and size control in batch crystallization: Comparing model predictive control with conventional operating policies.  Industrial & Engineering Chemistry Research 53, 5002–5014.
  * Madry et al. (2017) Madry, A., Makelov, A., Schmidt, L., Tsipras, D., Vladu, A., 2017.  Towards deep learning models resistant to adversarial attacks.  arXiv preprint arXiv:1706.06083 .
  * Maurer (2016) Maurer, A., 2016.  A vector-contraction inequality for rademacher complexities, in: Algorithmic Learning Theory: 27th International Conference, ALT 2016, Bari, Italy, October 19-21, 2016, Proceedings 27, Springer. pp. 3–17.
  * Mohri et al. (2018) Mohri, M., Rostamizadeh, A., Talwalkar, A., 2018.  Foundations of machine learning.  MIT press.
  * Moore and DeNero (2011) Moore, R.C., DeNero, J., 2011\.  L1 and L2 regularization for multiclass hinge loss models, in: Symposium on Machine Learning in Speech and Natural Language Processing.
  * Neyshabur et al. (2015) Neyshabur, B., Tomioka, R., Srebro, N., 2015.  Norm-based capacity control in neural networks, in: Conference on learning theory, PMLR. pp. 1376–1401.
  * Sabiri et al. (2022) Sabiri, B., El Asri, B., Rhanoui, M., 2022.  Mechanism of overfitting avoidance techniques for training deep neural networks., in: ICEIS (1), pp. 418–427.
  * Serrurier et al. (2021) Serrurier, M., Mamalet, F., González-Sanz, A., Boissin, T., Loubes, J.M., Del Barrio, E., 2021\.  Achieving robustness in classification using optimal transport with hinge regularization, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 505–514.
  * Shaham et al. (2015) Shaham, U., Yamada, Y., Negahban, S., 2015.  Understanding adversarial training: Increasing local stability of neural nets through robust optimization.  arXiv preprint arXiv:1511.05432 .
  * Sheela et al. (2013) Sheela, K.G., Deepa, S.N., et al., 2013\.  Review on methods to fix number of hidden neurons in neural networks.  Mathematical problems in engineering 2013\.
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., Salakhutdinov, R., 2014.  Dropout: a simple way to prevent neural networks from overfitting.  The journal of machine learning research 15, 1929–1958.
  * Szegedy et al. (2013) Szegedy, C., Zaremba, W., Sutskever, I., Bruna, J., Erhan, D., Goodfellow, I., Fergus, R., 2013.  Intriguing properties of neural networks.  arXiv preprint arXiv:1312.6199 .
  * Wu (2001) Wu, F., 2001.  Lmi-based robust model predictive control and its application to an industrial CSTR problem.  Journal of process control 11, 649–659.
  * Wu et al. (2022) Wu, Z., Alnajdi, A., Gu, Q., Christofides, P.D., 2022.  Statistical machine-learning–based predictive control of uncertain nonlinear processes.  AIChE Journal 68, e17642.
  * Wu et al. (2021a) Wu, Z., Luo, J., Rincon, D., Christofides, P.D., 2021a.  Machine learning-based predictive control using noisy data: evaluating performance and robustness via a large-scale process simulator.  Chemical Engineering Research and Design 168, 275–287.
  * Wu et al. (2021b) Wu, Z., Rincon, D., Gu, Q., Christofides, P.D., 2021b.  Statistical machine learning in model predictive control of nonlinear processes.  Mathematics 9, 1912\.
  * Wu et al. (2019a) Wu, Z., Tran, A., Rincon, D., Christofides, P.D., 2019a.  Machine learning-based predictive control of nonlinear processes. part i: theory.  AIChE Journal 65, e16729.
  * Wu et al. (2019b) Wu, Z., Tran, A., Rincon, D., Christofides, P.D., 2019b.  Machine-learning-based predictive control of nonlinear processes. part ii: Computational implementation.  AIChE Journal 65, e16734.
  * Ying (2019) Ying, X., 2019.  An overview of overfitting and its solutions.  Journal of Physics: Conference Series 1168, 022022.

Figure 1: Closed-loop state profiles of $T$ and $C_{A}$, and the evolution of the Lyapunov function $V(x)$ under the LMPCs using LCNN (red solid line) and first-principles model (green dashed line), respectively. Figure 2: Closed-loop manipulated input profiles of $C_{A0}$ and $Q$ under the LMPCs using LCNN (red solid line) and first-principles model (green dashed line), respectively. Figure 3: Closed-loop state profiles of $T$ and $C_{A}$, and the evolution of the Lyapunov function $V(x)$ under the LMPCs using LCNN (red solid line), Dense FNN (orange dotted line), and first-principles model (green dashed line), respectively, where noisy training data is used to develop LCNN and Dense FNN. Figure 4: Closed-loop manipulated input profiles of $C_{A0}$ and $Q$ under the LMPCs using LCNN (red solid line), Dense FNN (orange dotted line), and first-principles model (green dashed line), respectively, where noisy training data is used to develop LCNN and Dense FNN.  Hidden Layers | Noise SD | LCNN Testing Error | Dense Testing Error | Improvement Factor
---|---|---|---|---
(640,640) | 0.1 | $3.617\times 10^{-5}$ | $4.573\times 10^{-3}$ | $1.264\times 10^{2}$
(1280,1280) | 0.1 | $3.850\times 10^{-5}$ | $6.928\times 10^{-3}$ | $1.799\times 10^{2}$
(640,640) | 0.2 | $2.088\times 10^{-5}$ | $1.692\times 10^{-2}$ | $8.102\times 10^{1}$
(1280,1280) | 0.2 | $1.393\times 10^{-5}$ | $3.095\times 10^{-2}$ | $2.220\times 10^{2}$
Table 1: Comparison of the testing errors for various hidden layer architectures and standard deviation (SD) of noise introduced into the training dataset. The improvement factor is the ratio between the two testing errors.  Hidden Layers | Noise SD | Dense Lipschitz constant | LCNN Lipschitz constant
---|---|---|---
(640,640) | 0.1 | $2.447\times 10^{2}$ | 1.119
(1280,1280) | 0.1 | $2.682\times 10^{2}$ | 1.116
(640,640) | 0.2 | $2.511\times 10^{2}$ | 1.114
(1280,1280) | 0.2 | $9.285\times 10^{2}$ | 1.107
Table 2: Lipschitz constants for the various hidden layer architectures and
standard deviation (SD) of noise introduced into the training dataset.