Robust Machine Learning Modeling for Predictive Control Using
  Lipschitz-Constrained Neural Networks
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.
Wallace Tan Gian Yion, Zhe Wu
# Robust Machine Learning Modeling for Predictive Control Using Lipschitz-Constrained Neural Networks

###### 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 +
Footnote †: journal: Journal of Computer Vision

## 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. (2021), 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

### Notations

\(\|W\|_{F}\) and \(\|W\|_{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}\rightarrow\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 \(\|x\|\) denote the Euclidean norm of \(x\). A metric space is a set \(X\) equipped with a metric function \(d_{X}:X\times X\rightarrow\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.
### 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 \tag{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}\rightarrow\mathbb{R}^{n}\), \(g:\mathbb{R}^{n}\rightarrow\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\rightarrow\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\mid 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}\):

\[\|F(x,u)\|\leq K_{F} \tag{2a}\] \[\|F(x^{\prime},u)-F(x,u)\|\leq L_{x}\|x^{\prime}-x\| \tag{2b}\]

### 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-Lindelof 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 \tag{3}\]
We can then define \(\tilde{F}:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\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

\[\|f(x)-f(y)\|\leq L\cdot\|x-y\|\]

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.

### 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},\ldots,x_{m})=(\sigma^{\prime}(x_{1}),\sigma^{\prime}(x_{2 }),\ldots,\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}\rightarrow\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} \tag{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} \tag{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 \(\|W\|_{2}\) is also 1. The function \(\sigma:\mathbb{R}^{m}\rightarrow\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}_{n}^{m}\) be the class of Lipschitz-constrained neural networks (LCNNs) as follows:

\[\begin{array}{l}\mathcal{LN}_{n}^{m}:=\{\,f\ |\ f:\mathbb{R}^{n}\rightarrow \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 }\|W_{i}\|_{2}=1,i=1,...,j\,\} \end{array} \tag{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}\|J_{f}(x)\|_{2} \tag{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.

### 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}_{n}^{m}\), 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}_{n}^{m}\) 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\mid f:\mathbb{R}^{n}\rightarrow\mathbb{R},\, \exists j\in\mathbb{N},s.t.\,f=f_{j}\circ...f_{2}\circ f_{1},f_{i}=\sigma(W_{i}x +b),\|W_{i}\|_{2}=1,i=1,...,j\} \tag{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\rightarrow\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\mid|c|\leq k,\;f\in A\} \tag{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 \tag{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}|f(x)-\tilde{f}(x)|<\epsilon \tag{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\|x-y\|\), 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 \(\|w\|=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\|w\|\cdot\|y-x\|=\|y-x\| \tag{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},\ldots,W_{d_{f}}^{f},b_{d_{f}}^{f}]~{}~ {}[W_{1}^{g},b_{1}^{g},W_{2}^{g},b_{2}^{g},\ldots,W_{d_{g}}^{g},b_{d_{g}}^{g}] \tag{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 \(\|W_{i}^{h}\|_{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}] \tag{13}\]

where \(c\) is a positive constant chosen based on \(W_{1}^{f}\) and \(W_{1}^{g}\) to ensure that \(\|W_{1}^{h}\|_{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} \tag{14}\]
To prove that \(\|W_{i}^{h}\|_{2}=1\), we use the singular value decomposition.

\[W_{i}^{h} =\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} \tag{15a}\] \[=\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) \[=UDV^{T} \tag{15c}\]

It is noted that \(\|W_{i}^{h}\|_{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 \(\|W_{i}^{h}\|_{2}=1\). We also take the following.

\[b_{i}^{h}=c[b_{i}^{f},b_{i}^{g}] \tag{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} \tag{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}\rightarrow\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.

### 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 \(\|x\|\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}\rightarrow\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}\rightarrow\mathcal{Y}\) be a function that represents a neural network model (termed hypothesis) in a hypothesis class.
_Definition 9_.: We define **generalization error** as

\[\mathop{\mathbb{E}}_{(x,y)\sim\mathcal{D}}\left[L(h(x),y)\right]=\int_{\mathcal{X }\times\mathcal{Y}}L(h(x),y)\,\mathbb{P}(dx\times dy) \tag{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}) \tag{19}\]

The generalization error is a measure of how well a hypothesis \(h:\mathcal{X}\rightarrow\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})=\|y-y^{\prime}\|^{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.

### 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 **empiricial Rademacher complexity** (ERC) of \(S\) with respect to \(H\), denoted by \(\mathcal{R}_{S}(H)\), is defined as:

\[\mathcal{R}_{S}(H):=\mathop{\mathbb{E}}_{\epsilon}\sup_{h\in H}\frac{1}{m} \sum_{i=1}^{m}\epsilon_{i}h(x_{i}) \tag{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}}\rightarrow\mathcal{Y}\subset\mathbb{R}^{d _{y}}\) and \(L:\mathcal{Y}\times\mathcal{Y}\rightarrow\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}\} \tag{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}}{\mathrm{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}} \tag{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})=\mathbb{E}\sup_{\epsilon}\frac{1}{n \in\mathcal{H}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}L(h(x_{i}),y)\leq\sqrt{2} L_{r}\mathbb{E}\sup_{\epsilon}\frac{1}{m}\sum_{i=1}^{m}\sum_{k=1}^{d_{y}} \epsilon_{ik}h_{k}(x_{i}) \tag{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:

\[\mathbb{E}\sup_{\epsilon}\frac{1}{m}\sum_{i=1}^{m}\sum_{k=1}^{d_{y}}\epsilon_ {ik}h_{k}(x_{i})\leq\sum_{k=1}^{d_{y}}\mathbb{E}\sup_{\epsilon}\frac{1}{m} \sum_{i=1}^{m}\epsilon_{ik}h_{k}(x_{i}) \tag{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))) \tag{25}\]

