\section{Introduction}
Over the past decades, model-based iterative algorithms have been widely being used for solving linear inverse problems, such as super-resolution \cite{superres}, denoising \cite{denoise}, compressed sensing \cite{irls}, deblurring \cite{deblure}, X-ray computed tomography \cite{ct}.  As opposed to traditional iterative algorithms, deep neural networks (DNNs) have become very popular for tackling such tasks \cite{fdnn}, \cite{wavenet}, \cite{dr2,dagan}, \cite{learndeblur}, \cite{lct}, since DNNs have proven to significantly reduce the time complexity and increase the quality of the reconstruction. A new line of research lies on merging DNNs and iterative algorithms, leading to the so-called \emph{deep unfolding/unrolling} \cite{unfolding,unrolling}. The latter pertains to unfolding the iterations of well-known iterative algorithms into layers of a DNN, which reconstructs the signals of interest.\\
In this paper, we aim at interpreting an optimization-based algorithm for \emph{Compressed Sensing} (CS) \cite{cs} as a deep unfolding network. CS is a modern technique to recover signals of interest $x\in\mathbb{R}^n$ from few linear and possibly corrupted observations $y=Ax+e\in\mathbb{R}^m$, $m<n$. The applications of CS vary among Radar Imaging \cite{radar}, Cryptography \cite{crypto}, Telecommunications \cite{comm}, Magnetic Resonance Imaging \cite{mri}. CS heavily relies on the sparsity of $x$. We call a signal sparse, if it has very few --compared to its dimension-- nonzero entries. Most real-world signals are not naturally sparse, but are considered to be sparse when represented by the elements (\emph{atoms}) of a basis/frame, constituting a so-called \emph{dictionary} \cite{elad}. The choice of the appropriate sparsifying dictionary plays a key role in the reconstruction quality of CS and the optimal sampling rate \cite{genzel}, and depends on the application for which CS is employed. For example, cartoon-like images are proven to be optimally sparsely represented by shearlets \cite{shearlets}, while the sparse structure of audio signals is better captured by time-frequency transforms \cite{multipliers}. Since it is nontrivial to select a sparsifying dictionary for CS, an interesting idea is to combine a dictionary learning technique with a corresponding unfolded network \cite{istagen}.
\subsection{Related work}
Deep unfolding networks performing sparse recovery have gained much attention in the last few years \cite{lista}, \cite{wcs}, \cite{bertocchi}, because of some advantages they have compared to traditional DNNs: they are interpretable, integrate prior knowledge about the signal structure \cite{deligiannis}, and have a relatively small number of trainable parameters \cite{adm}. Especially in the case of CS, many unfolding networks have proven to work particularly well. For example, \cite{amp,ampnet}, \cite{adablocklista,ista-net,tista,tfrnn}, \cite{admm-net}, \cite{fpc} interpret the iterations of approximate message passing (AMP) \cite{ampalgo}, iterative soft-thresholding (ISTA) \cite{istalgo}, alternating direction method of multipliers (ADMM)  \cite{boyd:5} and fixed-point continuation (FPC) \cite{fpcalgo} algorithms, respectively, as layers of a corresponding neural network, which learns a \emph{decoder}; the latter is a function that reconstructs $x$ from $y$. Additionally, the different variants of these networks may jointly learn one, or more than one, from the following: a) step-sizes and/or thresholds used by the original iterative schemes b) the measurement matrix $A$ c) a transform that sparsely represents $x$. The sparsifying transform may either be a square matrix \cite{admm-net}, further constrained to be orthogonal \cite{istagen} -- integrating that way a dictionary learning technique -- or a nonlinear transform, e.g. a combination of linear convolutional operators, separated by a rectified linear unit (ReLU) \cite{ista-net}. Dictionary learning has proven itself very useful when combined with model-based methods for CS \cite{imagedl,iot,1bitcs}, so it looks natural to employ it in deep unfolding networks as well.\\
Although deep unfolding networks have experimentally shown very promising results in sparse recovery, research community focuses lately on the study of their mathematical properties. For example, \cite{theoreticalblocklista,theoreticalista} provide convergence guarantees regarding variants of learned ISTA networks, while \cite{bertocchi} studies the robustness of an unfolding network, produced by a forward-backward proximal interior point method. Moreover,  \cite{istagen,listagen,unfoldrnn} present generalization\footnote{Intuitively speaking, generalization pertains to the ability of a neural network to perform well on unseen data} error bounds, in terms of the Rademacher complexity of the hypothesis class consisting of the functions that a learnable ISTA network can implement. The concept of generalization error bounds is widely known in the field of statistical learning theory and studied by different -- albeit connected -- complexity terms, such as Rademacher complexity \cite{resrad}, Vapnik-Chervonenkis (VC) dimension \cite{vapnik}, stability \cite{stable} and robustness \cite{robust}.
\subsection{Motivation}
Our work is inspired by \cite{istagen}, \cite{admm-net}, \cite{admmdad}. Both \cite{admm-net} and \cite{admmdad} interpret the iterations of ADMM as layers of a neural network, which learns a decoder reconstructing the signals of interest from their (noisy) measurements. In \cite{admm-net}, the decoder is jointly learned with a square sparsifying dictionary -- initialized as a discrete cosine transform -- and the thresholds employed in the original ADMM. Then, the proposed framework is tested on real-world MRI images. The ADMM-based decoder of \cite{admmdad} is jointly learned with a redundant sparsifying analysis operator -- thus employing analysis sparsity in CS \cite{kr} -- and is tested on real-world image and speech datasets. The unfolding network designed in \cite{istagen} jointly learns an ISTA-based decoder and an orthogonal\footnote{the orthogonality constraint is reached by adding a corresponding regularization term} sparsifying dictionary -- hence imposing a synthesis sparsity model in CS \cite{analvssyn}. Moreover, the authors provide a generalization error bound for the hypothesis class consisting of the functions that their ISTA-net learns. In the end, they evaluate their theoretical results, by testing their proposed decoder on synthetic signals and the MNIST dataset \cite{mnist}.\\
In a similar spirit, we derive a decoder for analysis-sparsity-based CS (from now on called \emph{analysis CS}) by interpreting the iterations of the analysis-$l_1$ algorithm of \cite{tfocs:6} as layers of a DNN, which we call Decoding Network (\emph{DECONET}). DECONET jointly learns a redundant sparsifying analysis operator, combining that way a dictionary learning technique. We prefer to employ analysis sparsity instead of synthesis sparsity in CS, due to some advantages the former has compared to the latter. For example, analysis sparsity provides flexibility in modelling sparse signals, by leveraging the redundancy of the involved analysis operators (we refer to Section 2.2 for a detailed comparison between the two models). Among other optimization methods handling analysis sparsity, our choice of \cite{tfocs:6} is attributed to its optimal performance. From the mathematical point of view, we study the generalization ability of DECONET. To that end, we estimate the generalization error achieved by DECONET, in terms of the empirical Rademacher complexity of the hypothesis class consisting of all the functions/decoders DECONET can implement. To the best of our knowledge, we are the first to study the generalization ability of an unfolding network, that jointly learns a CS decoder and a redundant sparsifier. In the end, we evaluate our theoretical results by testing our proposed framework on two real-world image datasets.

\subsection{Key contributions}
Our key contributions are listed below.
\begin{enumerate}
    \item We build a new deep unfolding network dubbed DECONET, which jointly learns a) a decoder that solves the analysis CS problem b) a redundant analysis operator $W\in\mathbb{R}^{N\times n}$ ($N>n$) -- shared across the layers of DECONET -- for sparsification.
    \item We introduce the hypothesis class -- parameterized by $W$ -- of all the decoders DECONET can realize and restrict $W$ to be bounded in this class. On one hand, the boundedness of the learnable $W$ imposes a realistic structural constraint for the operator itself, that facilitates the estimation of the generalization error. On the other hand, the weight sharing assumption for $W$ leads us to a recurrent neural network with a moderate number of weights.
    \item We estimate the empirical Rademacher complexity of the aforementioned hypothesis class and use this estimate to derive meaningful generalization error bounds for DECONET. Our results showcase that the redundancy of $W$ and the number of layers $L$ affect the generalization ability of DECONET; roughly speaking, the generalization error scales like $\sqrt{NL}$.
    \item We confirm the validity of our theoretical guarantees by testing DECONET on real-world image datasets. Furthermore, we compare DECONET to an ISTA-net baseline \cite{istagen}; the latter jointly learns a decoder and an orthogonal sparsifier. Our experiments demonstrate that a) the generalization error of DECONET scales correctly with our theoretical findings b) DECONET outperforms the baseline in terms of the generalization error. This behaviour confirms improved performance when learning a redundant sparsifying transform instead of an orthogonal one.
\end{enumerate}