where each of the weight matrices \(W_{i}\in\mathbb{R}^{m_{i}\times n_{i}}\) has a bounded Frobenius norm, that is, \(\|W_{i}\|_{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}_{d}^{SD}\) 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))) \tag{26}\]

where each of the weight matrices satisfies \(\|W_{i}\|_{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}_{d}^{SD}\) is that \(\mathcal{H}_{d}\) is defined for any FNN using GroupSort activation functions, while \(\mathcal{H}_{d}^{SD}\) is defined for LCNNs that use the GroupSort activation function and satisfy \(\|W_{i}\|_{2}=1\). Additionally, the Frobenius norm is used in \(\mathcal{H}_{d}\), while the spectral norm is used in \(\mathcal{H}_{d}^{SD}\). Despite the differences between \(\mathcal{H}_{d}\) and \(\mathcal{H}_{d}^{SD}\), it will be demonstrated from the following lemmas and theorems that the two classes \(\mathcal{H}_{d}\) and \(\mathcal{H}_{d}^{SD}\) 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}_{d}^{SD}\) 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. (2021, 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|) \tag{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))) \tag{28}\]

where each of the weight matrices \(W_{i}\in\mathbb{R}^{m_{i}\times n_{i}}\) has a bounded Frobenius norm, that is, \(\|W_{i}\|_{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 \(\|W_{d}\|_{F}\leq R_{d}\). For any dataset with \(m\) data points, we have the following inequality:

\[\operatorname*{\mathbb{E}}_{\epsilon}\sup_{h\in\tilde{\mathcal{H}}_{d}}\left\| \frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\|\leq 2R_{d}\operatorname*{ \mathbb{E}}_{\epsilon}\sup_{h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\left\| \sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\| \tag{29}\]

Proof.: Letting \(w_{1}^{T},w_{2}^{T},\cdots,w_{k}^{T}\) represent the rows of \(W_{d}\), we have

\[\operatorname*{\mathbb{E}}_{\epsilon}\sup_{h\in\tilde{\mathcal{H }}_{d}}\left\|\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\| =\operatorname*{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{F}\leq R_{ d},h\in\tilde{\mathcal{H}}_{d-1}}\left\|\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i} \sigma(W_{d}h(x_{i}))\right\| \tag{30a}\] \[=\operatorname*{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{F}\leq R_ {d},h\in\tilde{\mathcal{H}}_{d-1}}\left\|\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\| \tag{30b}\]

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

\[\operatorname*{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{F}\leq R_{d},h\in\tilde {\mathcal{H}}_{d-1}}\left\|\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})| \end{bmatrix}\right\| \tag{31}\]
We can bound Eq. 31 using the triangle inequality with \(A_{1}\) and \(A_{2}\) defined as follows.

\[A_{1}= \mathbb{E}\sup_{\epsilon}\limits_{\|W_{d}\|_{F}\leq R_{d},h\in \tilde{\mathcal{H}}_{d-1}}\left\|\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\|+\mathbb{E}\sup_{\epsilon}\limits_{\|W_{d }\|_{F}\leq R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\left\|\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\| \tag{32a}\] \[A_{2}=\mathbb{E}\sup_{\|W_{d}\|_{F}\leq R_{d},h\in\tilde{ \mathcal{H}}_{d-1}}\left\|\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})|\\ \end{bmatrix}\right\| \tag{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}=\mathbb{E}\sup_{\epsilon}\limits_{\|W_{d}\|_{F}\leq R_{d},h\in\tilde{ \mathcal{H}}_{d-1}}\left\|\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})\\ \end{bmatrix}\right\|=\mathbb{E}\sup_{\epsilon}\limits_{\|W_{d}\|_{F}\leq R_{d },h\in\tilde{\mathcal{H}}_{d-1}}\left\|\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i} W_{d}h(x_{i})\right\| \tag{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}\bigg{(}\sum_{i=1}^{m}\bigg{(} \epsilon_{i}\frac{w_{j}^{T}}{\lVert w_{j}\rVert}h(x_{i})\bigg{)}\bigg{)}^{2}} \tag{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

\[\mathbb{E}\sup_{\epsilon}\limits_{\|W_{d}\|_{F}\leq R_{d},h\in\tilde{ \mathcal{H}}_{d-1}}\left\|\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})\\ \end{bmatrix}\right\| =\mathbb{E}\sup_{\epsilon}\limits_{\|w\|\leq R_{d},h\in\tilde{ \mathcal{H}}_{d-1}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}w^{T}h(x_{i})\] \[\leq\mathbb{E}\sup_{\epsilon}\limits_{\|w\|\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{\|}\]
\[=R_{d}\mathop{\mathbb{E}}_{\epsilon}\sup_{h\in\widetilde{\mathcal{H}}_{d-1}} \frac{1}{m}\bigg{\|}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\bigg{\|} \tag{35}\]

The second term \(A_{2}\) can be rewritten in the following way.

\[A_{2}=\mathop{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{F}\leq R_{d},h\in \widetilde{\mathcal{H}}_{d-1}}\left\|\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})|\\ \end{bmatrix}\right\| \tag{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 \(\|W_{d}\|_{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\|w_{2j}+w_{2j-1}\|^{2}\bigg{(}\sum_{i=1}^{ m}\epsilon_{i}\Big{|}\frac{(w_{2j}+w_{2j-1})^{T}}{\|w_{2j}+w_{2j-1}\|}h(x_{i}) \bigg{|}\bigg{)}^{2}} \tag{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{(}\|w_{2j}\|^{2}+\|w_{2j-1}\|^{2} \bigg{)}\bigg{(}\sum_{i=1}^{m}\epsilon_{i}\bigg{|}\frac{(w_{2j}+w_{2j-1})^{T}} {\|w_{2j}+w_{2j-1}\|}h(x_{i})\bigg{|}\bigg{)}^{2}} \tag{38}\]

It is readily shown that Eq. 38 is maximized when \(\|w_{2j}\|^{2}+\|w_{2j-1}\|^{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}\bigg{(}\sum_{i=1}^{m}\epsilon_{i}\bigg{|}\frac{( w_{2j}+w_{2j-1})^{T}}{\|w_{2j}+w_{2j-1}\|}h(x_{i})\bigg{|}\bigg{)}^{2}}=\frac{R_{d} }{m}\bigg{|}\sum_{i=1}^{m}\epsilon_{i}\bigg{|}\frac{(w_{1}+w_{2})^{T}}{\|w_{1 }+w_{2}\|}h(x_{i})\bigg{|}\bigg{|} \tag{39}\]

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