\subsection{Organization of the paper}
The rest of the paper is outlined as follows. In Section 2, we briefly introduce CS, compare synthesis to analysis sparsity model and present example algorithms that solve CS under each sparsity model. In Section 3, we interpret a state-of-the-art iterative algorithm for analysis CS as a neural network coined DECONET, formulate its associated hypothesis class and define the empirical Rademacher complexity of the latter. Section 4 is dedicated to boundedness results, which are then used in Section 5 to deliver our main Theorems; the latter provide meaningful upper bounds on the generalization error of DECONET. Section 6 is devoted to numerical experiments, where we evaluate our framework under multiple settings. Finally, we conclude in Section 7, wrapping up and giving some potential future directions.
\subsection{Notation}
We denote the cardinality of a set $S$ by $|S|$. For $n\in\mathbb{N}$, we write $[n]$ for the set of indices $\{1,2,\dots,n\}$. We denote the set of real, positive numbers by $\mathbb{R}_{+}$. For a sequence $a_n$ that is upper bounded by $M>0$, we write $\{a_n\}\leq M$; similarly if it is lower bounded. The support of a vector $x\in\mathbb{R}^n$ is the index set of its nonzero entries and is denoted by $\mathrm{supp}(x)$, i.e. $\mathrm{supp}(x)=\{i\in[n]:\,x_i\neq0\}$. For a matrix $A\in\mathbb{R}^{n\times n}$, we write $\|A\|_{2\rightarrow2}$ for its operator/spectral norm and $\|A\|_F$ for its Frobenius norm. For a family of vectors $(\phi_i)_{i=1}^N$ in $\mathbb{R}^n$, its associated analysis operator is given by $\Phi f:=\{\langle f,\phi_i\rangle\}_{i=1}^N$, where $f\in\mathbb{R}^n$. Its synthesis operator is simply the adjoint $\Phi^T$. For matrices $A_1,A_2\in\mathbb{R}^{N\times N}$, we denote by $[A_1;A_2]\in\mathbb{R}^{2N\times N}$ their concatenation with respect to the first dimension, while we denote by $[A_1\,|\,A_2]\in\mathbb{R}^{N\times 2N}$ their concatenation with respect to the second dimension. We write $O_{N\times N}$ for a real-valued $N\times N$ matrix filled with zeros. We denote by $\mathrm{diag}(\alpha)$ the square diagonal matrix having $\alpha\in\mathbb{R}$ in its main diagonal and zero elsewhere. For $x\in\mathbb{R},\,\tau>0$, the soft thresholding operator $\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}_\tau:\mathbb{R}\mapsto\mathbb{R}$ is defined as
    \begin{align}
    \mathcal{S}}\newcommand{\sgn}{\mathrm{sign}_\tau(x)=\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}(x,\tau)=\left\{\begin{array}{cc}
    \sgn(x)(|x|-\tau),&\ |x|\geq\tau\\
        0,&\ \text{otherwise}
        \end{array}\right.,
    \end{align} or in closed form
    \begin{equation}
    \label{sto}
    \mathcal{S}}\newcommand{\sgn}{\mathrm{sign}(x,\tau)=\sgn(x)\max(0,|x|-\tau).
    \end{equation} For $x\in\mathbb{R}^n$, the soft thresholding operator acts componentwise, i.e.\linebreak $(\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}_\tau(x))_i=\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}_\tau(x_i)$. For $y\in\mathbb{R}^n,\,\tau>0$, the mapping
    \begin{equation}\label{proxmap}
        P_G(\tau;y)=\mathrm{argmin}_{x\in\mathbb{R}^n}\left\{\tau G(x)+\frac{1}{2}\|x-y\|_2^2\right\},
    \end{equation}
    is called the \textit{proximal mapping associated to the convex function G}. In fact, for $G=\|\cdot\|_1$, \eqref{proxmap} coincides with \eqref{sto}. For $x\in\mathbb{R},\,\tau>0$, the truncation operator $\mathcal{T}_\tau:\mathbb{R}\mapsto\mathbb{R}$ is defined as
    \begin{align}
    \label{trunc}
        \mathcal{T}_\tau(x)=\mathcal{T}(x,\tau)=\sgn(x)\min\{|x|,\tau\}=\left\{\begin{array}{cc}
        \tau\sgn(x),&\ |x|\geq\tau\\
        x,&\ \text{otherwise}
        \end{array}\right. .
    \end{align}
    For $x\in\mathbb{R}^n$, the truncation operator acts componentwise and is 1-Lipschitz. The epigraph of the $l_2$-norm is the set $\mathcal{L}_2^n=\{(x,t)\in\mathbb{R}^{n+1}:\,\|x\|_2\leq t\}$. For two functions $f,g:\mathbb{R}^n\mapsto\mathbb{R}^n$, we write their composition as $f\circ g:\mathbb{R}^n\mapsto\mathbb{R}^n$ and if there exists some constant $C>0$ such that $f(x)\leq Cg(x)$, then we write $f(x)\lesssim g(x)$. The covering number $\mathcal{N}(T,\|\cdot\|,t)$ of a space $T$, equipped with a norm $\|\cdot\|$, at level $t>0$, is defined as the smallest number of balls of radius $t$ with respect to $\|\cdot\|$, required to cover $T$. For the ball of radius $t>0$ in $\mathbb{R}^n$, we write $B_{\|\cdot\|_2}^n(t)$. For $\Lambda>0$, the set of all matrices $W\in\mathbb{R}^{N\times n}$ with bounded -- by $\Lambda$ -- spectral norm is defined as $\mathcal{B}_\Lambda=\{W\in\mathbb{R}^{N\times n}\,\vert\,\exists\Lambda>0\ \text{such that}\ \|W\|_{2\rightarrow2}\leq\Lambda\}$.
\section{Model-based Compressed Sensing}
\subsection{A brief review of CS}
As already mentioned in Section 1, the main idea of CS is to reconstruct a vector $x\in\mathbb{R}^n$ from $y=Ax+e\in\mathbb{R}^m$, $m<n$, where $A$ is the so-called \emph{measurement matrix} \cite{rf} and $e\in\mathbb{R}^m$, with $\|e\|\leq\varepsilon$, corresponds to noise\footnote{for CS, we usually consider zero-mean Gaussian noise}, with $\varepsilon>0$ being an estimate on the noise level. In terms of classical linear algebra, the system of equations defined by $y=Ax+e$ is underdetermined, so there may be infinitely many solutions (provided of course that there exists at least one). In order to ensure exact/approximate reconstruction of $x$, we rely on two facts. First, $A$ must meet some conditions, for example the restricted isometry property or the null space property \cite{rf}. Second, we impose a \emph{sparse data model} to $x$ \cite{elad}. According to the latter, we consider $x$ to be $k$-sparse, that is, $\|x\|_0=\vert\mathrm{supp}(x)\vert\leq k$. Overall, we are called to solve the \emph{$l_0$-minimization} problem
\begin{equation}
    \min_{x\in\mathbb{R}^n}\|x\|_0\quad\text{subject to}\quad\|y-Ax\|_2\leq\varepsilon.
\end{equation}
However, the latter is NP-hard. A well-studied tractable alternative is the \emph{$l_1$-minimization} approach
\begin{equation}\label{l1min}
    \min_{x\in\mathbb{R}^n}\|x\|_1\quad\text{subject to}\quad\|y-Ax\|_2\leq\varepsilon.
\end{equation}
\subsection{Synthesis vs analysis sparsity model in CS}
Unfortunately, the sparsity assumption is rarely satisfied for real-world signals. However, signals are considered to be sparse when \emph{synthesized} by a few column vectors taken from a large matrix (\emph{dictionary}) $D\in\mathbb{R}^{p\times n}$ ($p\leq n$). In other words, $x=Dz$, where the coefficient vector $z\in\mathbb{R}^p$ is $k$-sparse and called the \emph{sparse representation} of $x$. The aforementioned model for $x$ is the \emph{synthesis sparsity model} and it is by now very well studied \cite{fistashear,rf,adcock,eeg}. Moreover, $D$ is usually an orthogonal transform, e.g. a wavelet or DCT matrix $D\in\mathbb{R}^{n\times n}$. Under the synthesis sparsity model with $D\in\mathbb{R}^{n\times n}$, \eqref{l1min} is transformed into
\begin{equation}
    \min_{z\in\mathbb{R}^n}\|z\|_1\quad\text{subject to}\quad\|y-ADz\|_2\leq\varepsilon
\end{equation}
and the latter is equivalent to
\begin{equation}\label{l1syn}
    \min_{z\in\mathbb{R}^n}\frac{1}{2}\|y-ADz\|_2^2+\lambda\|z\|_1,
\end{equation}
where $\lambda>0$ is a regularization parameter.\\
Nevertheless, the synthesis sparsity model has a ``twin'' named \emph{analysis sparsity model} \cite{genzel,kr,cosparse} (also known as \emph{co-sparse model} \cite{cosparse,kaba}). In this case, we assume there exists an \emph{analysis operator} $W\in\mathbb{R}^{N\times n}$ ($N\geq n$) so that the \emph{analysis representation} $s=Wx$ of $x$ is $k$-sparse. The associated optimization problem for CS is the \emph{analysis $l_1$-minimization} problem
\begin{equation}\label{l1anal}
        \min_{x\in\mathbb{R}^n}\|W x\|_1\quad\text{subject to}\quad \|Ax-y\|_2\leq\varepsilon
\end{equation}
and the latter is equivalent to
\begin{equation}\label{genlasso}
    \min_{x\in\mathbb{R}^n}\frac{1}{2}\|y-Ax\|^2_2+\lambda\|Wx\|_1.
\end{equation}
From now on, whenever we speak about the redundancy of such an analysis operator, we mean the number of its rows $N$.\\
Analysis sparsity has gained the attention of research community, due to some benefits it has compared to its synthesis counterpart. As we have already mentioned in Section 1.2, analysis sparsity model is more flexible in representing sparse signals. Moreover, it is computationally more appealing to solve the optimization algorithm of analysis CS, since the actual optimization takes place in the ambient space \cite{star} and the algorithm may need less measurements for perfect reconstruction, if one uses a redundant transform instead of an orthogonal one \cite{genzel}. As one may notice, the two models coincide precisely when $D,W$ are non-singular matrices, i.e. $D^{-1}=W$ (analogously $z=s$).
\subsection{Optimization-based algorithms for CS}
Most optimization-based algorithms for CS consist of an iterative scheme that incorporates a proximal mapping and synthesis sparsity. After a number of iterations and under certain conditions, the algorithm converges to a minimizer $\hat{x}$ of \eqref{l1syn}. For example, ISTA uses the proximal mapping\footnote{we remind at this point that for the convex penalty function $\|\cdot\|_1$, the proximal mapping with respect to $\|\cdot\|_1$ coincides to the soft thresholding operator \eqref{sto}} \eqref{proxmap} to yield the following iterative scheme
\begin{align}
    z_{k+1}=&\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}_{\tau\lambda}(z_k+\tau(AD)^T(y-ADz))\notag\\
    z_{k+1}=&\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}_{\tau\lambda}((I-D^TA^TAD)z+\tau(AD)^Tz),
\end{align}
for $k=0,1,\dots$, where $z_0=0$ and $\tau>0$ is a parameter of the algorithm. If $\tau\|AD\|_{2\rightarrow2}^2\leq1$ \cite{daubechies}, $z_k$ converges to a minimizer $\hat{z}$ of \eqref{l1syn}, so that the reconstructed $\hat{x}$ is simply given by $\hat{x}=D\hat{z}$\\
The drawback when using such a thresholding algorithm, is that it cannot handle the analysis sparsity model. In other words, when we have a penalty function of the form $\|W(\cdot)\|_1$, for $W\in\mathbb{R}^{N\times n}$ being a redundant transform ($N>n$), the proximal mapping associated to $\|W(\cdot)\|_1$ does not have a closed-form type. Hence, if we want to perform analysis CS, we have to come up with another iterative method. To tackle this issue, we choose the state-of-the-art $l_1$-analysis algorithm described in \cite{tfocs:6}, which employs a conic formulation methodology. For reasons of convenience, we will refer to the aforementioned algorithm as \emph{analysis conic form} (ACF) from now on. We will briefly describe the steps leading to the derivation of ACF, as these are stated in \cite{tfocs:6}. First, the authors of \cite{tfocs:6} determine an equivalent to \eqref{l1anal} smoothed conic formulation, i.e. 
\begin{equation}\label{l1conic}
    \min_{x\in\mathbb{R}^n}\|Wx\|_1+\frac{\mu}{2}\|x-x_0\|_2^2\quad\text{subject to}\quad(y-Ax,\varepsilon)\in\mathcal{L}_2^m,
\end{equation}
where $\mu\in\mathbb{R}_{+}$ is an adequate smoothing parameter, $x_0\in\mathbb{R}^n$ is an initial guess on $x$ and $\mathcal{L}_2^m$ is the epigraph of the $l_2$ norm. Second, they determine the dual of \eqref{l1conic} to be
\begin{equation}\label{dual}
    \begin{split}
        \text{maximize}\quad&\langle y,z^2\rangle-\varepsilon\|z^2\|_2\\
\text{subject to}\quad &A^Tz^2-W^Tz^1=0\\
&\|z^1\|_{\infty}\leq1,
    \end{split}
\end{equation}
where $z^1\in\mathbb{R}^N$, $z^2\in\mathbb{R}^m$ are dual variables. After presenting a collection of arguments and computations, they end up with Algorithm \ref{l1analysis} -- being a variant of an optimal first-order method -- stated below, with a step size multiplier $0<\{\theta_k\}\leq1$.

\begin{algorithm}[H]
\label{l1analysis}
\SetAlgoLined
\SetKwData{Left}{left}\SetKwData{This}{this}\SetKwData{Up}{up}\SetKwFunction{Union}{Union}\SetKwFunction{FindCompress}{FindCompress}\SetKwInOut{Input}{Input}\SetKwInOut{Output}{Output}
\Input{$x_0\in\mathbb{R}^n,\,z_0^1\in\mathbb{R}^N,\,z_0^2\in\mathbb{R}^m$, $\mu\in\mathbb{R}_+$, step sizes $\{t_k^1\},\{\,t_k^2\}$}
\Output{solution $\hat{x}_{\mu}$ of \eqref{l1conic}}
$\theta_0\leftarrow1,\,u_0^1=z_0^1,\,u_0^2=z_0^2$\;
\For{$\mathrm{iterations} \ k=0,1,\dots$}{
$x_k\leftarrow x_0+\mu^{-1}((1-\theta_k)W^Tu_k^1+\theta_kW^Tz_k^1-(1-\theta_k)A^Tu_k^2-\theta_kA^Tz_k^2)$\;
$z_{k+1}^1\leftarrow\mathcal{T}((1-\theta_k)u_k^1+\theta_kz_k^1-\theta_k^{-1}t_k^1Wx_k,\theta_k^{-1}t_k^1)$\;
$z_{k+1}^2\leftarrow\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}((1-\theta_k)u_k^2+\theta_kz_k^2-\theta_k^{-1}t_k^2(y-Ax_k),\theta_k^{-1}t_k^2\varepsilon)$\;
$u_{k+1}^1\leftarrow(1-\theta_k)u_k^1+\theta_kz_{k+1}^1$\;
$u_{k+1}^2\leftarrow(1-\theta_k)u_k^2+\theta_kz_{k+1}^2$\;
$\theta_{k+1}\leftarrow2/(1+(1+4/(\theta_k)^2)^{1/2})$\;
}
\caption{ACF}
\end{algorithm}

The dual function $g_\mu$ corresponding to \eqref{dual} has a Lipschitz continuous gradient, hence ACF converges \cite{tfocs:6} to a solution $\hat{x}_\mu$ of \eqref{l1conic}, for which we have $\hat{x}_{\mu}\overset{\mu\rightarrow0}{\longrightarrow}\hat{x}$, where $\hat{x}$ is an optimal solution of \eqref{l1analysis}. The authors of \cite{tfocs:6} clarify that when they speak about the optimal solution $\hat{x}$, they refer to this uniquely determined value. Additionally, they argue that there are situations where $\hat{x}$ and $\hat{x}_\mu$ coincide. Henceforward, we stick to their formulation and simply speak about the solution $\hat{x}$.\\
We will see in the next Section how ACF may be interpreted as a neural network with $L$ layers/iterations.
\section{A deep unfolding network for Compressed Sensing}
\subsection{Neural network formulation of the iterative algorithm}
We consider a typical scenario for the ACF\footnote{in terms of the associated MATLAB software that implements variants of ACF}, where $z_0^1=u_0^1=0$, $z_0^2=u_0^2=0$, $t^1_0=t^2_0=\theta_0=1$, $0<\{t_k^1\},\{\,t_k^2\},\{\theta_k\}\leq1$, $x_0=A^Ty$. We substitute first $x-$update into $z^1-$ and $z^2-$updates and second $z^1-$ and $z^2-$ into $u^1-$ and $u^2-$updates, respectively, concatenate $z_k^1,z_{k}^2,u_{k}^1,u_{k}^2$ in one vector $v_{k}$, i.e.
\begin{equation}
    v_k=\begin{pmatrix} z_{k}^1\\z_{k}^2\\u_{k}^1\\u_{k}^2\end{pmatrix}\in\mathbb{R}^{(2N+2m)\times1}\quad\text{for}\quad k\geq0,
\end{equation}
with $v_0=0$, and do the calculations, so that
\begin{align}
\label{finalv1}
v_k&=D_{k-1}v_{k-1}+\Theta_{k-1}\begin{pmatrix*}[l]
\mathcal{T}(G_{k-1}^1v_{k-1}-b_{k-1}^1,\theta_{k-1}^{-1}t_{k-1}^1)\\
\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}(G_{k-1}^2v_{k-1}-b_{k-1}^2,\theta_{k-1}^{-1}t_{k-1}^2\varepsilon)\\
\mathcal{T}(G_{k-1}^1v_{k-1}-b_{k-1}^1,\theta_{k-1}^{-1}t_{k-1}^1)\\
\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}(G_{k-1}^2v_{k-1}-b_{k-1}^2,\theta_{k-1}^{-1}t_{k-1}^2\varepsilon)\\
\end{pmatrix*},
\end{align}
where
\begin{align}
\label{dk}
    D_k=&\mathrm{diag}(\underbrace{1-\theta_0,\dots,1-\theta_0}_{N+m\, \text{times}},\underbrace{1-\theta_k,\dots,1-\theta_k}_{N+m\, \text{times}})\in\mathbb{R}^{(2N+2m)\times(2N+2m)}\\
    \label{thetak}
    \Theta_k=&\mathrm{diag}(\underbrace{\theta_0,\dots,\theta_0}_{N+m\, \text{times}},\underbrace{\theta_k,\dots,\theta_k}_{N+m\, \text{times}})\in\mathbb{R}^{(2N+2m)\times(2N+2m)}\\