\[A_{2} =\mathop{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{F}\leq R_{d},h\in \widetilde{\mathcal{H}}_{d-1}}\left\|\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})|\\ \end{bmatrix}\right\|\] \[\leq\mathop{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{F}\leq R_{d},h \in\widetilde{\mathcal{H}}_{d-1}}\frac{R_{d}}{m}\bigg{|}\sum_{i=1}^{m}\epsilon _{i}\bigg{|}\frac{(w_{1}+w_{2})^{T}}{\|w_{1}+w_{2}\|}h(x_{i})\bigg{|}\bigg{|}\] \[\leq R_{d}\mathop{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{F}\leq R_ {d},h\in\widetilde{\mathcal{H}}_{d-1}}\frac{1}{m}\bigg{|}\sum_{i=1}^{m}\epsilon _{i}\frac{(w_{1}+w_{2})^{T}}{\|w_{1}+w_{2}\|}h(x_{i})\bigg{|}\qquad\text{by Talagrand's Contraction Lemma}\]
\[\leq R_{d}\mathop{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{F}\leq R_{d},h\in\tilde{\mathcal{H}}_{d-1}}\frac{1}{m}\left\|\frac{(w_{1}+w_{2})^{T}}{\|w_{ 1}+w_{2}\|}\right\|\left\|\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\|\] \[\leq R_{d}\mathop{\mathbb{E}}_{\epsilon}\sup_{h\in\tilde{\mathcal{ H}}_{d-1}}\frac{1}{m}\left\|\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\| \tag{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.

\[\mathop{\mathbb{E}}_{\epsilon}\sup_{h\in\tilde{\mathcal{H}}_{d}} \left\|\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\| \leq A_{1}+A_{2}\] \[\leq R_{d}\mathop{\mathbb{E}}_{\epsilon}\sup_{h\in\tilde{ \mathcal{H}}_{d-1}}\frac{1}{m}\left\|\sum_{i=1}^{m}\epsilon_{i}h(x_{i}) \right\|+R_{d}\mathop{\mathbb{E}}_{\epsilon}\sup_{h\in\tilde{\mathcal{H}}_{d- 1}}\frac{1}{m}\left\|\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\|\] \[=2R_{d}\mathop{\mathbb{E}}_{\epsilon}\sup_{h\in\tilde{\mathcal{H} }_{d-1}}\frac{1}{m}\left\|\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\| \tag{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},\|x\|\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} \tag{42}\]

_where \(\|W_{i}\|_{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})=\mathop{\mathbb{E}}_{\epsilon}\sup_{h\in \mathcal{H}_{d}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})=\mathop{ \mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|\leq R_{d},h\in\tilde{\mathcal{H}}_{d-1}} \frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}W_{d}h(x_{i}) \tag{43}\]

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

\[\mathop{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{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})\leq \mathop{\mathbb{E}}_{\epsilon}\sup_{\|W_{d}\|_{2}\leq R_{d},h\in\tilde{ \mathcal{H}}_{d-1}}\|W_{d}\|_{2}\left\|\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x _{i})\right\|\]
\[\leq R_{d}\mathop{\mathbb{E}}\limits_{\epsilon}\sup_{h\in\tilde{\mathcal{H}}_{d-1} }\left\|\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\| \tag{44}\]

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

\[R_{d}\mathop{\mathbb{E}}\limits_{\epsilon}\sup_{h\in\tilde{\mathcal{H}}_{d-1} }\left\|\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}h(x_{i})\right\|\leq(2^{d-1}\Pi_{ i=1}^{d}R_{i})\mathop{\mathbb{E}}\limits_{\epsilon}\left\|\frac{1}{m}\sum_{i=1}^{m} \epsilon_{i}x_{i}\right\| \tag{45}\]

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

\[\mathop{\mathbb{E}}\limits_{\epsilon}\left\|\frac{1}{m}\sum_{i=1 }^{m}\epsilon_{i}x_{i}\right\| \leq\frac{1}{m}\sqrt{\mathop{\mathbb{E}}\limits_{\epsilon}\left\| \sum_{i=1}^{m}\epsilon_{i}x_{i}\right\|^{2}}\] \[=\frac{1}{m}\sqrt{\mathop{\mathbb{E}}\limits_{\epsilon}\sum_{i= 1}^{m}\sum_{j=1}^{m}\epsilon_{i}\epsilon_{j}x_{j}^{T}x_{i}}\] \[=\frac{1}{m}\sqrt{\mathop{\mathbb{E}}\limits_{\epsilon}m\|x\|^{ 2}}=\frac{B}{\sqrt{m}} \tag{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} \tag{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:

\[\|W\|_{2}\leq\|W\|_{F}\leq\min(m,n)\|W\|_{2} \tag{48}\]

The proof can be readily obtained based on the fact that \(\|W\|_{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}_{d}^{SD})\), where \(\mathcal{H}_{d}^{SD}\) is the hypothesis class of LCNNs using SpectralDense layers of depth \(d\).
**Corollary 1**.: _Let \(\mathcal{H}_{d}^{SD}\) 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},\|x\|\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}_{d}^{SD})\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}) \tag{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 \(\|W_{i}\|_{2}=1\) for all functions in \(\mathcal{H}_{d}^{SD}\) and Lemma 7 (i.e., \(\|W_{i}\|_{F}\leq\min(m_{i},n_{i})\|W_{i}\|_{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}_{d}^{SD}\). 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},\|x\|\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}_{d}^{SD}\) for \(k=1,...,d_{y}\), and \(L:\mathcal{Y}\times\mathcal{Y}\rightarrow\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}) \tag{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.
\[\mathcal{R}_{S}(\mathcal{G})=\mathbb{E}\sup_{\epsilon}\frac{1}{n} \sum_{h\in\mathcal{H}}^{m}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}L(h(x_{i}),y) \leq\sqrt{2}L_{r}\mathbb{E}\sup_{\epsilon}\frac{1}{m}\sum_{i=1}^{m} \sum_{k=1}^{d_{y}}\epsilon_{ik}h_{k}(x_{i})\] \[\leq\sqrt{2}L_{r}\sum_{k=1}^{d_{y}}\mathbb{E}\sup_{h_{k}\in \mathcal{H}_{d}^{SD}}\frac{1}{m}\sum_{i=1}^{m}\epsilon_{ik}h_{k}(x_{i})\] \[=\sqrt{2}L_{r}d_{y}\,\mathbb{E}\sup_{h_{k}\in\mathcal{H}_{d}^{SD} }\frac{1}{m}\sum_{i=1}^{m}\epsilon_{ik}h_{k}(x_{i})\] \[=\frac{\sqrt{2}d_{y}L_{r}B}{\sqrt{m}}2^{d-1}\Pi_{i=1}^{d-1}\min( m_{i},n_{i}) \tag{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},\|x\|\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}_{d}^{SD}\) for \(k=1,...,d_{y}\), and \(L:\mathcal{Y}\times\mathcal{Y}\rightarrow\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}} \tag{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. 

### 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})}{\sqrt{m}}\, \Pi_{i=1}^{d}R_{i} \tag{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.

### Chemical Process Description

The chemical process in consideration, which was first developed by Wu et al. (2019), 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} \tag{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. (2019). The CSTR of Eq. 54 has an unstable steady-state \([C_{As},T_{s},C_{A0_{s}},Q_{s}]=[1.95\text{ mol }/\text{ dm}^{3},402\text{ K},4\text{ mol }/\text{ dm}^{3},0\text{ kJ }/\text{ 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.

### 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}\) hr was used to obtain the value of \(\tilde{F}(x_{0},u)\), that is, to deduce the state after \(\Delta=10^{-3}\) 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.

### 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.

### 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. (2019):

\[\mathcal{J}=\min_{\tilde{u}}\sum_{i=1}^{N}L(\tilde{x}_{k+i},u_{k+ i-1})\] (55a) s.t. \[\tilde{x}_{t+1}=\tilde{F}_{nn}(\tilde{x}_{t},u_{t}),\ \ \forall t\in[k,k+N) \tag{55b}\] \[u_{t}\in U,\ \ \forall t\in[k,k+N)\] (55c) \[\tilde{x}_{k}=x_{k} \tag{55d}\]
\[V(\tilde{F}_{nn}(x_{k},u_{k}))\leq V(\tilde{F}_{nn}(x_{k},\Phi(x_{k}))), \quad\text{if}\;\;x_{k}\in\Omega_{\rho}\setminus\Omega_{\rho_{nn}} \tag{55e}\] \[V(x_{t})\leq\rho_{nn}\;\;\forall t\in[k,k+N],\;\;\text{if}\;\;x_{k} \in\Omega_{\rho_{nn}} \tag{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. (2019) 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. (2019), 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]\).

### 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\text{ K}\;,\;-1.65\;\text{ kmol}\;/\;\text{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.

### 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.

### 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.