G_k^1=&(\theta_k(I-\theta_k^{-1}t_k^1\mu^{-1}WW^T)\mid t_k^1\mu^{-1}WA^T\mid(1-\theta_k)\notag\\
    \label{g1}
    &\cdot(I-\theta_k^{-1}t_k^1\mu^{-1}WW^T)\mid(1-\theta_k)\theta_k^{-1}t_k^1\mu^{-1}WA^T)\in\mathbb{R}^{N\times(2N+2m)}\\
    G_k^2=&(t_k^2\mu^{-1}AW^T\mid\theta_k(I-\theta_k^{-1}t_k^2\mu^{-1}AA^T)\mid(1-\theta_k)\theta_k^{-1}t_k^2\mu^{-1}AW^T\notag\\
    \label{g2}
    &\mid(1-\theta_k)(I-\theta_k^{-1}t_k^2\mu^{-1}AA^T))\in\mathbb{R}^{m\times(2N+2m)}\\
    \label{bias1}
    b_k^1=&\theta_k^{-1}t_k^1Wx_0\in\mathbb{R}^N\\
    \label{bias2}
    b_k^2=&\theta_k^{-1}t_k^2(y-Ax_0)\in\mathbb{R}^m.
\end{align}


We observe that \eqref{finalv1} can be interpreted as a layer of a neural network, with weights $G^1, G^2$, biases $b^1,b^2$ and activation functions $\mathcal{T},\mathcal{S}}\newcommand{\sgn}{\mathrm{sign}$. Nevertheless, this interpretation of ACF as a DNN does not account for any trainable parameters. We cope with this issue by considering $W$ to be a) unknown b) bounded with respect to the operator norm, i.e. $W\in\mathcal{B}_\Lambda$, for some $\Lambda>0$ c) learned from a training sequence $\mathcal{S}=\{(x_i,y_i)\}_{i=1}^s$ with i.i.d. samples drawn from an unknown distribution\footnote{formally speaking, this is a distribution over the $x_i$ and then $y_i=Ax_i+e$, with fixed $A,e$} $\mathcal{D}^s$. Hence, the trainable parameters are the entries of $W$.\\
Now, based on \eqref{finalv1}, we formulate ACF as a neural network with $L$ layers/iterations, defined as
\begin{align}
   
   
   
    \label{layer1}
    f_1(y)&=\sigma(y)\\
    \label{layerk}
   
    f_k(v)&=D_{k-1}v+\Theta_{k-1}\sigma(v),\quad k=2,\dots,L,
\end{align}
where
\begin{equation}
    \begin{split}
        \sigma(y)^T=(&\mathcal{T}(-t_0^1Wx_0,t_0^1),\mathcal{S}(t_0^2(y-Ax_0),t_0^2\varepsilon),\\
        &\mathcal{T}(-t_0^1Wx_0,t_0^1),\mathcal{S}(t_0^2(y-Ax_0),t_0^2\varepsilon))^T,
    \end{split}
\end{equation}
\begin{equation}\label{vsigma}
    \begin{split}
    \sigma(v)^T=(&\mathcal{T}(G_k^1v-b_k^1),\mathcal{S}(G_k^2v-b_k^2),\\ &\mathcal{T}(G^1_kv-b^1_k),\mathcal{S}(G_k^2v-b_k^2))^T,\quad k=2,\dots,L.
    \end{split}
\end{equation}
We denote the concatenation of $L$ such layers (all having the same $W$) as
\begin{equation}
\label{composed}
    f_{W}^L(y)=f_L\circ f_{L-1}\circ\dots\circ f_1(y).
\end{equation}
The latter constitutes the realization of a neural network with $L$ layers, that reconstructs $v$ from $y$. Thus, we call \eqref{composed} \emph{dual decoder}, since it is not the final form of the decoder we want to derive. Towards this end, in order to get the solution $\hat{x}$, we apply an affine map $\phi:\mathbb{R}^{(2N+2m)\times 1}\mapsto\mathbb{R}^{n\times1}$ after the last layer $L$, so that
\begin{equation}
    \hat{x}:=\phi(v)=\Phi v+x_0,
\end{equation}
where
\begin{equation}
    \begin{split}
    \label{phi}
        \Phi=&(\mu^{-1}\theta_LW^T\vert-\mu^{-1}\theta_LA^T\\
        &\vert\ \mu^{-1}(1-\theta_L)W^T\vert-\mu^{-1}(1-\theta_L)A^T)\in\mathbb{R}^{n\times (2N+2m)}.
    \end{split}
\end{equation}
Moreover, in order to clip the output $\phi(f_W^L(y))$ in case its norm falls out of a reasonable range, we add an extra function $\psi : \mathbb{R}^n \to \mathbb{R}^n$ after the application of the affine map $\phi$ and define it as
\begin{equation}\label{psi}
    \psi(x) = \left\{\begin{array}{cc}
    x,&\ \|x\|_2 \leq B_{\mathrm{out}}\\
    B_{\mathrm{out}}\frac{x}{\|x\|_2},&\ \text{otherwise}
    \end{array}\right.,
\end{equation} for some fixed constant $B_{\mathrm{out}}>0$.
Now, for a fixed number of layers $L$, the desired learned decoder is written as
\begin{equation}\label{decoder}
    \text{dec}_W^L(y)=\psi(\phi(f^L_W(y))).
\end{equation}
We call DECONET (DECOding NETwork) the neural network that implements such a decoder. Notice that the latter is parameterized by $W$, since $W$ is shared across the layers of DECONET.

\subsection{Defining the hypothesis class and the Rademacher complexity}
We introduce the hypothesis class
\begin{equation}\label{hypothesis}
    \mathcal{H}^L=\{ h:\mathbb{R}^m\mapsto\mathbb{R}^n:\, h(y)=\psi(\phi(f^L_{W}(y))),\,W\in\mathcal{B}_\Lambda\},
\end{equation}
parameterized by $W$ and consisting of all the functions/decoders DECONET can implement. Given \eqref{hypothesis} and the training set $\mathcal{S}$, DECONET yields a function $h_{\mathcal{S}}\in\mathcal{H}^L$ that aims at reconstructing $x$ from $y$. For a loss function $\ell:\mathcal{H}^L\times\mathbb{R}^n\times\mathbb{R}^m\mapsto\mathbb{R}_+$, the empirical loss of a hypothesis $h_\mathcal{S}$ is the reconstruction error on the training set, i.e.
\begin{equation}
    \mathcal{\hat{L}}_{train}(h_\mathcal{S})=\frac{1}{s}\sum_{i=1}^s\ell(h_\mathcal{S},x_i,y_i).
\end{equation}
In this paper, we choose as loss function $\ell$ the squared $l_2$-norm, which is considered a typical measure of reconstruction error in regression-style problems like CS. Thus, the empirical loss takes the form of the training \emph{mean-squared error} (MSE)
\begin{equation}\label{trainmse}
    \mathcal{\hat{L}}_{train}(h_{\mathcal{S}})=\frac{1}{s}\sum_{j=1}^s\|h_{\mathcal{S}}(y_j)-x_j)\|_2^2.
\end{equation}
The true loss is
\begin{equation}
    \mathcal{L}(h_{\mathcal{S}})=\mathbb{E}_{(x,y)\sim\mathcal{D}}(\|h_{\mathcal{S}}(y)-x\|_2^2).
\end{equation}
Now, the generalization error is given as the difference\footnote{some of the existing literature denotes the true loss as the generalization error, but the definition we give in \eqref{generror} is more convenient for our purposes} between the empirical and true loss
\begin{equation}\label{generror}
    GE(h_{\mathcal{S}})=|\mathcal{\hat{L}}_{train}(h_{\mathcal{S}})-\mathcal{L}(h_{\mathcal{S}})|.
\end{equation}

\begin{remark}
The interested reader may wonder why we do not choose a loss function $\ell(\cdot)$ like $\|\cdot\|_2$, which is easier to plug and play in mathematical computations, due to the fact that it has a Lipschitz constant equal to one. The reason is we prefer to be consistent with the forthcoming numerical experiments, where we train DECONET with respect to $\|\cdot\|_2^2$.
\end{remark}

A typical way to estimate \eqref{generror} consists in upper bounding it in terms of the \emph{Rademacher complexity}. The \emph{empirical Rademacher complexity} is defined as \begin{equation}
    \label{erc}
    \mathcal{R}_{\mathcal{S}}(\ell\circ\mathcal{H}^L)=\mathbb{E}\sup_{h\in\mathcal{H}^L}\frac{1}{s}\sum_{i=1}^s\epsilon_i\|h(y_i)-x_i\|_2^2,
\end{equation}
where $\epsilon$ is a Rademacher vector, that is, a vector with entries taking the values $\pm1$ with equal probability. Then, the Rademacher complexity is defined as
\begin{equation}\label{rademacher}
    \mathcal{R}_s(\ell\circ\mathcal{H}^L)=\mathbb{E}_{{\mathcal{S}}\sim\mathcal{D}^s}(\mathcal{R}_{\mathcal{S}}(\ell\circ\mathcal{H}^L)).
\end{equation}
In this paper, we solely work with \eqref{erc}. We rely on the following Theorem that estimates \eqref{generror} in terms of \eqref{erc}.
\begin{theorem}\label{radem}
Let $\mathcal{H}$ be a family of functions, $\mathcal{S}$ the training set drawn from $\mathcal{D}^s$, and $\ell$ a real-valued bounded loss function satisfying $|\ell(h,z)|\leq c$, for all $h\in\mathcal{H}, z\in Z$. Then, for $\delta\in(0,1)$, with probability at least $1-\delta$, we have for all $h\in\mathcal{H}$
\begin{equation}
    \mathcal{L}(h)\leq\mathcal{\hat{L}}(h)+2\mathcal{R}_{\mathcal{S}}(\ell\circ\mathcal{H})+4c\sqrt{\frac{2\log(4\delta)}{s}}.
\end{equation}
\end{theorem}
In order to use the previous Theorem, the loss function must be bounded. Towards this end, we make two reasonable (for the machine learning literature) assumptions regarding the training set $\mathcal{S}=\{(x_i,y_i)\}_{i=1}^s$. Let us suppose that  with overwhelming probability we have
\begin{equation}
    \label{ybound}
    \|y_i\|_2\leq\mathrm{B_{in}},
\end{equation}
for some constant $\mathrm{B_{in}}>0$, $i=1,\dots,s$. Moreover, we assume that for any $h\in\mathcal{H}^L$, with overwhelming probability over $y_i$ chosen from $\mathcal{D}$, the following holds
\begin{equation}\label{hbound}
    \|h(y_i)\|_2\leq B_{\text{out}},
\end{equation}
by definition of $\psi$, for some constant $B_{\text{out}}>0$, for all $i=1,\dots,s$. Hence, the loss function is bounded as $\|h(y_i)-x_i\|_2^2\leq (B_{\text{in}}+B_{\text{out}})^2$, for all $i=1,\dots,s$. Following the previous assumptions, it is easy to check that $\|\cdot\|^2_2$ is a Lipschitz function, with Lipschitz constant $\mathrm{Lip}_{\|\cdot\|^2_2}=2\mathrm{B_{in}}+2\mathrm{B_{out}}$.\\
The Lipschitzness of $\|\cdot\|_2^2$ allows us to remove the loss function in \eqref{erc} and study $\mathcal{R}_S(\mathcal{H})$ alone. To do so, we employ the so-called (vector-valued) contraction principle \cite{contraction}.
\begin{lemma}\label{contraction}
Let $\mathcal{H}$ be a set of function $h:\mathcal{X}\mapsto\mathbb{R}^n$, $f:\mathbb{R}^n\mapsto\mathbb{R}^n$ a $K$-Lipschitz function and $\mathcal{S}=\{x_i\}_{i=1}^s$. Then
\begin{equation}
    \mathbb{E}\sup_{h\in\mathcal{H}}\sum_{i=1}^s\epsilon_if\circ h(x_i)\leq\sqrt{2}K\mathbb{E}\sup_{h\in\mathcal{H}}\sum_{i=1}^s\sum_{k=1}^n\epsilon_{ik}h_k(x_i),
\end{equation}
where $(\epsilon_i),(\epsilon_{ik})$ are both Rademacher sequences.
\end{lemma}
We apply the previous Lemma to \eqref{erc}, yielding
\begin{align}
    \mathcal{R}_s(l\circ\mathcal{H}^L)\leq&\sqrt{2}\mathrm{Lip}_{\|\cdot\|^2_2}\mathcal{R}_s(\mathcal{H}^L)=\sqrt{2}\mathrm{Lip}_{\|\cdot\|^2_2}\mathbb{E}\sup_{h\in\mathcal{H}^L}\sum_{i=1}^s\sum_{k=1}^n\epsilon_{ik}h_k(x_i)\notag\\\label{contradem}
    =&\sqrt{2}(2\mathrm{B_{in}}+2\mathrm{B_{out}})\mathbb{E}\sup_{h\in\mathcal{H}^L}\sum_{i=1}^s\sum_{k=1}^n\epsilon_{ik}h_k(x_i).
\end{align}
We will come back to the estimation of \eqref{contradem} in Section 5, after presenting the adequate mathematical tools in Section 4. Towards that end, we loosely follow the mathematical strategy described in \cite{istagen}.
\section{Boundedness results}
We present a series of boundedness results that are essential for the estimation of the generalization error bounds.
\subsection{Bounding outputs}
We take into account the number of training samples and pass to matrix notation. Due to \eqref{ybound}, \eqref{hbound} and the Cauchy-Schwartz inequality, we get
\begin{align}
\label{boundedin}
\|Y\|_F&\leq\sqrt{s}\mathrm{B}_{\text{in}}\\
\label{boundedout}
\|h(Y)\|_F&=\|\psi(\phi(f_{W}^L(Y)))\|_F\leq\sqrt{s}\mathrm{B}_{\text{out}}.
\end{align}
We will make wide use of the following inequality, so we state it below as a Lemma.

\begin{lemma}\label{g12bound}
Let $k\geq0$. For any $W\in\mathcal{B}_\Lambda$, step sizes $\{t_k^1\},\{t_k^2\}>0$ with $t_0^1=t_0^2=1$, step size multiplier $0<\{\theta_k\}\leq1$ with $\theta_0=1$, and smoothing parameter $\mu>0$, the following holds for the matrices $G^1_k,\,G_k^2$ defined in \eqref{g1}, \eqref{g2}, respectively:
\begin{equation}
    \begin{split}\label{g1+2}
        2\|G_k^1\|_{2\rightarrow2}+2\|G_k^2\|_{2\rightarrow2}+1\leq\Gamma_k,
    \end{split}
\end{equation}
where
\begin{equation}\label{gammak}
    \Gamma_k=2\big[c_{1,k}\Lambda^2+c_{2,k}\|A\|_{2\rightarrow2}^2+2\|A\|_{2\rightarrow2}\Lambda(c_{1,k}+c_{2,k})\big]+1,
\end{equation}
with $\Gamma_k\leq6(\max\{\Lambda,\|A\|_{2\rightarrow2})^2(c_{1,k}+c_{2,k})\}+1$ and $c_{1,k}=\theta^{-1}_{k}\mu^{-1}t^1_{k}$, $c_{2,k}=\theta^{-1}_{k}\mu^{-1}t^2_{k}$, for all $k\geq0$. Moreover, if  $\{c_{1,k}\},\{c_{1,k}\}\leq1$, then 
\begin{equation}
    \Gamma_k\leq\gamma,
\end{equation}
where\footnote{thus, for any $\Lambda>0$ and $\|A\|_{2\rightarrow2}\geq0$, we always have $\gamma>0$.} $\gamma=2(\Lambda^2+\|A\|_{2\rightarrow2}^2+4\|A\|_{2\rightarrow2}\Lambda)+1>1$.
\end{lemma}
\begin{proof}
Based on \eqref{g1}, \eqref{g2}, we get
\begin{gather*}
   \|G^1_{k}\|_{2\rightarrow2}+\|G^2_{k}\|_{2\rightarrow2}\leq\big[\theta_{k}\|I-\theta^{-1}_{k}\mu^{-1}t^1_{k}WW^T\|_{2\rightarrow2}\\
   +\mu^{-1}t^1_{k}\|W\|_{2\rightarrow2}\|A\|_{2\rightarrow2}+(1-\theta_{k})\|I-\theta^{-1}_{k}\mu^{-1}t^1_{k}WW^T\|_{2\rightarrow2}\notag\\
   +(1-\theta_{k})\theta^{-1}_{k}\mu^{-1}t^1_{k}\|W\|_{2\rightarrow2}\|A\|_{2\rightarrow2}\big]\\
   +\big[\mu^{-1}t^2_{k}\|A\|_{2\rightarrow2}\|W\|_{2\rightarrow2}+\theta_{k}\|I-\theta^{-1}_{k}\mu^{-1}t^2_{k}AA^T\|_{2\rightarrow2}+(1-\theta_{k})\\
   \cdot\theta^{-1}_{k}\mu^{-1}t^2_{k}\|A\|_{2\rightarrow2}\|W\|_{2\rightarrow2}+(1-\theta_{k})\|I-\theta^{-1}_{k}\mu^{-1}t^2_{k}AA^T\|_{2\rightarrow2}\big]\notag\\
   =[\|I-\theta^{-1}_{k}\mu^{-1}t^1_{k}WW^T\|_{2\rightarrow2}+\|W\|_{2\rightarrow2}\|A\|_{2\rightarrow2}(\theta^{-1}_{k}\mu^{-1}t_{k}^1+\mu^{-1}t_{k}^1)]\\
   +[\|I-\theta^{-1}_{k}\mu^{-1}t^2_{k}AA^T\|_{2\rightarrow2}+\|W\|_{2\rightarrow2}\|A\|_{2\rightarrow2}(\theta^{-1}_{k}\mu^{-1}t^2_{k}+\mu^{-1}t^2_{k})].
\end{gather*}
Now, we define $c_{1,k}=\theta^{-1}_{k}\mu^{-1}t^1_{k}$, $c_{2,k}=\theta^{-1}_{k}\mu^{-1}t^2_{k}$, use the boundedness assumption of $W$ and the fact that $\|I-c_{1,k}WW^T\|_{2\rightarrow2}\leq\max\{{c_{1,k}\Lambda}-1,1\}\leq{c_{1,k}\Lambda}$ to get
\begin{gather}
    2\|G^1_{k}\|_{2\rightarrow2}+2\|G^2_{k}\|_{2\rightarrow2}+1\notag\\
    \leq2\big[c_{1,k}\Lambda^2+c_{2,k}\|A\|_{2\rightarrow2}^2   +\|A\|_{2\rightarrow2}\Lambda(c_{1,k}+\theta_k c_{1,k}+c_{2,k}+\theta_k c_{2,k})\big]+1\notag\\
    \label{g12proved}
    \overset{\theta_k\leq1}{\leq}2\big[c_{1,k}\Lambda^2+c_{2,k}\|A\|_{2\rightarrow2}^2
    +2\|A\|_{2\rightarrow2}\Lambda(c_{1,k}+c_{2,k})\big]+1.
\end{gather}

Now, we set $\Gamma_k=2\big[c_{1,k}\Lambda^2+c_{2,k}\|A\|_{2\rightarrow2}^2+2\|A\|_{2\rightarrow2}\Lambda(c_{1,k}+c_{2,k})\big]+1$. If we take the maximum between $\|A\|_{2\rightarrow2}^2$ and $\Lambda^2$, then
\begin{align}
    \Gamma_k&\leq2\big[3c_{1,k}(\max\{\Lambda,\|A\|_{2\rightarrow2}\})^2+3c_{2,k}(\max\{\Lambda,\|A\|_{2\rightarrow2}\})^2\big]+1\notag\\
    &=6(\max\{\Lambda,\|A\|_{2\rightarrow2})^2(c_{1,k}+c_{2,k})+1.
\end{align}
Combining the latter with the fact that $\Gamma_k\geq1$, for any $k\geq0$, gives us the restriction of $\Gamma_k$ in the desired interval. Moreover, if $\{c_{1,k}\},\{c_{2,k}\}\leq1$ for any $k\geq0$, then we have the following simplified upper bound for $\Gamma_k$, that does not depend on $k$:
\begin{gather}
    \Gamma_k
    \leq2(\Lambda^2+\|A\|_{2\rightarrow2}^2+4\|A\|_{2\rightarrow2}\Lambda)+1\notag.
\end{gather}
We set $\gamma:=2(\Lambda^2+\|A\|_{2\rightarrow2}^2+4\|A\|_{2\rightarrow2}\Lambda)+1$ (obviously $\gamma>1$) and the proof follows.
\end{proof}

\begin{remark}
As one may notice, we could simply set $\Gamma_k=2\big[c_{1,k}\Lambda^2+c_{2,k}\|A\|_{2\rightarrow2}^2+2\|A\|_{2\rightarrow2}\Lambda(c_{1,k}+c_{2,k})\big]$, so that $2\|G_k^1\|_{2\rightarrow2}+2\|G_k^2\|_{2\rightarrow2}\leq\Gamma_k$. However, as it will become apparent later on, the quantity $2\|G_k^1\|_{2\rightarrow2}+2\|G_k^2\|_{2\rightarrow2}+1$ is easier to handle for our mathematical purposes, so we preferred to estimate $2\|G_k^1\|_{2\rightarrow2}+2\|G_k^2\|_{2\rightarrow2}+1$ instead of $2\|G_k^1\|_{2\rightarrow2}+2\|G_k^2\|_{2\rightarrow2}$.
\end{remark}



Apart from the boundedness assumptions we have presented so far, we can upper-bound the output $f^k_W(Y)$ with respect to the Frobenius norm, after any number of layers $k$ and especially for $k<L$, so that $\phi$ and $\psi$ are not applied after the final layer $L$. The following Lemma will be needed later on, when we prove that $f_W^L(Y)$ is Lipschitz continuous with respect to $W$.
\begin{lemma}\label{fbound}
Let $k\in\mathbb{N}$. For any $W\in\mathcal{B}_\Lambda$, step sizes $\{t_k^1\}$, $\{t_k^2\}>0$ with $t_0^1=t_0^2=1$, $t_{-1}^1=t_{-1}^2=0$, step size multiplier $0<\{\theta_k\}\leq1$ with $\theta_0=\theta_{-1}=1$, and smoothing parameter $\mu>0$, the following holds for the output of the function $f_W^k$ defined in \eqref{layerk}:
\begin{align}
\|f^k_W(Y)\|_F\leq&2\mu\|Y\|_F\sum_{i=0}^{k-1}\Bigg(\bigg(\|A\|_{2\rightarrow2}(c_{1,i-1}\Lambda+c_{2,i-1}\|A\|_{2\rightarrow2})+ c_{2,i-1}\bigg)\prod_{j=i}^{k-1}\Gamma_j\Bigg)\notag\\
\label{finductive}
&+2\mu\|Y\|_F\bigg(\|A\|_{2\rightarrow2}(c_{1,k-1}\Lambda+c_{2,k-1}\|A\|_{2\rightarrow2})+ c_{2,k-1}\bigg),
\end{align}
where $\{\Gamma_k\}_{k\geq0}$, $\{c_{1,k}\}_{k\geq0}$, $\{c_{2,k}\}_{k\geq0}$ are defined as in Lemma \ref{g12bound} and $c_{1,-1}=c_{2,-1}=0$. In particular, if $\{c_{1,k}\},\{c_{2,k}\}\leq1$, then we have the simplified upper bound
\begin{align}
\label{finductive2}
    \|f^k_W(Y)\|_F\leq&2\mu\|Y\|_F(\|A\|_{2\rightarrow2}(\Lambda+\|A\|_{2\rightarrow2})+1)(\zeta_k+1),
\end{align}
where
\begin{equation}\label{zeta}
\zeta_k=\frac{\gamma^k-1}{\gamma-1}
\end{equation}
and $\gamma$ is defined as in Lemma \ref{g12bound}.
\end{lemma}

\begin{proof}
First, we notice that both $\mathcal{T}(\cdot)$ and $\mathcal{S}(\cdot)$ are 1-Lipschitz functions. Second, by definition of the matrices $D_k$ and $\Theta_k$ in \eqref{dk} and \eqref{thetak} respectively, we have $\|D_k\|_{2\rightarrow2}\leq1$ and $\|\Theta_k\|_{2\rightarrow2}=1$, for any $k\geq1$, since $0<\{\theta_k\}\leq1$ and $\theta_0=1$. Now we use the previous statements, along with \eqref{layer1} and \eqref{layerk}, to prove \eqref{finductive} via induction. For $k=1$, we have
\begin{align*}
    \|f_W^1(Y)\|_F&\leq2t_0^1\Lambda\|X_0\|_F+2t_0^2(\|Y\|_F+\|A\|_{2\rightarrow2}\|X_0\|_F)\\
    &=2\mu c_{1,0}\Lambda\|A\|_{2\rightarrow2}\|Y\|_F+2\mu c_{2,0}(\|Y\|_F+\|A\|^2_{2\rightarrow2}\|Y\|_F)\\
    &=2\mu\|Y\|_F\bigg(\|A\|_{2\rightarrow2}(c_{1,0}\Lambda+c_{2,0}\|A\|_{2\rightarrow2})+c_{2,0}\bigg).
\end{align*}
Suppose \eqref{finductive} holds for $k$. Then, for $k+1$:
\begin{gather*}
    \|f_W^{k+1}(Y)\|_F<\|f_W^k(Y)\|_F+2\|G_k^1f_W^k(Y)-B_k^1\|_F+2\|G_k^2f_W^k(Y)-B_k^2\|_F\\
    \leq\|f_W^k(Y)\|_F(2\|G_k^1\|_{2\rightarrow2}+2\|G_k^2\|_{2\rightarrow2}+1)+2\|B_k^1\|_F+2\|B_k^2\|_F\\
    \overset{Lemma\  \ref{g12bound}}{\leq}\Gamma_k\|f_W^k(Y)\|_F+2\mu\|X_0\|_F(c_{1,k}\Lambda+c_{2,k}\|A\|_{2\rightarrow2})+2\mu c_{2,k}\|Y\|_F\\
    \leq\Gamma_k2\mu\|Y\|_F\sum_{i=0}^{k-1}\Bigg(\bigg(\|A\|_{2\rightarrow2}(c_{1,i-1}\Lambda+c_{2,i-1}\|A\|_{2\rightarrow2})+ c_{2,i-1}\bigg)\prod_{j=i}^{k-1}\Gamma_j\Bigg)\notag\\
    +\Gamma_k2\mu\|Y\|_F\bigg(\|A\|_{2\rightarrow2}(c_{1,k-1}\Lambda+c_{2,k-1}\|A\|_{2\rightarrow2})+ c_{2,k-1}\bigg)\\
    +2\mu\|Y\|_F\bigg(\|A\|_{2\rightarrow2}(c_{1,k}\Lambda+c_{2,k}\|A\|_{2\rightarrow2})+ c_{2,k}\bigg)\\
    =2\mu\|Y\|_F\sum_{i=0}^{k}\Bigg(\bigg(\|A\|_{2\rightarrow2}(c_{1,i-1}\Lambda+c_{2,i-1}\|A\|_{2\rightarrow2})+ c_{2,i-1}\bigg)\prod_{j=i}^{k}\Gamma_j\Bigg)\notag\\
    +2\mu\|Y\|_F\bigg(\|A\|_{2\rightarrow2}(c_{1,k}\Lambda+c_{2,k}\|A\|_{2\rightarrow2})+ c_{2,k}\bigg)
\end{gather*}
Therefore, we proved that \eqref{finductive} holds for any $k\in\mathbb{N}$. Now, under the additional assumptions $\{c_{1,k}\}_{k\geq0},\{c_{2,k}\}_{k\geq0}\leq1$, we may apply Lemma \ref{g12bound} on \eqref{finductive}, yielding
\begin{align}
    \|f_W^k(Y)\|_F&\leq2\mu\|Y\|_F\bigg(\|A\|_{2\rightarrow2}(\Lambda+\|A\|_{2\rightarrow2})+1\bigg)\Bigg(\sum_{i=0}^{k-1}\prod_{j=i}^{k-1}\gamma+1\Bigg)\notag\\
    &=2\mu\|Y\|_F\bigg(\|A\|_{2\rightarrow2}(\Lambda+\|A\|_{2\rightarrow2})+1\bigg)\Bigg(\sum_{i=1}^{k}\gamma^i+1\Bigg)\notag\\
    &=2\mu\|Y\|_F\bigg(\|A\|_{2\rightarrow2}(\Lambda+\|A\|_{2\rightarrow2})+1\bigg)(\zeta_k+1)\notag,
\end{align}
where
\begin{equation}
    \zeta_k=\frac{\gamma^k-1}{\gamma-1}
\end{equation}
and $\gamma$ is defined as in Lemma \ref{g12bound}.
\end{proof}

\subsection{Lipschitzness results}
We first prove that the dual decoder defined in \eqref{composed} is Lipschitz continuous with respect to any $W\in\mathcal{B}_\Lambda$.
\begin{theorem}\label{perturbation}
Let $f^L_W$ defined as in \eqref{composed}, $L\geq2$, dictionary $W\in\mathcal{B}_\Lambda$, step sizes $\{t_k^1\}_{k\geq-1},\{t_k^2\}_{k\geq-1}>0$ with $t_0^1=t_0^2=1$, $t_{-1}^1=t_{-1}^2=0$, step size multiplier $0<\{\theta_k\}_{k\geq-1}\leq1$ with $\theta_0=\theta_{-1}=1$, and smoothing parameter $\mu>0$. Then, for any $W_1,W_2\in\mathcal{B}_\Lambda$, we have
\begin{equation}\label{fpert}
    \|f^L_{W_1}(Y)-f^L_{W_2}(Y)\|_F\leq K_L\|W_1-W_2\|_{2\rightarrow2},
\end{equation} where
\begin{equation}
    \label{kl}
    \begin{split}
    K_L=&2\mu\|Y\|_F\Bigg[\mu^{-1}\|A\|_{2\rightarrow2}+\sum_{k=2}^L\Bigg(\left(\max_{0\leq l\leq L-1}\Gamma_l\right)^{L-k}2\Bigg[\sum_{i=0}^{k-2}\Bigg(\bigg(\|A\|_{2\rightarrow2}\\
    &\cdot(c_{1,i-1}\Lambda+c_{2,i-1}\|A\|_{2\rightarrow2})+c_{2,i-1}\bigg)\prod_{j=i}^{k-2}\Gamma_j\Bigg)\\
    &+\bigg(\|A\|_{2\rightarrow2}(c_{1,k-2}\Lambda+c_{2,k-2}\|A\|_{2\rightarrow2})+c_{2,k-2}\bigg)\Bigg]\\
    &\cdot(2\Lambda c_{1,k-1}+\|A\|_{2\rightarrow2}(c_{1,k-1}+c_{2,k-1}))+c_{1,k-1}\|A\|_{2\rightarrow2}\Bigg)\Bigg],
    \end{split}
\end{equation}
with $\{\Gamma_k\}_{k\geq0}$, $\{c_{1,k}\}_{k\geq0}$, $\{c_{2,k}\}_{k\geq0}$ defined as in Lemma \ref{g12bound} and $c_{1,-1}=c_{2,-1}=0$. Moreover, if  $\{c_{1,k}\},\{c_{2,k}\}\leq1$, then we have the simplified upper bound
\begin{align}
\label{klsimple}
    K_L<&2\mu\|Y\|_F\bigg[\|A\|_{2\rightarrow2}(\mu^{-1}+(L-1))\notag\\
    &+4\|A\|_{2\rightarrow2}((\Lambda+\|A\|_{2\rightarrow2})+1)(\Lambda+\|A\|_{2\rightarrow2})\kappa_L\bigg],
\end{align}
with
\begin{align}
\label{kappal}
    \kappa_L=\frac{\gamma^L((L-2)+\gamma^2(\gamma-2))-\gamma^2(\gamma-2)}{\gamma-1}
\end{align}
and $\gamma$ as in Lemma \ref{g12bound}.
\end{theorem}
\begin{proof}
First, we set $f^0_{W_1}(Y)=f^0_{W_2}(Y)=Y$ for a uniform treatment of all layers. Then, we write $\{G_{1,k}^i\}_{i=1,2},\,\{G_{2,k}^i\}_{i=1,2}$ (similarly for $\{B^i_{1,k}\}_{i=1,2}$, $\{B^i_{2,k}\}_{i=1,2}$) to denote the dependency on $W_1,W_2$, respectively. By the 1-Lipschitzness of $\mathcal{T}(\cdot),\mathcal{S}(\cdot)$, the estimates $\|D_k\|_{2\rightarrow2}\leq1$ and $\|\Theta_k\|_{2\rightarrow2}=1$ that hold for any $k\geq0$, and the introduction of mixed terms, we get
\begin{gather*}
    \|f^k_{W_1}(Y)-f^k_{W_2}(Y)\|_F\\
    \leq\|D_{k-1}f^{k-1}_{W_1}(Y)+\Theta_{k-1}\sigma(f_{W_1}^{k-1})-D_{k-1}f^{k-1}_{W_2}(Y)-\Theta_{k-1}\sigma(f_{W_2}^{k-1})\|_F\\
    \leq\|D_{k-1}\|_{2\rightarrow2}\|f^{k-1}_{W_1}(Y)-f^{k-1}_{W_2}(Y)\|_F+\|\Theta_{k-1}\|_{2\rightarrow2}\|\sigma(f_{W_1}^{k-1})-\sigma(f_{W_2}^{k-1})\|_F\\
    \leq\|f^{k-1}_{W_1}(Y)-f^{k-1}_{W_2}(Y)\|_F+2\|G_{1,k-1}^1f_{W_1}^{k-1}(Y)-B_{1,k-1}^1-G_{2,k-1}^1f_{W_2}^{k-1}(Y)\\
    +B_{2,k-1}^1\|_F+2\|G_{1,k-1}^2f_{W_1}^{k-1}(Y)-B_{1,k-1}^2-G_{2,k-1}^2f_{W_2}^{k-1}(Y)+B_{2,k-1}^2\|_F\\
    =\|f^{k-1}_{W_1}(Y)-f^{k-1}_{W_2}(Y)\|_F+2\|G_{1,k-1}^1f_{W_1}^{k-1}(Y)-B_{1,k-1}^1+G_{1,k-1}^1f_{W_2}^{k-1}(Y)\\
    -G_{1,k-1}^1f_{W_2}^{k-1}(Y)-G_{2,k-1}^1f_{W_2}^{k-1}(Y)+B_{2,k-1}^1\|_F+2\|G_{1,k-1}^2f_{W_1}^{k-1}(Y)-B_{1,k-1}^2\\
    +G_{1,k-1}^2f_{W_2}^{k-1}(Y)-G_{1,k-1}^2f_{W_2}^{k-1}(Y)-G_{2,k-1}^2f_{W_2}^{k-1}(Y)+B_{2,k-1}^2\|_F\\
    \leq\|f^{k-1}_{W_1}(Y)-f^{k-1}_{W_2}(Y)\|_F+2(\|G_{1,k-1}^1\|_{2\rightarrow2}\|f^{k-1}_{W_1}(Y)-f^{k-1}_{W_2}(Y)\|_F\\
    +\|f_{W_2}^{k-1}(Y)\|_F\|G_{2,k-1}^1-G_{1,k-1}^1\|_{2\rightarrow2}+\|B_{2,k-1}^1-B_{1,k-1}^1\|_F)\\
    +2(\|G_{1,k-1}^2\|_{2\rightarrow2}\|f^{k-1}_{W_1}(Y)-f^{k-1}_{W_2}(Y)\|_F+\|f_{W_2}^{k-1}(Y)\|_F\|G_{2,k-1}^2-G_{1,k-1}^2\|_{2\rightarrow2}\\
    +\underbrace{\|B_{2,k-1}^2-B_{1,k-1}^2\|_F}_{=0\text{, since it is not parameter-dependent}}),
\end{gather*}
Consequently,
\begin{gather}
    \|f^k_{W_1}(Y)-f^k_{W_2}(Y)\|_F\leq\|f^{k-1}_{W_1}(Y)-f^{k-1}_{W_2}(Y)\|_F(2\|G_{1,k-1}^1\|_{2\rightarrow2}\notag\\
    \label{est1}
    +2\|G_{1,k-1}^2\|_{2\rightarrow2}+1)+2\|f^{k-1}_{W_2}(Y)\|_F(\|G_{1,k-1}^2-G_{2,k-1}^2\|_{2\rightarrow2}\\
    +\|G_{1,k-1}^1-G_{2,k-1}^1\|_{2\rightarrow2})+2\|B_{2,k-1}^1-B_{1,k-1}^1\|_F.\notag
\end{gather}
For simplification, we will treat separately some terms in \eqref{est1}.
\smallskip

\underline{Estimation of $2\|G_{1,k-1}^1\|_{2\rightarrow2}+2\|G_{1,k-1}^2\|_{2\rightarrow2}+1$}. Due to Lemma \ref{g12bound} and the assumption that $W_1,W_2\in\mathcal{B}_\Lambda$, we have
\begin{equation}\label{p1}
    2\|G_{1,k-1}^1\|_{2\rightarrow2}
    +2\|G_{1,k-1}^2\|_{2\rightarrow2}+1\leq\Gamma_{k-1},
\end{equation}
with $\{\Gamma_k\}_{k\geq0}$ defined as in Lemma \ref{g12bound}.
\medskip

\underline{Estimation of $\|G^1_{2,k-1}-G^1_{1,k-1}\|_{2\rightarrow2}+\|G^2_{2,k-1}-G^2_{1,k-1}\|_{2\rightarrow2}$}.
\begin{gather}
    \|G^1_{2,k-1}-G^1_{1,k-1}\|_{2\rightarrow2}+\|G^2_{2,k-1}-G^2_{1,k-1}\|_{2\rightarrow2}\notag\\
    \leq\big[\mu^{-1}t^1_{k-1}\|W_2W_2^T-W_1W^T_1\|_{2\rightarrow2}+\mu^{-1}t^1_{k-1}\|A\|_{2\rightarrow2}\|W_2-W_1\|_{2\rightarrow2}\notag\\
    +(1-\theta_{k-1})\theta^{-1}_{k-1}\mu^{-1}t^1_{k-1}\|W_2W_2^T-W_1W^T_1\|_{2\rightarrow2}\notag\\
    +(1-\theta_{k-1})\theta^{-1}_{k-1}\mu^{-1}t^1_{k-1}\|A\|_{2\rightarrow2}\|W_2-W_1\|_{2\rightarrow2}\big]\notag\\
    +\big[\mu^{-1}t^2_{k-1}\|A\|_{2\rightarrow2}\|W_2-W_1\|_{2\rightarrow2}\notag\\
    \label{gprosum}
    +(1-\theta_{k-1})\theta^{-1}_{k-1}\mu^{-1}t^2_{k-1}\|A\|_{2\rightarrow2}\|W_2-W_1\|_{2\rightarrow2}\big]
\end{gather}
We define $c_{1,k}=\theta^{-1}_{k}\mu^{-1}t^1_{k}$, $c_{2,k}=\theta^{-1}_{k}\mu^{-1}t^2_{k}$, for $k\geq-1$, and $c_{1,-1}=c_{2,-1}=0$ (this holds due to the assumption that $t_{-1}^1=t_{-1}^2=0$). The previous statements and the assumption that $0<\{\theta_k\}_{k\geq-1}\leq1$, yield for \eqref{gprosum}
\begin{gather}
    \label{gsum}
    \|G^1_{2,k-1}-G^1_{1,k-1}\|_{2\rightarrow2}+\|G^2_{2,k-1}-G^2_{1,k-1}\|_{2\rightarrow2}\\
    \leq c_{1,k-1}\|W_2W_2^T-W_1W^T_1\|_{2\rightarrow2}+\|A\|_{2\rightarrow2}\|W_2-W_1\|_{2\rightarrow2}(c_{1,k-1}+c_{2,k-1}).\notag
\end{gather}
Furthermore,
\begin{gather}
    \|W_2W_2^T-W_1W^T_1\|_{2\rightarrow2}=\|W_2W_2^T-W_1W_2^T+W_1W_2^T-W_1W^T_1\|_{2\rightarrow2}\notag\\
    \label{estw}
    \leq\|W_2\|_{2\rightarrow2}\|W_2-W_1\|_{2\rightarrow2}+\|W_1\|_{2\rightarrow2}\|W_2-W_1\|_{2\rightarrow2}\implies\notag
\end{gather}
\begin{equation}
    \|W_2W_2^T-W_1W^T_1\|_{2\rightarrow2}\leq2\Lambda\|W_2-W_1\|_{2\rightarrow2}.
\end{equation}
Hence, we substitute the latter into \eqref{gsum}, yielding
\begin{gather}
    \|G^1_{2,k-1}-G^1_{1,k-1}\|_{2\rightarrow2}+\|G^2_{2,k-1}-G^2_{1,k-1}\|_{2\rightarrow2}\notag\\
    \leq2\Lambda c_{1,k-1}\|W_2-W_1\|_{2\rightarrow2}+\|A\|_{2\rightarrow2}\|W_2-W_1\|_{2\rightarrow2}(c_{1,k-1}+c_{2,k-1})\notag\implies
\end{gather}
\begin{equation}\label{gsumfinal}
    \begin{split}
        \|&G^1_{2,k-1}-G^1_{1,k-1}\|_{2\rightarrow2}+\|G^2_{2,k-1}-G^2_{1,k-1}\|_{2\rightarrow2}\\
    &\leq(2\Lambda c_{1,k-1}+\|A\|_{2\rightarrow2}(c_{1,k-1}+c_{2,k-1}))\|W_2-W_1\|_{2\rightarrow2}.
    \end{split}
\end{equation}
\medskip

\underline{Estimation of $\|B_{2,k-1}^1-B_{1,k-1}^1\|_F$}.
\begin{align}
    \|B_{2,k-1}^1-B_{1,k-1}^1\|_F&\leq\theta^{-1}_{k-1}t^1_{k-1}\|X_0\|_F\|W_2-W_1\|_{2\rightarrow2}\notag\\
    \label{estb}
    &\leq\mu c_{1,k-1}\|A\|_{2\rightarrow2}\|Y\|_F\|W_2-W_1\|_{2\rightarrow2}.
\end{align}
Now, we substitute \eqref{p1}, \eqref{gsumfinal}, \eqref{estb} into \eqref{est1} and apply Lemma \ref{fbound}, to get
\begin{gather*}
    \|f^k_{W_2}(Y)-f^k_{W_1}(Y)\|_F\leq\Gamma_{k-1}\|f^{k-1}_{W_2}(V)-f^{k-1}_{W_1}(V)\|_F\\
    +2\|f_{W_2}^{k-1}(V)\|_F\big[2\Lambda c_{1,k-1}+\|A\|_{2\rightarrow2}(c_{1,k-1}+c_{2,k-1})\big]\|W_2-W_1\|_{2\rightarrow2}\\
    +2\mu c_{1,k-1}\|X_0\|_F\|W_2-W_1\|_{2\rightarrow2}\\
    \leq\Gamma_{k-1}\|f^{k-1}_{W_2}(V)-f^{k-1}_{W_1}(V)\|_F+\big[2\Delta_{k-1}(2\Lambda c_{1,k-1}+\|A\|_{2\rightarrow2}(c_{1,k-1}+c_{2,k-1}))\\
    +2\mu c_{1,k-1}\|A\|_{2\rightarrow2}\|Y\|_F\big]\|W_2-W_1\|_{2\rightarrow2},
\end{gather*}
where
\begin{equation}
    \begin{split}
        \Delta_k=&2\mu\|Y||_F\Bigg[\sum_{i=0}^{k-1}\Bigg(\bigg(\|A\|_{2\rightarrow2}(c_{1,i-1}\Lambda+c_{2,i-1}\|A\|_{2\rightarrow2})+c_{2,i-1}\bigg)\\
        &\cdot\prod_{j=i}^{k-1}\Gamma_j\Bigg)+\|A\|_{2\rightarrow2}(c_{1,k-1}\Lambda+c_{2,k-1}\|A\|_{2\rightarrow2})+ c_{2,k-1}\Bigg],\\
        \Delta_0=&0.
    \end{split}
\end{equation}
Now, we set
\begin{equation*}
    E_k=2\Delta_{k-1}(2\Lambda c_{1,k-1}+\|A\|_{2\rightarrow2}(c_{1,k-1}+c_{2,k-1}))+2\mu c_{1,k-1}\|A\|_{2\rightarrow2}\|Y\|_F,
\end{equation*}
thus
\begin{equation}\label{gamma-e}
    \begin{split}
        \|f^k_{W_2}(Y)-f^k_{W_1}(Y)\|_F\leq&\Gamma_{k-1}\|f^{k-1}_{W_2}(V)-f^{k-1}_{W_1}(V)\|_F\\
        &+E_k\|W_2-W_1\|_{2\rightarrow2}
    \end{split}
\end{equation}
Using the abbreviations defined by $\Gamma_k,\,\Delta_k,\,E_k$, the general formula for $K_L$ in \eqref{fpert} is
\begin{align}\label{klgeneral}
    K_L=\sum_{k=1}^L\left(\max_{0\leq i\leq L-1}\Gamma_i\right)^{L-k}E_k\qquad\text{for}\quad L\geq1.
    \end{align}
Based on \eqref{gamma-e}, we prove via induction that \eqref{fpert} holds for any number of layers $L\geq1$, with $K_{L}$ given by \eqref{klgeneral}. For $L=1$, we can directly calculate $K_1$:
\begin{gather*}
    \|f^1_{W_2}(Y)-f^1_{W_1}(Y)\|_F\leq2t_0^1\|A\|_{2\rightarrow2}\|Y\|_F\|W_2-W_1\|_{2\rightarrow2}\\
    =2\underbrace{\theta_0^{-1}t_0^1}_{=1}\|A\|_{2\rightarrow2}\|Y\|_F\|W_2-W_1\|_{2\rightarrow2}=2\mu c_{1,0}\|A\|_{2\rightarrow2}\|Y\|_F\|W_2-W_1\|_{2\rightarrow2},
\end{gather*}
so that $2\mu c_{1,0}\|A\|_{2\rightarrow2}\|Y\|_F=E_1=K_1$ as claimed in \eqref{klgeneral}.
Now, let us assume \eqref{fpert} holds for some $L\in\mathbb{N}$. Then, applying the estimate that appears in \eqref{gamma-e} for $L+1$:
\begin{gather*}
    \|f^{L+1}_{W_2}(Y)-f^{L+1}_{W_1}(Y)\|_F\leq\Gamma_L\|f^L_{W_2}(Y)-f^L_{W_1}(Y)\|_F+E_{L+1}\|W_2-W_1\|_F\\
    \leq(\Gamma_LK_L+E_{L+1})\|W_2-W_1\|_F\leq[(\max_{0\leq i\leq L}\Gamma_i)K_L+E_{L+1}]\|W_2-W_1\|_F\\
    =[(\max_{0\leq i\leq L}\Gamma_i)\sum_{k=1}^L(\max_{0\leq i\leq L-1}\Gamma_i)^{L-k}E_{k}+E_{L+1}]\|W_2-W_1\|_F\\
    \leq\left(\sum_{k=1}^L(\max_{0\leq i\leq L}\Gamma_i)^{k}E_{k}+(\max_{0\leq i\leq L}\Gamma_i)^0E_{L+1}\right)\|W_2-W_1\|_F\\
    =\left(\sum_{k=1}^{L+1}(\max_{0\leq i\leq L}\Gamma_i)^{L+1-k}E_{k}\right)\|W_2-W_1\|_F=K_{L+1}\|W_2-W_1\|_F.
\end{gather*}
We successfully calculated the desired $K_L$. Now, under the additional assumptions $\{c_{1,k}\}_{k\geq0},\{c_{2,k}\}_{k\geq0}\leq1$, we may apply the ``in particular'' part of Lemma \ref{fbound} on \eqref{kl}, to obtain a simplified upper bound on $K_L$. Thus, for $\gamma$ and $\zeta_k$ defined in Lemmata \ref{g12bound} and \ref{fbound} respectively, we have
\begin{align}
    K_L\leq&2\mu\|Y\|_F\Bigg[\mu^{-1}\|A\|_{2\rightarrow2}+\sum_{k=2}^L\Bigg(\gamma^{L-k}\bigg(4\|A\|_{2\rightarrow2}((\Lambda+\|A\|_{2\rightarrow2})+1)\notag\\
    &\cdot(\Lambda+\|A\|_{2\rightarrow2})\bigg)(\zeta_{k-1}+1)+\|A\|_{2\rightarrow2}\Bigg)\Bigg]\notag.
\end{align}
In the latter, we set $Z=4\|A\|_{2\rightarrow2}((\Lambda+\|A\|_{2\rightarrow2})+1)(\Lambda+\|A\|_{2\rightarrow2})$ for the sake of brevity, so that
\begin{align}
    K_L\leq2\mu\|Y\|_F\Bigg[\|A\|_{2\rightarrow2}(\mu^{-1}+(L-1))+Z\gamma^L\sum_{k=2}^L\bigg(\gamma^{-k}(\zeta_{k-1}+1)\bigg)\Bigg].\notag
\end{align}
Now, plugging \eqref{zeta} in the latter yields
\begin{align}
    K_L\leq&2\mu\|Y\|_F\Bigg[\|A\|_{2\rightarrow2}(\mu^{-1}+(L-1))+Z\gamma^L\sum_{k=2}^L\gamma^{-k}\left(\frac{\gamma^{k-1}-1}{\gamma-1}+1\right)\Bigg]\notag\\
    =&2\mu\|Y\|_F\Bigg[\|A\|_{2\rightarrow2}(\mu^{-1}+(L-1))+\frac{Z\gamma^L}{\gamma-1}\sum_{k=2}^L\left(\frac{1}{\gamma}-\frac{2}{\gamma^k}+\frac{\gamma}{\gamma^k}\right)\Bigg]\notag\\
    =&2\mu\|Y\|_F\Bigg[\|A\|_{2\rightarrow2}(\mu^{-1}+(L-1))+\frac{Z\gamma^L}{\gamma-1}\left(\frac{L-2}{\gamma}+(\gamma-2)\sum_{k=2}^L\frac{1}{\gamma^k}\right)\Bigg]\notag\\
    =&2\mu\|Y\|_F\Bigg[\|A\|_{2\rightarrow2}(\mu^{-1}+(L-1))+\frac{Z\gamma^L}{\gamma-1}\left(\frac{L-2}{\gamma}+(\gamma-2)\frac{\gamma^{1-L}-1}{\gamma^{-1}-1}\right)\Bigg]\notag\\
    =&2\mu\|Y\|_F\Bigg[\|A\|_{2\rightarrow2}(\mu^{-1}+(L-1))\notag\\
    &+\frac{Z}{\gamma-1}\left(\frac{\gamma^L(L-2)}{\gamma}+\frac{\gamma(\gamma-2)(\gamma^L-\gamma)}{\gamma-1}\right)\Bigg]\notag\\
    <&2\mu\|Y\|_F\Bigg[\|A\|_{2\rightarrow2}(\mu^{-1}+(L-1)\notag\\
    &+Z\frac{\Bigg(\gamma^L\bigg((L-2)+\gamma^2(\gamma-2)\bigg)-\gamma^2(\gamma-2)\Bigg)}{\gamma-1}\Bigg].
\end{align}
The proof follows.
\end{proof}

It is time to prove the Lipschitzness of the decoder defined in \eqref{decoder}.

\begin{corollary}
Let $h\in\mathcal{H}^L$ defined as in \eqref{hypothesis} with $L\geq2$ and dictionary $W\in\mathcal{B}_\Lambda$. Then, for any $W_1,W_2\in\mathcal{B}_\Lambda$, we have:
\begin{equation}
    \|\psi(\phi(f_{W_2}^L(Y)))-\psi(\phi(f_{W_1}^L(Y)))\|_F\leq\mu^{-1}(\Lambda+\|A\|_{2\rightarrow2})K_L\|W_2-W_1\|_F,
\end{equation}
with $K_L$ as in Theorem \ref{perturbation}.
\end{corollary}

\begin{proof}
By definition, $\psi$ is a 1-Lipschitz function. Moreover, as an affine map, $\phi$ is Lipschitz continuous with Lipschitz constant $\text{Lip}_{\phi}=\|\Phi\|_{2\rightarrow2}$, where $\Phi$ is defined in \eqref{phi}. We evaluate $\|\Phi\|_{2\rightarrow2}$:
\begin{gather*}
    \|\Phi\|_{2\rightarrow2}\leq\mu^{-1}\theta_L\|W\|_{2\rightarrow2}+\mu^{-1}\theta_L\|A\|_{2\rightarrow2}+\mu^{-1}(1-\theta_L)\|W\|_{2\rightarrow2}\\
    +\mu^{-1}(1-\theta_L)\|A\|_{2\rightarrow2}\leq\mu^{-1}(\Lambda+\|A\|_{2\rightarrow2}).
\end{gather*}
Combining the previous estimates, we get
\begin{gather*}
    \|\psi(\phi(f_{W_2}^L(Y)))-\psi(\phi(f_{W_1}^L(Y)))\|_F\leq\|\phi(f_{W_2}^L(Y))-\phi(f_{W_1}^L(Y))\|_F\\
    \leq\|\Phi\|_{2\rightarrow2}\|f_{W_2}^L(Y)-f_{W_1}^L(Y)\|_F\leq\mu^{-1}(\Lambda+\|A\|_{2\rightarrow2})K_L\|W_2-W_1\|_F,
\end{gather*}
where in the last inequality we used Theorem \ref{perturbation}.
\end{proof}
\section{Main Results}
In this Section, we estimate the desired generalization error. To do so, we combine the tools we developed so far with a chaining technique.
\subsection{Covering numbers and Dudley's inequality}
For a fixed number of layers $L\in\mathbb{N}$, we define the set $\mathcal{M}\subset\mathbb{R}^{n\times s}$ corresponding to the hypothesis class $\mathcal{H}^L$:
\begin{equation}
    \begin{split}         \mathcal{M}&:=\{(h(y_1)|h(y_2)|\dots|h(y_s))\in\mathbb{R}^{n\times s}:\ h\in \mathcal{H}^{L}\}\\
    &=\{\psi(\phi((f^L_W(Y)))\in\mathbb{R}^{n\times s}:\ W\in\mathcal{B}_\Lambda\}.
    \end{split}
\end{equation}
The column elements of each matrix in $\mathcal{M}$ are the reconstructions given by a decoder $h\in\mathcal{H}^L$ when applied to the measurements $y_i$. Since $\mathcal{M}$ is parameterized by $W$ like $\mathcal{H}^L$ is, we may rewrite \eqref{contradem} as
\begin{align}\label{rlh}
    \mathcal{R}_s(l\circ\mathcal{H}^L)&\leq\sqrt{2}(2\mathrm{B_{in}}+2\mathrm{B_{out}})\mathcal{R}_s(\mathcal{H}^L)=\sqrt{2}(2\mathrm{B_{in}}+2\mathrm{B_{out}})\mathcal{R}_s(\mathcal{M})\notag\\
    &=\sqrt{2}(2\mathrm{B_{in}}+2\mathrm{B_{out}})\mathbb{E}\sup_{M\in\mathcal{M}}\frac{1}{s}\sum_{i=1}^{s}\sum_{k=1}^n\epsilon_{ik}M_{ik}.
\end{align}
Thus, we are left with estimating the Rademacher process
\begin{equation}
    \mathbb{E}\sup_{M\in\mathcal{M}}\frac{1}{s}\sum_{i=1}^{s}\sum_{k=1}^n\epsilon_{ik}M_{ik}.
\end{equation}
The latter has subgaussian increments, hence we use Dudley's inequality \cite{rf} to upper bound it in terms of the covering numbers of the set $\mathcal{M}$.
\begin{theorem}[Dudley's inequality]\label{dudley}
Let $(X_t)_{t\in T}$ be a centered subgaussian process with radius $\Delta(T)=\sup_{t\in T}\sqrt{\mathbb{E}|X_t|^2}$. Then, 
\begin{equation}
    \mathbb{E}\sup_{t\in T}X_t\leq4\sqrt{2}\int_0^{\Delta(T)/2}\sqrt{\log(\mathcal{N}(T,d,u))}du.
\end{equation}
\end{theorem}
\noindent Therefore, we first calculate the radius of $\mathcal{M}$, i.e.
\begin{align}
    \Delta(\mathcal{M})&=\sup_{h\in\mathcal{H}^L}\sqrt{\mathbb{E}\left(\sum_{i=1}^{s}\sum_{k=1}^n\epsilon_{ik}h_k(y_i)\right)^2}\leq\sup_{h\in\mathcal{H}^L}\sqrt{\mathbb{E}\sum_{i=1}^{s}\sum_{k=1}^n\epsilon_{ik}(h_k(y_i))^2}\notag\\\label{radius}
    &\leq\sup_{h\in\mathcal{H}^L}\sqrt{\sum_{i=1}^{s}\|h(y_i)\|_2^2}\overset{\eqref{boundedout}}{\leq}\sqrt{s}B_{\mathrm{out}}.
\end{align}
Now, applying Theorem \ref{dudley} to \eqref{rlh} yields
\begin{equation}\label{dudleys}
    \mathcal{R}_s(l\circ\mathcal{H}^L)\leq\frac{16(\mathrm{B_{in}}+\mathrm{B_{out}})}{s}\int_0^{\sqrt{s}B_{\mathrm{out}}/2}\sqrt{\log\mathcal{N}(\mathcal{M},\|\cdot\|_F,\varepsilon)}d\varepsilon.
\end{equation}
According to \eqref{dudleys}, we need to estimate $\mathcal{N}(\mathcal{M},\|\cdot\|_F,\varepsilon)$. Towards that end, we state the following Lemma.
\begin{lemma}\label{cover}
The covering number of $B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(\Lambda)=\{X\in\mathbb{R}^{N\times n}:\,\|X\|_{2\rightarrow2}\leq \Lambda,\,\Lambda>0\}$ satisfies the following for any $\varepsilon>0$:
\begin{equation}
    \mathcal{N}(B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(\Lambda),\|\cdot\|_{2\rightarrow2},\varepsilon)\leq\left(1+\frac{2\Lambda}{\varepsilon}\right)^{Nn}.
\end{equation}
\end{lemma}
\begin{proof}
For $|\cdot|$ denoting the volume in $\mathbb{R}^{N\times n}$, the following is an adaptation of a well-known result (Proposition 4.2.12 of \cite{vershynin}) connecting covering numbers and volume in $\mathbb{R}^{N\times n}$:
\begin{align*}
    \mathcal{N}(B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(\Lambda),\|\cdot\|_{2\rightarrow2},\varepsilon)&\leq\frac{|B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(\Lambda)+(\varepsilon/2)B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(1)|}{|(\varepsilon/2)B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(1)|}\\
    &=\frac{|\Lambda\cdot B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(1)+(\varepsilon/2)B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(1)|}{|(\varepsilon/2)B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(1)|}\\
    &=\frac{|(\Lambda+\varepsilon/2)B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(1)|}{|(\varepsilon/2)B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(1)|}.
\end{align*}
Hence,
\begin{align*}
    \mathcal{N}(B_{\|\cdot\|_{2\rightarrow2}}^{N\times n}(\Lambda),\|\cdot\|_{2\rightarrow2},\varepsilon)&\leq\left(1+\frac{2\Lambda}{\varepsilon}\right)^{Nn}.\qedhere
\end{align*}
\end{proof}

\begin{proposition}\label{coveringm}
The following estimate holds for the covering numbers of $\mathcal{M}$:
\begin{equation}\label{coverestim}
    \mathcal{N}(\mathcal{M},\|\cdot\|_F,\varepsilon)\leq\left(1+\frac{2\Lambda \mu^{-1}(\Lambda+\|A\|_{2\rightarrow2})K_L}{\varepsilon}\right)^{Nn}.
\end{equation}
\end{proposition}
\begin{proof}
We first consider the set $\Omega=\{W:\,W\in\mathcal{B}_\Lambda\}\subset\mathbb{R}^{N\times n}$. Then, due to Lemma \ref{cover}, we can upper bound the covering numbers of $\Omega$ as follows:
\begin{equation}
    \mathcal{N}(\Omega,\|\cdot\|_{2\rightarrow2},\varepsilon)\leq\left(1+\frac{2\Lambda }{\varepsilon}\right)^{Nn}.
\end{equation}
Now, for the covering numbers of $\mathcal{M}$ we have
\begin{align*}
    \mathcal{N}(\mathcal{M},\|\cdot\|_F,\varepsilon)&\leq\mathcal{N}(\mu^{-1}(\Lambda+\|A\|_{2\rightarrow2})K_L\Omega,\|\cdot\|_{2\rightarrow2},\varepsilon)\\
    &=\mathcal{N}(\Omega,\|\cdot\|_{2\rightarrow2},\varepsilon/(\mu^{-1}(\Lambda+\|A\|_{2\rightarrow2})K_L))\\
    &\leq\left(1+\frac{2\Lambda \mu^{-1}(\Lambda+\|A\|_{2\rightarrow2})K_L}{\varepsilon}\right)^{Nn}.\qedhere
\end{align*}
\end{proof}

\subsection{Generalization error bounds for DECONET}
We are now in position to prove our main results, that estimate the generalization error of DECONET.
\begin{theorem}\label{gentheorem}
Let $\mathcal{H}^L$ be the hypothesis class defined in \eqref{hypothesis}. With probability at least $1-\delta$, for all $h\in\mathcal{H}^L$, the generalization error is bounded as
\begin{align}
    \mathcal{L}(h)\leq\mathcal{\hat{L}}(h)+&8(B_{\mathrm{in}}+B_{\mathrm{out}})B_{\mathrm{out}}\sqrt{\frac{Nn}{s}}\sqrt{\log\left(e\left(1+\frac{4\mu^{-1}\Lambda(\Lambda+\|A\|_{2\rightarrow2})K_L}{\sqrt{s}B_{\mathrm{out}}}\right)\right)}\notag\\ 
    +&4(B_{\mathrm{in}}+B_{\mathrm{out}})^2\sqrt{\frac{2\log(4/\delta)}{s}},
\end{align}
with $K_L$ defined in \eqref{kl}.
\end{theorem}
\begin{proof}
We apply Proposition \ref{coveringm} to \eqref{dudleys}, yielding:
\begin{align*}
    \mathcal{R}_s(l\circ\mathcal{H}^L)&\leq\frac{16(\mathrm{B_{in}}+\mathrm{B_{out}})}{s}\int_0^{\frac{\sqrt{s}B_{\mathrm{out}}}{2}}\sqrt{\log\mathcal{N}(\mathcal{M},\|\cdot\|_F,\varepsilon)}d\varepsilon\\
    &\leq\frac{16(\mathrm{B_{in}}+\mathrm{B_{out}})}{s}\int_0^{\frac{\sqrt{s}B_{\mathrm{out}}}{2}}\sqrt{Nn\log\left(1+\frac{2\Lambda\mu^{-1}(\Lambda+\|A\|_{2\rightarrow2})K_L}{\varepsilon}\right)}d\varepsilon\\
    &\leq8(\mathrm{B_{in}}+\mathrm{B_{out}})B_{\mathrm{out}}\sqrt{\frac{Nn}{s}}\sqrt{\log\left(e\left(1+\frac{4\Lambda\mu^{-1}(\Lambda+\|A\|_{2\rightarrow2})K_L}{\sqrt{s}B_{\mathrm{out}}}\right)\right)},
\end{align*}
where in the last step we used the inequality
\begin{equation*}
    \int_0^a\sqrt{\log\left(1+\frac{b}{t}\right)}dt\leq a\sqrt{\log(e(1+b/a))},\qquad a,b>0.
\end{equation*}
Now, we use Theorem \ref{radem} with the upper bound $c=(B_{\mathrm{in}}+B_{\mathrm{out}})^2$ for the loss function $\|\cdot\|_2^2$, and the proof follows. 
\end{proof}

Similarly to Section 4, we may further assume that the ratios $\{t^1_k/\mu\theta_k\}_{k\geq0}$, $\{t^2_k/\mu\theta_k\}_{k\geq0}$ are upper bounded by 1, so we obtain

\begin{corollary}\label{gencor}
Let $\mathcal{H}^L$ be the hypothesis class defined in \eqref{hypothesis} and assume that $\{t^1_k/\mu\theta_k\}_{k\geq0}$, $\{t^2_k/\mu\theta_k\}_{k\geq0}\leq1$. With probability at least $1-\delta$, for all $h\in\mathcal{H}^L$, the generalization error is bounded as
\begin{align}
    \mathcal{L}(h)<&\mathcal{\hat{L}}(h)+8(B_{\mathrm{in}}+B_{\mathrm{out}})\Bigg(B_{\mathrm{out}}\sqrt{\frac{Nn}{s}}\\
    &\cdot\sqrt{\log\left(e\left(1+\frac{\|Y\|_F(p+q(L-1))+r\kappa_L)}{\sqrt{s}B_{\mathrm{out}}}\right)\right)}+\sqrt{\frac{2\log(4/\delta)}{s}}\Bigg),\notag
\end{align}
with $\kappa_L$ as in Theorem \ref{perturbation} and $p,q,r>0$ constants depending on $\|A\|_{2\rightarrow2}, \Lambda, \mu$. 
\end{corollary}
\begin{proof}
The estimate easily follows from Theorems \ref{perturbation} - \ref{gentheorem}, if we set 
\begin{align}
    p:=&8\mu^{-1}\Lambda(\Lambda+\|A\|_{2\rightarrow2})\|A\|_{2\rightarrow2}\\
    q:=&8\Lambda(\Lambda+\|A\|_{2\rightarrow2})\\
    r:=&32\Lambda(\Lambda+\|A\|_{2\rightarrow2})\|A\|_{2\rightarrow2}((\Lambda+\|A\|_{2\rightarrow2})+1)(\Lambda+\|A\|_{2\rightarrow2}).
\end{align}
\end{proof}

The last and most important result of this section is given below.

\begin{theorem}\label{gengentheorem}
Let $\mathcal{H}^L$ be the hypothesis class defined in \eqref{hypothesis}. Assume there exist pair-samples $\{(x_i,y_i)\}_{i=1}^s$, with $y_i=Ax_i+e$, $\|e\|_2\leq\varepsilon$, for some $\varepsilon>0$, that are drawn i.i.d. according to an unknown distribution $\mathcal{D}$, and that it holds $\|y_i\|_2\leq\mathrm{B_{in}}$ almost surely with $\mathrm{B_{in}}=\mathrm{B_{out}}$ in \eqref{psi}. Let us further assume that for step sizes $\{t^1_k\}_{k\geq0}$, $\{t^2_k\}_{k\geq0}>0$, step size multiplier $0<\{\theta_k\}_{k\geq0}\leq1$ and smoothing parameter $\mu>0$, we have $\{t^1_k/\mu\theta_k\}_{k\geq0}$, $\{t^2_k/\mu\theta_k\}_{k\geq0}\leq1$. Then with probability at least $1-\delta$, for all $h\in\mathcal{H}^L$, the generalization error is bounded as
\begin{align}
    \mathcal{L}(h)<&\mathcal{\hat{L}}(h)+16B^2_{\mathrm{out}}\sqrt{\frac{Nn}{s}}\sqrt{\log\left(e\left(1+p+q(L-1))+r\kappa_L\right)\right)}\notag\\
    \label{genbound}
    &+16B_{\mathrm{out}}\sqrt{\frac{2\log(4/\delta)}{s}}.
\end{align}
\end{theorem}
\begin{proof}
The result is obtained by applying the previous Corollary and \eqref{boundedin}.
\end{proof}

\begin{remark}\label{asymptoticL}
If we consider the dependence of the generalization error bound \eqref{genbound} only on $L,N,s$ and treat all other terms as constants, then according to Theorem \ref{perturbation} for $\kappa_L$, we have
\begin{align}\label{genbigoh}
    |\mathcal{L}(h)-\mathcal{\hat{L}}(h)|\lesssim\sqrt{\frac{NL}{s}}.
\end{align}

\end{remark}

\section{Numerical Experiments}
In this Section, we are interested in examining whether our theory regarding the generalization error of DECONET is consistent with real-world CS paradigms.
\subsection{Settings}
We train and test DECONET on two real-world image datasets: MNIST (60000 training and 10000 test $28\times28$ image examples) and CIFAR10 \cite{cifar} (50000 training and 10000 test $32\times32$ coloured image examples). For the CIFAR10 dataset, we transform the images into grayscale ones. We consider the vectorized form of the images. This means that for the ambient dimension $n$ of an arbitrary vectorized image $x\in\mathbb{R}^n$ we have $n=28^2=784$ (if $x$ is a $28\times28$ MNIST image), or $n=32^2=1024$ (if $x$ is a $32\times32$ CIFAR10 image). We examine DECONET with varying number of layers $L$. We consider two CS ratios, i.e. $m/n=25\%$ and $m/n=50\%$. We choose a random Gaussian measurement matrix $A\in\mathbb{R}^{m\times n}$ and appropriately normalize it, i.e. $\Tilde{A}=A/\sqrt{m}$. We add zero-mean Gaussian noise $e$ with standard deviation $\mathrm{std}=10^{-4}$ to the measurements $y$, so that $y=\tilde{A}x+e$. We set $\varepsilon=\|y-\tilde{A}x\|_2$ and $x_0=A^Ty$, which are standard algorithmic setups. We take different values of $N$ and perform He (normal) initialization for $W\in\mathbb{R}^{N\times n}$. We set $\mu=100$, initial step sizes $t_0^1=t_0^2=1$ and step size multiplier $\theta_0=1$. The authors in \cite{tfocs:6} define $t_k^2=\alpha t_k^1$, where $\alpha=\frac{\|W\|_2^2}{\|A\|_2^2}$, but in the current we extend their result by adapting separately $t_k^1$, $t_k^2$, with $t_{k}^1=\alpha t_{k-1}^1$, $t_k^2=\beta t_{k-1}^2$, $k=1,\dots,L$, respectively, where $(\alpha,\beta)\in(0,1)\times(0,1)$. We set $(\alpha,\beta)=(0.5,0.3)$ and $(\alpha,\beta)=(0.7,0.5)$ for the MNIST and CIFAR10 datasets, respectively. Moreover, we consider the following update rule for the step size multiplier: $\theta_k=\theta_{k-1}\cdot\theta_s$, $k=1,\dots,L$, where $\theta_s=\frac{1-\sqrt{\mu/\Tilde{L}}}{1+\sqrt{\mu/\Tilde{L}}}$ and $\Tilde{L}$ is an upper bound on the smoothing parameter $\mu$; we set $\Tilde{L}=1000$. For MNIST and CIFAR10, we train all networks with learning rate $\eta=10^{-2}$ and $\eta=10^{-3}$, respectively. All networks are implemented in PyTorch \cite{pytorch} and trained using the \emph{Adam} optimizer \cite{adam}, with batch size $128$. For our experiments, we report the \emph{test MSE} defined by
\begin{equation}\label{testmse}
    \mathcal{L}_{test}=\frac{1}{d}\sum_{i=1}^d\|h(\tilde{y}_i)-\tilde{x}_i\|_2^2,
\end{equation}
where $\mathcal{D}=\{(\tilde{y}_i,\tilde{x}_i)\}_{i=1}^d$ 
is a set of $d$ test data, not used in the training phase, and the \emph{empirical generalization error} (EGE) defined by
\begin{equation}\label{genmse}
    \mathcal{L}_{gen}=|\mathcal{L}_{test}-\mathcal{L}_{train}|,
\end{equation}
where $\mathcal{L}_{train}$ is the train MSE defined in \eqref{trainmse}. Due to the fact that test MSE approximates the true loss, we use \eqref{genmse} -- which can be explicitly computed -- to approximate \eqref{generror}. We train all networks, on all datasets, employing an early stopping technique \cite{earlystop} with respect to \eqref{genmse}. Training and test of each DECONET's variant are repeated at least 10 times and the results are averaged over the runs. Finally, we compare the EGE of DECONET to the EGE of another deep unfolding network, that is ISTA-net, proposed in \cite{istagen}. ISTA-net jointly learns a decoder for CS and an orthogonal dictionary for sparsification. Under this setting, ISTA-net solves the CS problem employing synthesis sparsity instead of analysis sparsity. We aim to see how the EGE is affected by each of the two sparsity models. For ISTA-net, we set the best hyper-parameters proposed by the original authors.
\subsection{Experiments}
We test DECONET on MNIST and CIFAR10 datasets under multiple experimental scenarios.

\subsubsection{Fixed CS ratio for 10- and 50-layer DECONET with varying $N/n$}
We examine the performance of 10- and 50-layer DECONET for a fixed $25\%$ CS ratio, redundancy ratio $N/n$ varying in the set $\{10,15,20,25, 30\}$ and report the results in Fig.~\ref{mnist25} (MNIST) and \ref{cifar25} (CIFAR10). The top subplots of Fig.~\ref{mnist25} and \ref{cifar25} illustrate how the test MSEs, achieved by 10- and 50-layer DECONET, drop as $L$ and $N/n$ increase. The decays seem reasonable, if one considers a standard analysis CS scenario: a) the performance and reconstruction quality provided by the analysis-$l_1$ algorithm typically benefit from the (high) redundancy offered by the involved analysis operator b) more iterations/layers result to a higher reconstruction quality. The bottom subplots of Fig.~\ref{mnist25} and \ref{cifar25} demonstrate the increment of the EGEs, achieved by 10- and 50-layer DECONET, as both $L$ and $N/n$ increase. This behaviour confirms our theoretical result depicted in Theorem \ref{gengentheorem} (and Remark \ref{asymptoticL}).

\begin{figure*}[h!
\centering
\begin{subfigure}
{.8\textwidth}
\centering
\includegraphics[width=\textwidth]{mnist_25.jpeg}
\captionsetup{justification=centering}
\caption{MNIST}
\label{mnist25}
\end{subfigure}\vfill
\begin{subfigure}
{.8\textwidth}
\centering
\includegraphics[width=\textwidth]{cifar_25.jpeg}
\captionsetup{justification=centering}
\caption{CIFAR10}
\label{cifar25}
\end{subfigure}
\caption{Average test MSEs (top subplots) and empirical generalization errors (bottom subplots) for 10- and 50-layer DECONET, with $25\%$ CS ratio, tested on (a) MNIST and (b) CIFAR10 datasets.}
\label{redplots}
\end{figure*}

\subsubsection{Fixed CS ratio with varying $L$ and $N$}

We examine the generalization ability of DECONET for $m=n/2$, with increasing number of layers $L$, under different choices of $N$. Inspired by frames with redundancy ratio $N/n\notin\mathbb{N}$ \cite{shearletsframes}, we consider $N$ of the form
\begin{equation}\label{Nformula}
    N=pn^2+q,\qquad p,q\in\mathbb{N}.
\end{equation}
We report the results in Fig.~\ref{mnist50} (MNIST) and \ref{cifar50} (CIFAR10). Similarly to Section 6.2.1, we observe that the empirical generalization error increases in $L$ and $N$, for both datasets. Even though the upper bound in \eqref{genbound} depends on many terms, the empirical generalization error appears to grow at the rate of $\sqrt{N}$. The behaviour of DECONET again conforms with our theoretical results presented in Theorem~\ref{gengentheorem}. One may also notice that -- in general -- we choose different $N$ for each of the two datasets; this is simply due to \eqref{Nformula}, i.e. $N$ depends on the vectorized ambient dimension $n$.

\begin{figure*}[h!
\centering
\begin{subfigure}
{.8\textwidth}
\centering
\includegraphics[width=\textwidth]{mnist50.jpeg}
\captionsetup{justification=centering}
\caption{Empirical generalization error for $m=n/2$ measurements on MNIST, with alternating $N$}
\label{mnist50}
\end{subfigure}\vfill
\begin{subfigure}
{.8\textwidth}
\centering
\includegraphics[width=\textwidth]{cifar50.jpeg}
\captionsetup{justification=centering}
\caption{Empirical generalization error for $m=n/2$ measurements on CIFAR10, with alternating $N$}
\label{cifar50}
\end{subfigure}
\caption{Performance plots for DECONET with $50\%$ CS ratio, tested on MNIST (top) and CIFAR10 (bottom) datasets.}
\label{Lplots}
\end{figure*}

\subsubsection{DECONET vs ISTA-net}
In this set of experiments, we aim to see how analysis and synthesis sparsity models affect the generalization ability of CS-oriented unfolding networks. Towards this end, we compare the proposed DECONET's decoder to ISTA-net's decoder, for 10, 30 and 50 layers, with $25\%$ and $50\%$ CS ratio, and fixed $N=37145$ for our learnable analysis operator, on all datasets. We report the corresponding empirical generalization errors in Table~\ref{decoders}. First of all, we see that our proposed decoder outperforms the ISTA-net's decoder, consistently for both datasets. This behaviour indicates that learning a redundant sparsifier instead of an orthogonal one, improves the performance of a CS-oriented unfolding network. Second, albeit the authors of \cite{istagen} prove that the generalization error of ISTA-net increases in $L$, we notice that this is not the case for either of the two datasets and either of the two CS ratios. Our theoretical result, instead, seems to align with the experiments, since the EGE of DECONET increases as $L$ also increases. Moreover, for the CIFAR10 dataset, the EGE of DECONET decreases as the number of measurements $m$ increases, but this decay is not explained by our theoretical result, which does not account for $m$. Hence, this observation is in need of future mathematical explanation. 

\begin{table*}[h!]
    \centering
    \scalebox{0.7}{\begin{tabular}{||c|c|c|c|c|c|c|c|c||}\hline
         &  \multicolumn{8}{|c|}{$25\%$ CS ratio}\\
         \hline
         Dataset & \multicolumn{3}{|c|}{MNIST} & \multicolumn{3}{|c|}{CIFAR10} \\
         \hline
         \diagbox{Decoder}{Layers} & $L=10$ & $L=30$ & $L=50$ & $L=10$ & $L=30$ & $L=50$ \\ \hline
         DECONET & \bf0.000033 & \bf0.000314 & \bf0.000446 & \bf0.000066 & \bf0.000100 & \bf0.0000128 \\
         \hline
         LISTA & 0.013408 & 0.007874 & 0.005937 & 0.007165 & 0.004120 & 0.003085\\
         \hline
    \end{tabular}}
    
    \medskip
    
    \scalebox{0.7}{\begin{tabular}{||c|c|c|c|c|c|c|c|c||}\hline
         &  \multicolumn{7}{|c|}{$50\%$ CS ratio}\\
         \hline
         Dataset & \multicolumn{3}{|c|}{MNIST} & \multicolumn{3}{|c|}{CIFAR10} \\
         \hline
         \diagbox{Decoder}{Layers} & $L=10$ & $L=30$ & $L=50$ & $L=10$ & $L=30$ & $L=50$ \\ \hline
         DECONET & \bf0.000131 & \bf0.000241 & \bf0.000465 & \bf0.000024 & \bf0.000063 & \bf0.000111\\
         \hline
         LISTA & 0.016547 & 0.009311 & 0.007157 & 0.009274 & 0.005576 & 0.004377\\
         \hline
    \end{tabular}}
    \caption{Empirical generalization error for 10-, 30- and 50-layer decoders (all datasets), with fixed $N=37145$. Bold letters indicate the best performance between the two decoders.}
    \label{decoders}
\end{table*}

\section{Conclusion and Future Work}

In the present paper, we derived a new deep unfolding network dubbed DECONET, based on a well-known optimization algorithm solving analysis-sparsity-based Compressed Sensing. DECONET jointly learns a CS decoder and a redundant analysis operator serving as a sparsifying transform. We introduced the hypothesis class consisting of all the functions/decoders that DECONET can implement and upper bounded its corresponding Rademacher complexity using chaining techniques. In the end, we derived generalization error bounds for DECONET, in terms of the aforementioned Rademacher complexity estimate. Our generalization error bounds depend on the number of layers and the redundancy of the learned analysis operator. To the best of our knowledge, we are the first to explain the generalization ability of an unfolded network that jointly learns a CS decoder and a redundant sparsifying transform. Moreover, an important aspect of our derived generalization error bounds is that they scale like the square root of the redundancy of the learnable sparsifier. From the experimental perspective, we evaluated our derived theory by testing DECONET on two real-world  image datasets and compared it to a state-of-the-art learnable synthesis-sparsity-based decoder. Our experiments confirmed that the generalization errors achieved by DECONET on both datasets scaled like our theoretical generalization error bounds. Moreover, our analysis decoder outperforms the baseline in terms of the generalization error, consistently for both datasets and different choices of layers and CS ratios. As a future direction, we could enlarge the hypothesis class, i.e., we could also learn some of the parameters of the original iterative algorithm. Furthermore, it would be interesting to check if the redundancy of the learnable sparsifier introduces some kind of implicit regularization into DECONET.

\section*{Acknowledgements}
The author would like to thank H. Rauhut for his valuable advice and inspiring conversations around the subject covered in this paper.

\printbibliography
\end{document}