\section{Introduction}\label{sec:introduction}
Recent experimental and theoretical results \citep{zhang2021understanding, belkin2019reconciling} have demonstrated that the classical understanding of overparameterized machine learning (ML) models requires further examination. One model that has been studied extensively is the random features (RF) mode \citep{rahimi2007random}, which is closely related to overparameterized neural networks \citep{daniely2016toward, daniely2017sgd, jacot2018neural, liu2021random, bach2017equivalence}. In this paper, we examine an extension of the RF model, which we call the deep RF (DRF) model, being equivalent to a deep NN, but only trained in the output layer. We consider the asymptotic regime, where the number of data points, model parameters, and input dimension grow infinite at constant ratio \citep{belkin2020two, hastie2019surprises, bartlettLongLugosiTsigler_2020, Bartlett30063} and give exact expressions that characterize the deep RF model in terms of training and generalization error.

Our analysis consists of two key steps. First, we prove universality, i.e. we demonstrate that the DRF model is asymptotically equivalent to a deep Gaussian surrogate model, matching the original model in the first and second moments, at each layer \citep{panahi2017universal, oymak2018universality}. Universality for the 1-layer RF model has previously been proven, e.g. in \citep{HuUniversalityLaws}. We make use of a different proof technique to extend these results to arbitrary many layers and introduce a new Gaussian surrogate model for DRF. This universality result alleviates the general difficulty of analyzing RF or DRF models, as the non Gaussian features are in general not amenable to stardard analysis techniques such as comparison theorem \citep{gordon1985some, thrampoulidis2015gaussian}, Gaussian widths \citep{chandrasekaran2012convex} or replica methods \citep{mezard1987spin}.

Having established universality, we then make use of the Convex Gaussian Min Max Theorem (CGMT) \citep{thrampoulidis2015gaussian} to study DRFs. This theorem allows us to consider an alternative optimization problem with the same asymptotic statistics, and is a popular tool in the analysis of the asymptotic regime \citep{bosch2021double, chang2020provable, dhifallah2020precise, thrampoulidis2015regularized, Loureiro2021Learning, Bosch2022DoubleDescent}. We make use of a recursive application of the CGMT \citep{Bosch2022DoubleDescent} to obtain asymptotic expressions for square loss functions with arbitrary convex regularization for $L$-layer DRF models.

\section{Related Works}\label{sec:relatedworks}

The random features (RF) \citep{rahimi2007random} model has been extensively examined in the asymptotic regime, under a multitude of conditions. For an incomplete list see \citep{hastie2019surprises, mei2019generalization, montanari2019generalization, goldt2020gaussian, goldt2019modeling, gerace2020generalisation, dhifallah2020precise, ghorbani2021linearized, Bosch2022DoubleDescent}. In the case of ridge regression \citep{louart2018random, mei2019generalization} exact expression for the training and generalization error can be established. In other cases, exact analysis is difficult. It was observed by many authors \citep{mei2019generalization, hastie2019surprises, goldt2019modeling, gerace2020generalisation, goldt2020gaussian} that a Gaussian surrogate model that matched the first and second moments had asymptotically equivalent statistics. A concrete proof of RF universality is given in \cite{hu2022universality}.
We utilize Lindeberg's approach \citep{lindeberg1922neue} to demonstrate universality of DRF. This approach has been used to prove universality results in many other optimization problems \citep{korada2011applications,panahi2017universal, montanari2017universality,oymak2018universality, abbasi2019universality}. \cite{hu2022universality} prove a central limit theorem between random features and their Gaussian equivalent features as a key step in demonstrating universality. We make use of a different proof technique, by instead considering the problem in a dual space, where we may directly bound the difference between the leave one out iterates.

Beside RF,  universality has been demonstrated for many other models \citep{goldt2020gaussian, seddik2020random, dhifallah2021inherent, loureiro2021capturing, gerace2022gaussian}. Recently, \cite{montanari2022universality} gave a proof for the universality of empirical risk minimization for not necessarily convex loss and regularization functions. Their result also assume that a central limit theorem similar to \citep{hu2022universality} holds.

Subject to Gaussian features, the CGMT \citep{gordon1985some, thrampoulidis2015gaussian} is a powerful tool in the determination of the asymptotic performance \citep{Loureiro2021Learning, dhifallah2020precise, thrampoulidis2015regularized, chang2020provable, bosch2021double, Bosch2022DoubleDescent}. The CGMT determines an alternative, asymptotically equivalent optimization problem in statistical properties. In the case of correlated features, such as in the RF or DRF model, the alternative optimization still remains intractable. This issue is resolved in \citep{Bosch2022DoubleDescent} by applying the CGMT twice. Relying on the particular structure of the DRF covariance matrices, we extend the method of \citep{Bosch2022DoubleDescent} where the CGMT is applied    recursively to determine a nested scalar optimization that is asymptotically equivalent to the DRF model.

The covariance matrices for the Gaussian surrogate model that we obtain are similar in structure to the kernel matrices given in \cite{lee2017deep}. The authors demonstrate an exact equivalence between an infinitely wide deep NN and a Gaussian Process with covariance kernels that are recursively defined in a similar manner to the ones discussed in this paper. However, \citep{lee2017deep} consider networks of fixed size but infinite width, while we consider the asymptotic regime, where the number of data points and the input dimensions grow as well, hence maintaining a relatively narrower network.

\subsection{Paper Outline}\label{sec:introduction:subsec:PaperOutline}
In section \ref{sec:setup}, we introduce the DRF problem and its Gaussian surrogate, and express the necessary assumptions for our results to hold.
In section \ref{sec:Universality}, we prove the main universality theorem of this paper. Our proof takes two steps, first proving universality of a single layer, and subsequently using an inductive argument to extend this result to a full DRF problem.
In section \ref{sec:CGMTanalysis}, we give an alterative scalar optimization problem derived by means of the CGMT, that is asymptotically equivalent to the DRF problem subject to square loss and arbitrary, strongly convex regularization. We demonstrate experimentally the veracity of the determined expressions.



\section{Setup and Assumptions}\label{sec:setup}
\subsection{Random Feature Model and Preliminaries}
We consider a supervised learning setup with a dataset $\mathcal{D} = \lbrace(\mathbf{x}_k, y_k) \in \mathbb{R}^{d}\times \mathbb{R}\rbrace_{k=1}^n$. To find a relationship between the data points $\mathbf{x}_k$ and the labels $y_k$, we consider function of the following form 
\begin{eqnarray}\label{eq:GenericFeatureMap}
    Y_{\bm{\theta}}(\mathbf{x}) = \frac{1}{\sqrt{p}}\bm{\theta}^T \mathcal{F}(\mathbf{x}), \qquad \bm{\theta}\in\mathbb{R}^p,
\end{eqnarray}
where $\mathcal{F}:\mathbb{R}^d \rightarrow \mathbb{R}^p$ is a given mapping of the data, called a \textit{feature map}. We note that $Y_{\bm{\theta}}$ is dependent upon the choice of the vector $\bm{\theta}$ by a linear relation. This shows the main advantage of \eqref{eq:GenericFeatureMap}: while $Y_{\bm{\theta}}$ can represent nonlinear functions, selecting $\bm{\theta}$ amounts to a linear regression task. To find the optimal value of $\bm{\theta}$, denote $\mathbf{f}_i=\mathcal{F}(\mathbf{x}_i)$ and take $\bm{F}$ as a matrix with $\mathbf{f}_i$ as columns. We consider the empirical risk minimization framework and the following optimization problem:
\begin{equation} \label{eq:OptimalPoint}
    \hat\bm{\theta}=\hat{\bm{\theta}}(\bm{F}) = \arg\min_{\bm{\theta} \in \mathbb{R}^p} \frac{1}{n}\sum_{i = 1}^n \ell\left(\frac{1}{\sqrt{p}}\bm{\theta}^T\mathbf{f}_i, y_i\right) + R(\bm{\theta}).
\end{equation}
 Here, $\ell(\cdot, \cdot)$ is a loss function, and $R(\bm{\theta})$ is a regularization function. To measure the performance of $\hat{\bm{\theta}}$ we make use of the two common metric for supervised learning, that being the training error $\mathcal{E}_{train}(\bm{F})$, i.e the optimal value in \eqref{eq:OptimalPoint}, and
the generalization error
\begin{equation}\label{eq:GeneralizationLoss}
    \mathcal{E}_{gen}(\bm{F}) = \mathbb{E}\left[\ell\left(  \frac{1}{\sqrt{p}}\hat{\bm{\theta}}^T\mathcal{F}(\mathbf{x}_{new}),y_{new}\right)\right],
\end{equation}
where the expectation is taken over $(\mathbf{x}_{new}, y_{new})$, a new datapoint drawn from the same distribution as the dataset $\mathcal{D}$. 

The main purpose of this paper is to obtain exact asymptotic expressions for the supervised learning metrics $\mathcal{E}_{train},\mathcal{E}_{gen}$ and other properties of $\hat{\bm{\theta}}$, when the feature map $\mathcal{F}$ is a deep random feature, generalizing the \textit{random features maps} \citep{rahimi2007random}. To define the deep features, we remind the (shallow) random features are given by $\phi(\mathbf{x}, \mathbf{w}_j):= \sigma(\mathbf{w}_j^T\mathbf{x})$, for $j=1,2,\ldots,p$, where $\sigma$ is an activation function, and $\mathbf{w}_j \sim \mathcal{N}(0, \frac{1}{d}\mathbf{I}_d)$ are a set of random weights. In vector form, we express these relations as $\bm{\phi}(\mathbf{x}, \mathbf{W}) = (\sigma(\mathbf{w}_j^T\mathbf{x}))_{j=1}^p$,
where the matrix $\mathbf{W}$ has rows $\mathbf{w}_j$. Then, the deep random features are
given through the following recursion: For $p_0 = d, p_1, p_2, \ldots,p_L \in \mathbb{N}$, we define the matrices $\mathbf{W}^{(l)} \in \mathbb{R}^{p_l\times p_{l-1}}$ for $l= 1, \ldots, L$, each having independent rows $\mathbf{w}^{(l)}_j \sim \mathcal{N}(0, \frac{1}{p_{l-1}}\mathbf{I})$. Letting $\mathbf{x}^{(0)}:= \mathbf{x}$ we define
\begin{equation}\label{eq:RecursiveRandomFeatures}
    \mathbf{x}^{(l)} = \bm{\phi}(\mathbf{x}^{(l-1)}, \mathbf{W}^{(l)}) = (\sigma(\mathbf{w}_j^{(l)T}\mathbf{x}^{(l-1)}))_{j=1}^{p_l}, \qquad l = 1, \ldots, L,
\end{equation}
\subsection{Necessary Assumptions}
Our results rely on the following assumptions:

\begin{enumerate}
    \item[A1] For some universal positive constants $\mu, M$, the regularization function $R$ is $\mu$-strongly convex and $M$-smooth with $M$-bounded third derivative in tensor (operator) norm. Moreover $R(\bm{\theta})$ is minimized at $\bm{\theta} = \bm{0}$
    \item[A2] $\ell$ is a $\frac{1}{Cn}-$ strongly convex function in the first argument and its third derivative with respect to the first argument is bounded by $Cn$ for some constant $C$. Moreover, there exists a vector $\bm{\alpha} = (\alpha_k)$ called \textit{isolated predictions} satisfying: $ \alpha_k \in \arg\min_{\alpha} \ell(\alpha, y_k)$, and $\norm{\bm{\alpha}}_2 \leq C\sqrt{n}$ for a fixed constant $C$.
    \item[A3] The activation function $\sigma$ is an odd function applied element wise, with bounded derivatives. Furthermore, let $g_1, g_2$ be Gaussian variables distributed as
    \begin{equation}
        \begin{bmatrix}
        g_1\\
        g_2
        \end{bmatrix} \sim \mathcal{N}\left(\bm{0}, \begin{bmatrix}
            \alpha_1 & \rho \\
            \rho & \alpha_2
        \end{bmatrix}\right).
    \end{equation}
    Let the functions $\eta_1(\alpha_1, \alpha_2, \rho) = \mathbb{E}[\sigma(g_1)\sigma(g_2)]$ and $\eta_2(\alpha_1) = \mathbb{E}[\sigma^2(g_1)]$. Then $\eta_1, \eta_2$ should be thrice differentiable at $\alpha_1 = \alpha_2 = 1$ and $\rho = 0$
    \item[A4] The dimensions of the number of data points $n$, the size of the input $p_0 = d$ and the size of subsequent layers $p_l$, where $l=1,\ldots, L$ all grow to infinity at fixed ratios. We denote this by $n\sim p_0\sim\cdots\sim p_L$, where $a \sim b$ is defined to mean that $\frac{a}{b} \xrightarrow[a, b\rightarrow\infty]{} C$ for some constant $C$.
    \item[A5] For each layer, $l$, the weight matrix $\mathbf{W}^{(l)}\in\mathbb{R}^{p_l \times p_{l-1}} = [\mathbf{w}_1^{(l)}\ \mathbf{w}_2^{(l)}\ \cdots\ \mathbf{w}_{p_l}^{(l)}]^T$ are independent Gaussian variables $\mathbf{w}_i^{(l)} \overset{i.i.d}{\sim}\mathcal{N}(0, \frac{1}{p_{l-1}}\mathbf{I}_{p_{l-1}})$ for $1\leq i \leq p_l$. Furthermore, $\mathbf{W}^{(l)}$ are independent of the input variables $\mathbf{x}$. 
\end{enumerate}

\begin{remark} For assumption 2, we note that the strong convexity assumption on the loss function becomes less restrictive as $n$ grows. In the asymptotic limit, the strong convexity is no longer a significant requirement. Furthermore, the isolated prediction vectors exist and the condition is satisfied immediately if the loss function is minimized at the labels, i.e. it is minimized at the point $\ell(y_k, y_k) < \infty$.
\end{remark}

\begin{remark}
For assumption 3, we note that the condition holds for the majority of loss function used in practice including $\tanh$ and the error function. Furthermore, if oddness is dropped, the assumption on the functions $\eta_1$ and $\eta_2$ are additionally satisfied for functions like ReLU, sigmoids, and Gaussian activations. However, we require oddness.   
\end{remark}


Finally we impose a condition upon the input vectors $\mathbf{x}_i$:

\begin{definition}\label{def:regularity}
Let $d\sim n$, we call a set $\lbrace \mathbf{x}_k\in\mathbb{R}^{d} \rbrace_{k=1}^n$ \textit{regular} if
\begin{enumerate}
    \item Letting $\mathbf{X} = [\mathbf{x}_1\ \mathbf{x}_2\ \cdots\ \mathbf{x}_n]$, there is a constant $c<\infty$ such that $\frac{1}{\sqrt{n}}\norm{\mathbf{X}}_{\mathrm{op}} < c$
    \item It holds that
    \begin{equation}
        \max_{i, j}\left|\frac{1}{d}\mathbf{x}_i^T\mathbf{x}_j - \delta_{ij} \right| \leq \frac{\mathrm{polylog}\  n}{\sqrt{n}},
    \end{equation}
    where $\delta_{ij}$ is the Kronecker delta.
\end{enumerate}
\end{definition}
Note that the first condition for regularity is trivially satisfied for finite $n$, however the condition must also hold for a fixed $c$ in the asymptotic limit. Further, note that regularity is exhibited by $\mathbf{x}$ being Gaussian with high probability.


\section{Universality}\label{sec:Universality}
In the case of a single layer, it has been proven \citep{HuUniversalityLaws} that the following Guassian feature map has asymptotically equivalent statistics to the random features given in section~\ref{sec:setup}
\begin{equation}
    \Tilde{\bm{\phi}}(\mathbf{x}, \mathbf{W}) = \rho_1\mathbf{W}\mathbf{x} + \rho_2\mathbf{g},
\end{equation}
where $\mathbf{g}$ is a standard normal vector, and $\rho_1, \rho_2$ are constants depending only on the activation function, given by
$
     \rho_1  = \mathbb{E}[\sigma'(z)], \ \rho_2 = \sqrt{\mathbb{E}[\sigma^2(z)] - \rho_0^2 + \rho_1^2},\  z\sim\mathcal{N}(0, 1).
$
Similarly we define a \textit{deep Gaussian equivalent feature map}, recursively. We define $\bm{\gamma}^{(0)} = \mathbf{x}$, and then define
\begin{equation}\label{eq:GaussianRandomFeatures}
    \bm{\gamma}^{(l)} = \Tilde{\bm{\phi}}_l(\bm{\gamma}^{(l-1)}, \mathbf{W}^{(l)}) := \rho_{1,l}\mathbf{W}^{(l)}\bm{\gamma}^{(l-1)} + \rho_{2,l}\mathbf{g}^{(l)}, \qquad l = 1, \ldots, L,
\end{equation}
where $\mathbf{g}^{(l)}$ is an independent standard normal vector of dimension $p_l$, and the constants $\rho_{1,l},\rho_{2,l}$ are recursively defined as
$    \rho_{1,l} = \mathbb{E}[\sigma'(\alpha_{l-1}z)], \ \rho_{2,l} = \sqrt{\mathbb{E}[\sigma^2(\alpha_{l-1}z)] - \alpha_{l-1}^2\rho_{1,l}^2}$, $ \ z\sim\mathcal{N}(0, 1)$.
Here, $\alpha_l$ are constants given by the following recursive definition:
 $       \alpha_0 = 1, \  \alpha_{l} = \sqrt{\rho_{1,l}^2\alpha_{l-1}^2 + \rho_{2,l}^2}$.
Now, we consider the following two feature mappings for an input vector $\mathbf{x}^{(0)}$:
\begin{eqnarray}\label{eq:featureMapsDef}
    \mathcal{F}(\mathbf{x}^{(0)}) = \mathbf{x}^{(l)},\qquad \mathcal{G}(\mathbf{x}^{(0)}) = \bm{\gamma}^{(l)},
\end{eqnarray}
where $\mathbf{x}^{(l)}$ and $\bm{\gamma}^{(l)}$ are given in \eqref{eq:RecursiveRandomFeatures} and  \eqref{eq:GaussianRandomFeatures}, respectively.


\subsection{Revisiting Universality of a Single Layer}\label{subsec:UniversalitySingleLayer}
The proof of universality of deep random features is a specific application of a universality theorem for a single layer, which we derive in this section. This result is more general than the previous studies such as \cite{hu2022universality}. In the subsequent section, we shall demonstrate how the universality of deep random features follows from these results.

Let $\phi_{j}:\mathbb{R}^{d}\times \Omega_j\rightarrow \mathbb{R}$ for $j=1,2,\ldots, p$ be a random feature map, where $\Omega_j$ is a sample space equipped with an arbitrary probability measure, such that $\phi_j(\cdot, \omega)$ for any $\omega\in \Omega_j$ is a particular realization of the feature map. Let $\Omega = \Omega_1\times \Omega_2 \times \cdots \Omega_p$ be a product space equipped with the product measure, and let $\bm{\phi}:\mathbb{R}^{d}\times \Omega\rightarrow \mathbb{R}^p$ represent the vector of random features, such that $\bm{\phi}(\mathbf{x}, \bm{\omega}) = (\phi_{j}(\mathbf{x}, \omega_j))_j$, where $\bm{\omega} = (\omega_j)\in \Omega$ is a realization. 

Next we consider a $m\times p$ matrix $\mathbf{D}$ with columns $\mathbf{d}_j\in\mathbb{R}^m$, which we call a \textit{synthesis dictionary}. We define the re-represented random feature vectors $\mathbf{f}:\mathbb{R}^d\times \Omega \rightarrow \mathbb{R}^{m}$ given by
\begin{equation}\label{eq:RerepresentedRF}
    \mathbf{f}(\mathbf{x}, \bm{\omega}):= \sum_{j=1}^p\mathbf{d}_{j}\phi_j(\mathbf{x}, \omega_j) = \mathbf{D}\bm{\phi}(\mathbf{x}, \omega).
\end{equation}
We note that if $m = p$ we can choose $\mathbf{D} = \mathbf{I}_p$ and retain the original set of random features. However, re-representing the features is necessary for the proof of the deep random features case. We will drop the argument $\bm{\omega}$ when there is no risk of confusion, and denote $\bm{\phi}(\mathbf{x})$, $\mathbf{f}(\mathbf{x})$ as the random features and their re-representation. Similarly let $\bm{\phi}_k = \bm{\phi}(\mathbf{x}_k)$ and $\mathbf{f}_k = \mathbf{f}(\mathbf{x}_k)$ for $k = 1,\ldots, n$ which are random vectors. Finally, let $\bm{\Phi}$ and $\bm{F}$ be the matrices with $(\bm{\phi}_k), (\mathbf{f}_k)$ as columns. We assume that the random features are centered:
\begin{equation}\label{eq:RFcentered}
    \mathbb{E}_{\bm{\omega}}[\bm{\phi}(\mathbf{x}_k, \bm{\omega})] = \bm{0}\qquad k =1, 2, \ldots, n.
\end{equation}
We further define the data kernel matrices $\bm{K}_{j} = (K_{j, kl})_{kl}$ where $\bm{K}_j$ is the covariance matrix of the $j$th row of $\bm{\Phi}$, given by
\begin{equation}\label{eq:RFcovarianceMatrices}
    K_{j, kl} = \mathbb{E}_{\omega_j}[\phi(\mathbf{x}_k, \omega_j)\phi(\mathbf{x}_l, \omega_j)].
\end{equation}

Next, we introduce a $p\times n$ Gaussian matrix $\bm{\Gamma}$ with independent rows, and where the $j$th row is distributed by $\mathcal{N}(0, \bm{K}'_j)$. We note that if $\bm{K}_j = \bm{K}_j'$ that $\bm{\Phi}$ and $\bm{\Gamma}$ have the same first and second moments amongst their elements. We then define $\mathbf{G} = \mathbf{D}\bm{\Gamma}$ and let $\mathbf{g}_k$ be the $k$th column of $\mathbf{G}$.

Before stating the main theorem for this section we state the conditions on the dataset and matrices $\bm{K}_j$ and $\mathbf{D}$ that must hold. We shall show in the next section that these conditions are satisfied in the case of deep random features. We remind the reader of the definition of a sub-Gaussian vector:

\begin{definition}\label{def:subgaussianity}
    We say that a random vector $\mathbf{u} = (u_k)\in \mathbb{R}^n$ is $\tau-$sub-Gaussian if for any unit vector $\mathbf{a} = (a_k)\in\mathbb{R}^n$ the variable $A = \mathbf{a}^T\mathbf{u}$ is sub-Gaussian with parameter $\tau$, i.e.
   
     $   \mathbb{E}\left[e^{\lambda A} \right] \leq e^{\frac{\tau^2\lambda^2}{2}}$
   
    for all $\lambda \in \mathbb{R}$.
\end{definition}
We state the following requisite conditions:
\begin{enumerate}
    \item[B1] There exists a positive constant $C$ such that for all $j$, it holds that $\norm{\bm{K}_j}_{\mathrm{op}}\leq C$ and the $j$th random feature vector $\bm{\phi}^j = \lbrace \phi(\mathbf{x}_k, \omega_j) \rbrace_{k}$ is $C$-sub-Gaussian
    \item[B2] There exists a positive constant $C$ such that $\norm{\mathbf{D}}_{\mathrm{op}} \leq C$.
\end{enumerate}
These assumptions must hold for all values of $n, d, p, m$ and must continue to hold when they grow asymptotically.
Subject to these conditions we state the following theorem that demonstrates universality.
\begin{theorem}\label{thm:GenericUniversality}
    Suppose that assumptions A1, A2, B1 and B2 hold, and that $n\sim p \sim m$. Then,
    \begin{enumerate}
    \item For any real function $\psi$ with bounded first, second, and third derivatives, there exists a constant $c < \infty$ such that
    \begin{equation}
        \left|\mathbb{E}\psi(\mathcal{E}_{train}(\bm{F})) - \mathbb{E}\psi(\mathcal{E}_{train}(\mathbf{G})) \right| \leq  \frac{c}{n}\sum_{j =1}^p\norm{\bm{K}_j - \bm{K}_j'}_{\mathrm{op}} + \frac{c}{\sqrt{n}}.
    \end{equation}
    \item Let $\hat{\bm{\theta}}_F$ and $\hat{\bm{\theta}}_G$ be the optimal points for the optimization \eqref{eq:OptimalPoint} for $\bm{F}$ and $\mathbf{G}$ respectively. For any bounded function $h:\mathbb{R}^{m}\rightarrow\mathbb{R}$ with bounded second and third derivatives (in tensor norm), where the bounds are constant in $n, p, m$. There exists a constant $c<\infty$ such that
    \begin{equation}
        \left|\mathbb{E} h\left(\hat{\bm{\theta}}_F\right) - \mathbb{E} h\left(\hat{\bm{\theta}}_G\right) \right| \leq \frac{c}{n}\sum_{j =1}^p\norm{\bm{K}_j - \bm{K}_j'}_{\mathrm{op}} + \frac{c}{\sqrt{n}}.
    \end{equation}
    \end{enumerate}
\end{theorem}

\subsubsection{Proof Sketch}\label{sec:UniversalityTheorem:subsec:proof}
The proof is based on an application of Lindebergs argument with respect to the random features $\mathbf{f}$ in a dual space. We consider the optimization problem given in \eqref{eq:OptimalPoint} for some generic map $\mathbf{Z}$ and note that by means of a splitting argument it may be expressed as
\begin{eqnarray}
    \mathcal{E}_{train}(\mathbf{Z}) = \min_{\bm{\theta}\in\mathbb{R}^{p}}\frac{1}{n}\sum_{k=1}^{n}\ell(\mathbf{z}_k^T\bm{\theta}, y_k) + R(\bm{\theta})\nonumber\\
    = \min_{\bm{\theta}\in\mathbb{R}^{p}, \bm{\alpha}\in\mathbb{R}^{n}}\max_{\mathbf{d}\in\mathbb{R}^{n}}\frac{1}{n}\left(\sum_{k=1}^n\ell(\alpha_k, y_k) + d_k(\alpha_k -\mathbf{z}_k^T\bm{\theta}) \right) + R(\bm{\theta})\nonumber\\
    = - \min_{\mathbf{d}\in\mathbb{R}^n}\frac{1}{n}\sum_{k=1}^n \ell^*(-d_k, y_k) + R^*\left(\frac{1}{n}\mathbf{Z}\mathbf{d}\right)
\end{eqnarray}
where $\ell^*$ and $R^*$ are the legendre transforms of $\ell$ and $R$ respectively. We note that by assumption A1 that $R^*$ is $\frac{1}{M}-$ strongly convex and $\frac{1}{\mu}$-smooth. We then proceed in defining a series of $\mathbf{Z}_r$ such that $\mathbf{Z}_0 = \bm{\Phi}$ and $\mathbf{Z}_p = \bm{\Gamma}$. We show that the difference of the optimal value in the dual space between $\psi(\mathcal{E}_{train}(\mathbf{Z}_r))$ and $\psi(\mathcal{E}_{train}(\mathbf{Z}_{r+1}))$ is bounded by the sum of a $O(\frac{1}{n^{3/2}})$ term and the difference in operator norm between $\frac cn\norm{\bm{K}_r - \bm{K}'_r}_{op}$, which allows us to bound the total difference as in the given result. 

For part two, we note that $R_{\epsilon}(\bm{\theta}) = R(\bm{\theta}) \pm \epsilon h(\bm{\theta})$ remains strongly convex for sufficiently small values of $\epsilon > 0$. As such, part 1 of the theorem holds for these cases. By bounding the difference in the values of $\mathcal{E}_{train}$ at $\epsilon>0$ and at $\epsilon = 0$ the bound on $h(\bm{\theta})$ may be obtained. The proof is given in full in appendix \ref{app:sec:UniversalityTheoremProof}.

\subsection{Multiple Layers}\label{subsec:multipleLayerUniversality}
In this section, we apply the results of the previous section to prove universality for DRF.
We shall consider the deep random features as given in eq \eqref{eq:featureMapsDef}. The proof of the equivalence relies on fixing all layers, except a single one, and demonstrating that the individual layer may be replaced by their Gaussian equivalent. This relies on an intermediate result, given in the following theorem, stating that the regularity of a dataset, as defined in definition \ref{def:regularity} is preserved under random feature mappings.
\begin{theorem}\label{thm:PreservationOfRegularity}
    Suppose that the set $\lbrace\mathbf{x}_i\in\mathbb{R}^d \rbrace_{i=1}^n$ is regular and assumption 3 holds. Then define $\mathbf{z}_i = \sigma(\mathbf{W}\mathbf{x}_i)$ where $\mathbf{W}$ is a $p\times d$ matrix and has independent rows distributed by $\mathcal{N}(0, \frac{1}{d}\mathbf{I})$. Then with probability higher than $1 - n^{-10}$ the set $\lbrace \mathbf{z}_i\rbrace_{i =1}^n$ is regular\footnote{The exponent of $n$ is arbitrary and can be replaced by any other number}.
\end{theorem}
The main consequence of this theorem is that for $L$ layers where $n\sim p_0\sim p_1\sim \cdots \sim p_L$, with a probability converging to 1, all dataset $\mathbf{X}^{(l)} = \lbrace \mathbf{x}_i^{(l)} \rbrace_{i=1}^n$ for $l=1,\ldots, L$ are regular\footnote{Here we assume that the numbers of layers $L$ is fixed, but it is simple to show that the argument also holds for $L=\mathrm{poly}(n)$}. Now, we can state the main result of this section, which demonstrates a slightly more generic version of universality for an $l$-layered deep random feature model.

\begin{theorem}\label{thm:DeepRFUniversality}
Suppose that $n \sim p_0 \sim \cdots \sim p_l$ and take $q = \mathcal{O}(n)$ and let assumption A1-A5 hold. For a fixed final layer $l$, define noise appended features $\Tilde{\mathbf{x}}^{(l)}_i, \tilde{\bm{\gamma}}_i^{(l)}$ as
\begin{equation}
    \tilde{\mathbf{x}}_i^{(l)} = \begin{bmatrix}
        \mathbf{x}_i^{(l)}\\
        \mathbf{v}_i^{(l)}
    \end{bmatrix}\qquad 
    \tilde{\bm{\gamma}}_i^{(l)} = \begin{bmatrix}
        \bm{\gamma}_i^{(l)}\\
        \mathbf{v}_i^{(l)}
    \end{bmatrix}
\end{equation}
where $\mathbf{v}_i^{(l)}\in\mathbb{R}^{q}$ are independent standard Gaussian vectors. Take a $m \times (q+p_l)$ dictionary $\mathbf{D}$, where $\norm{\mathbf{D}}_{\mathrm{op}} < c$ for some constant $c<\infty$ and define the re-represented features
\begin{equation}
    \mathbf{f}_i = \mathbf{D}\Tilde{\mathbf{x}}_i, \quad \mathbf{g}_i = \mathbf{D}\tilde{\bm{\gamma}}_i
\end{equation}
and let $\bm{F} = [\mathbf{f}_1\ \cdots\ \mathbf{f}_n]$ and let $\mathbf{G} = [\mathbf{g}_1\ \cdots\ \mathbf{g}_n]$ be their matrix representations. Then under the assumption that $\mathbf{X}$ is regular,
\begin{enumerate}
    \item For any real function $\psi$ with bounded first, second, and third derivatives, there exists a constant $c < \infty$ such that
    \begin{equation}
        \left|\mathbb{E}\psi(\mathcal{E}_{train}(\bm{F})) - \mathbb{E}\psi(\mathcal{E}_{train}(\mathbf{G})) \right| \leq \frac{\mathrm{polylog}\  n}{\sqrt{n}}
    \end{equation}
    \item Let $\hat{\bm{\theta}}_F$ and $\hat{\bm{\theta}}_G$ be the optimal solution of problem \eqref{eq:OptimalPoint} for $\bm{F}$ and $\mathbf{G}$. Then for any bounded function $h:\mathbb{R}^{p_{L}}\rightarrow\mathbb{R}$ with bounded second and third derivatives (in tensor norm), where the bounds are constant in $n, m, p_i$ for $0\leq i \leq l$. There exists a constant $c<\infty$ such that
    \begin{equation}
        \left|\mathbb{E} h\left(\hat{\bm{\theta}}_F\right) - \mathbb{E} h\left(\hat{\bm{\theta}}_G\right) \right| \leq \frac{\mathrm{polylog}\  n}{\sqrt{n}}
    \end{equation}
\end{enumerate}
\end{theorem}

Universality of the DRF problem follows directly from this theorem by choosing the final $L$th layer,  $q = 0$ and, hence adding no additional noise and $\mathbf{D} = \mathbf{I}_{p_L}$ such that no re-representation appears.

\subsubsection{Proof Sketch}\label{sec:RFuniversality:subsec:Proof}
The proof proceeds by means of induction. For the case that $l = 0$, ie a zero layer network the proof is immediate as $\mathbf{x}^{0} = \bm{\gamma}^{(0)}$. Assuming that the induction hypothesis holds for a layer $l-1$ we may consider layer $l$.

We make use of an intermediate results which may be found in the appendix. In theorem~\ref{thm:ReplacemetnOfCovarianceMatrix} we show that if the data set $\mathbf{x}^{(l)}$ is regular then covariance matrices of $\mathbf{x}^{(l)}$ and $\bm{\gamma}^{(l)}$ are bounded by $\frac{c\mathrm{polylog}\  n}{\sqrt{n}}$ for some constant $c$. Then, the proof proceeds in two steps: First, we consider an intermediate vector 
\begin{eqnarray}
    \bar{\bm{\gamma}}^{(l)}_i = \begin{bmatrix}
        \rho_{1, l}\mathbf{W}^{(l)}\mathbf{x}_i^{(l-1)} + \rho_{2,l}\mathbf{h}_i^{(l)}\\
        \mathbf{v}_i^{(l)}
    \end{bmatrix}.
\end{eqnarray}
We bound the performance difference ($\mathcal{E}_{train}$) between $\mathbf{x}^{(l)}$ and $\bar{\bm{\gamma}}^{(l)}$  by  theorem \ref{thm:DeepRFUniversality}. Second, we observe that the difference in performance  between $\bar{\bm{\gamma}}^{(l)}$ and $\bm{\gamma}^{(l)}$ depends only on the difference between $\mathbf{x}^{(l-1)}$ and $\bm{\gamma}^{(l-1)}$. As such, we may make use of the induction hypothesis to bound this difference. The full proof is given in appendix~\ref{app:sec:DeepRFUniversalityProof}.

\section{CGMT Analysis}\label{sec:CGMTanalysis}

Thanks to the universality results, we only require to analyze the deep Gaussian features $\bm{\gamma}_L$. Here, we present this analysis in one particular case where $\ell$ is the square loss, and the regularization function is generic. Additionally, we need to impose a model for the relationship between the labels $\mathbf{y}$ and the input variables $\mathbf{x}^{(0)}$, which we specifically assume to be independent standard normal vectors. For this we make the following definition
\begin{equation}\label{eq:ymodel}
    y_i = \mathbf{x}^{(L)T}_i\bm{\theta}^* + \nu_i,
\end{equation}
where $\bm{\theta}^*\in\mathbb{R}^{p_L}$ is the "true" relationship between the data and the parameters, $\nu_i\sim\mathcal{N}(0, \sigma^2_{\bm{\nu}}\mathbf{I})$ is noise, and $\mathbf{x}^{(L)}$ is defined in \eqref{eq:RecursiveRandomFeatures}. We let $\bm{\nu} = (\nu_i)_i$ and let $\mathbf{X}^{(L)} = [\mathbf{x}^{(L)}_1\ \mathbf{x}_2^{(L)}\ \cdots \mathbf{x}_n^{(L)}]$. Then, we consider the following optimization problem

\begin{equation}\label{eq:CGTMP1}
    P_1 = \min_{\bm{\theta}}\frac{1}{2n}\norm{\mathbf{y} - \mathbf{X}^{(L)}\bm{\theta}}_2^2 + R(\bm{\theta})=\min_{\mathbf{e}}\frac{1}{2n}\norm{\bm{\nu} - \mathbf{X}^{(L)}\mathbf{e}}_2^2 + R(\bm{\theta}^*+\mathbf{e}),
\end{equation}
where $\mathbf{e}=\bm{\theta}-\bm{\theta}^*$ and the optimal solutions are denoted by $\hat{\bm{\theta}}_1,\hat\mathbf{e}_1$. We similarly consider the Gaussian equivalent model defined in eq~\eqref{eq:GaussianRandomFeatures}. In this case, the data is  generated by
\begin{equation}\label{eq:ytildemodel}
    \Tilde{\mathbf{y}}_i = \bm{\gamma}^{(L)T}_i\bm{\theta}^* + \nu_i.
\end{equation}
Again, we let $\Tilde{\mathbf{X}}^{(L)} = [\bm{\gamma}_1^{(L)}\ \bm{\gamma}_2^{(L)}\ \cdots \bm{\gamma}_n^{(L)}]$ and define the Guassian equivalent optimization problem as
\begin{equation}\label{eq:CGTMP2}
    P_2 = \min_{\bm{\theta}}\frac{1}{2n}\norm{\tilde{\mathbf{y}} - \Tilde{\mathbf{X}}^{(L)}\bm{\theta}}_2^2 + R(\bm{\theta})= \min_{\mathbf{e}}\frac{1}{2n}\norm{\bm{\nu} - \Tilde{\mathbf{X}}^{(L)}\mathbf{e}}_2^2 + R(\bm{\theta}^*+\mathbf{e})
\end{equation}
with corresponding optimal solutions $\hat{\bm{\theta}}_2,\hat\mathbf{e}_2$. By applying theorem \ref{thm:DeepRFUniversality} to $\bm{\theta}\to\mathbf{e}$ and $\mathbf{y}\to\bm{\nu}$, we establish that  the statistics of $P_1$ and $P_2$ become weakly similar in the sense of their distributions.
Furthermore, for this particular choice of the relationship between the data and the labels the generalization error for the problem $P_1$ may be expressed as
\begin{eqnarray}
    \mathcal{E}_{gen}(\mathbf{e}) = \sigma_{\bm{\nu}^2} + \mathbf{e}^T\mathbb{E}[\tilde{\mathbf{x}}_{new}^{(L)T}\tilde{\mathbf{x}}_{new}^{(L)}]\mathbf{e}.
\end{eqnarray}
This function satisfies the conditions on the function $h(\mathbf{e})$. As such in this case, the generalization error is also universal.

For problem $P_2$, as the matrix $\Tilde{\mathbf{X}}$ is Gaussian, it may be analyzed by the CGMT, (see appendix theorem~\ref{app:thm:CGMT}), which gives an asymptotic equivalence to a second alternative problem is follows:

\begin{theorem}\label{thm:CGMTsquareloss}
Let $n\sim p_0\sim \cdots \sim p_L$ and let assumptions 1-5 hold true. Consider the following optimization problem
\begin{eqnarray}\label{eq:CGTMP3}
    P_3 =     \max_{\beta >0 }\min_{q} T_L + \max_{\xi_L > 0, \chi_L > 0} \min_{t_L > 0, k_L > 0} T_{L-1} + \min_{\mathbf{e}} \frac{a}{2pl}\norm{\mathbf{e}}^2 + \frac{b}{p_L}\mathbf{e}^T\mathbf{g} + R(\bm{\theta} + \bm{\theta}^*) + \nonumber\\   \max_{\xi_{L-1} > 0, \chi_{L-1} > 0} \min_{t_{L-1} > 0, k_{L-1} > 0} \cdots \max_{\xi_0 \geq 0, \chi_0 \geq 0} \min_{t_0 > 0, k_0 > 0} \sum_{i=1}^{L-2} T_l(\mathbf{e}) 
\end{eqnarray}
Where $T_L$ is a function of $\beta, q$; $T_{L-1}, a, b$ are functions of $\beta, q, \xi_L, \chi_L, t_L, k_L$;  $\mathbf{g}\in\mathbb{R}^{p_L}$ is a standard normal and $T_l$ are functions of $\mathbf{e} , \beta, q, \xi_i, \chi_i, t_i, k_i$ for $L \geq i \geq l$. The exact expressions for the functions $a, b, T_i$ are complicated and are given in the appendix equation \eqref{CGMTFullEquations}. 

Then, 
\begin{enumerate}
    \item Then the values of $P_2$ and $P_3$ become close, in sense that if $P_3$ converges to come value $c$ then $P_2$ will converge to the same value. 

    \item Let $\hat{\bm{\theta}}_3$ be the optimal point of $P_3$. Then for any bounded function $h:\mathbb{R}^{p_L}\rightarrow \mathbb{R}$ with bounded second and third derivatives (in tensor norm), where the bounds are constant in $n, p_i$. for $0\leq i\leq L$, then

    \begin{eqnarray}
        \Pr\left(|h(\hat{\mathbf{e}}_2) - h(\hat{\mathbf{e}}_3)| > \epsilon \right) \rightarrow 0 \quad \mathrm{as}\quad  n,p_0,\ldots,p_L \rightarrow\infty
    \end{eqnarray}

\end{enumerate}
\end{theorem}

 
A proof of this theorem may be found in the Appendix Section \ref{app:sec:CGMTproof}. Furthermore, if all layers, except the input have the same dimension $p$ the CGMT result can be simplified substantially, these results may be seen in theorem \ref{app:CGMT:allLayersSameSize}. It can be clearly seen that by the triangle inequality and the results of theorem \ref{thm:DeepRFUniversality} that $P_3$ and $P_1$ will similarly asymptotically become weakly similar; as will $h(\hat{\bm{\theta}}_3)$ and $h(\hat{\bm{\theta}}_1)$.

\subsection{Experimental Results}\label{sec:DepthBetter:subsec:ExperimentalResults}
We now demonstrate the validity of our results experimentally. We consider two regularization functions that satisfy assumption A1: the $\ell_2^2$ regularization and elastic net regularization, where $R(\bm{\theta}) = \lambda_1\norm{\bm{\theta}}_1 + \frac{1}{2}\lambda_2\norm{\bm{\theta}}_2^2$.

We consider standard Gaussian input of dimension $d$ and examine a 2-Layer RF model where both layers are of dimension $p$ and a 1 layer RF model with hidden layer of dimension $p$. The ratio $\frac{n}{d}$ was fixed to 1.5 for all experiments. The activation function was chosen to be $\tanh$.

In figure~\ref{fig:L2reglarization} we show the training and generalization error for $\ell_2^2$ regularization for 3 different regularization values as a function of the ratio $\frac{p}{n}$. We note that in the 1-Layer case $\frac{p}{n}$ is a measure of the under or overparameterization of the network. This relationship does not hold in the two layer case, however as may be seen from the figure this ratio is still useful in comparing the two models. In figure~\ref{fig:L2reglarization} the solid line represents the 2-layer case and the dashed line represents the 1-layer case. The triangles are our theoretical predictions for 2-layers, and squares similarly for 1-layer. For the Elastic net case we fix $\lambda_2$ to be $10^{-5}$ and vary only $\lambda_1$ these results are similarly shown in figure~\ref{fig:elacticNet}.

We note that in both types of regularization functions, for all values of $\frac{p}{n}$, the 2-layer deep RF model has consistently lower generalization error. With respect to training error the two layer case only outperforms 1-layer at large values of regularization. This suggests that even when training of the layer is not performed there can be a benefit to a deeper embedding of the input data. 
\begin{figure}[!hb]
    \centering
    
    \subfigure [Training Error] {
    \label{subfig:trainingErrorL2Synthetic}
    \resizebox{0.45 \textwidth}{!}{%
        \begin{tikzpicture}
        \begin{axis}[
          xlabel={$\gamma = \frac{p}{n}$},
          ylabel=Training Error,
          grid = both,
            minor tick num = 1,
            major grid style = {lightgray},
            minor grid style = {lightgray!25},
            ymode=log
          ]
        \addplot[ smooth, thin, red] table[ y=TrainActivation, x=gamma]{Data/SyntheticL2lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 1$}
        %
        \addplot[ smooth, thin, blue] table[ y=TrainActivation, x=gamma]{Data/SyntheticL2lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 10^{-2}$}
        %
        \addplot[ smooth, thin, green] table[ y=TrainActivation, x=gamma]{Data/SyntheticL2lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 10^{-5}$}
        \addplot[ dashed, thin, red] table[ y=TrainActivation, x=gamma]{Data/SyntheticL21Layerlam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        %
        \addplot[ dashed, thin, blue] table[ y=TrainActivation, x=gamma]{Data/SyntheticL21Layerlam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        %
        \addplot[ dashed, thin, green] table[ y=TrainActivation, x=gamma]{Data/SyntheticL21Layerlam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[red, mark = triangle, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]{Data/SyntheticL2lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
       \addplot[blue, mark = triangle, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]
        {Data/SyntheticL2lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[green, mark = triangle, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]
        {Data/SyntheticL2lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[red, mark = square, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]{Data/SyntheticL21Layerlam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
       \addplot[blue, mark = square, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]
        {Data/SyntheticL21Layerlam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[green, mark = square, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]
        {Data/SyntheticL21Layerlam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \end{axis}
    \end{tikzpicture}
        } 
    }
 \subfigure [Generalization Error] {
    \label{subfig:GeneralizationErrorL2Synthetic}
    \resizebox{0.45\textwidth}{!}{%
        \begin{tikzpicture}
        \begin{axis}[
          xlabel={$\gamma = \frac{p}{n}$},
          ylabel=Generalization Error,
          grid = both,
            minor tick num = 1,
            major grid style = {lightgray},
            minor grid style = {lightgray!25},
            ymode=log
          ]
        \addplot[ smooth, thin, red] table[ y=TestActivation, x=gamma]{Data/SyntheticL2lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 1$}
        %
        \addplot[ smooth, thin, blue] table[ y=TestActivation, x=gamma]{Data/SyntheticL2lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 10^{-2}$}
        %
        \addplot[ smooth, thin, green] table[ y=TestActivation, x=gamma]{Data/SyntheticL2lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 10^{-5}$}
        \addplot[ dashed, thin, red] table[ y=TestActivation, x=gamma]{Data/SyntheticL21Layerlam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        %
        \addplot[ dashed, thin, blue] table[ y=TestActivation, x=gamma]{Data/SyntheticL21Layerlam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        %
        \addplot[ dashed, thin, green] table[ y=TestActivation, x=gamma]{Data/SyntheticL21Layerlam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[red, mark = triangle, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]{Data/SyntheticL2lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
       \addplot[blue, mark = triangle, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]
        {Data/SyntheticL2lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[green, mark = triangle, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]
        {Data/SyntheticL2lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[red, mark = square, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]{Data/SyntheticL21Layerlam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
       \addplot[blue, mark = square, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]
        {Data/SyntheticL21Layerlam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[green, mark = square, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]
        {Data/SyntheticL21Layerlam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \end{axis}
    \end{tikzpicture}
        } 
    }
       \caption{Comparison of 1-Layer and 2-Layer RFs, with square loss function, $\ell_2^2$ regularization with regularization strength $\lambda$. Solid lines represent 2 layer and dashed lines 1-Layer. Triangles are the CGMT results for 2-layers and squares for 1-layer}
    \label{fig:L2reglarization}
\end{figure}


\begin{figure}[!hb]
    \centering
    
    \subfigure [Training Error] {
    \label{subfig:trainingErrorElasticSynthetic}
    \resizebox{0.45 \textwidth}{!}{%
        \begin{tikzpicture}
        \begin{axis}[
          xlabel={$\gamma = \frac{p}{n}$},
          ylabel=Training Error,
          grid = both,
            minor tick num = 1,
            major grid style = {lightgray},
            minor grid style = {lightgray!25},
            ymode=log
          ]
        \addplot[ smooth, thin, red] table[ y=TrainActivation, x=gamma]{Data/SyntheticElasticL21e-5_L1lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 1$}
        %
        \addplot[ smooth, thin, blue] table[ y=TrainActivation, x=gamma]{Data/SyntheticElasticL21e-5_L1lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 10^{-2}$}
        \addplot[ smooth, thin, green] table[ y=TrainActivation, x=gamma]{Data/SyntheticElasticL21e-5_L1lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 10^{-5}$}
        \addplot[ dashed, thin, red] table[ y=TrainActivation, x=gamma]{Data/SyntheticElastic1LayerL21e-5_L1lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[ dashed, thin, blue] table[ y=TrainActivation, x=gamma]{Data/SyntheticElastic1LayerL21e-5_L1lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[ dashed, thin, green] table[ y=TrainActivation, x=gamma]{Data/SyntheticElastic1LayerL21e-5_L1lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[red, mark = triangle, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]{Data/SyntheticElasticL21e-5_L1lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
       \addplot[blue, mark = triangle, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]
        {Data/SyntheticElasticL21e-5_L1lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[green, mark = triangle, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]
        {Data/SyntheticElasticL21e-5_L1lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[red, mark = square, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]{Data/SyntheticElastic1LayerL21e-5_L1lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
       \addplot[blue, mark = square, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]
        {Data/SyntheticElastic1LayerL21e-5_L1lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[green, mark = square, mark size = 1pt, only marks] table[y=TrainGaussian, x=gamma]
        {Data/SyntheticElastic1LayerL21e-5_L1lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \end{axis}
    \end{tikzpicture}
        } 
    }
 \subfigure [Generalization Error] {
    \label{subfig:GeneralizationErrorElasticSynthetic}
    \resizebox{0.45\textwidth}{!}{%
        \begin{tikzpicture}
        \begin{axis}[
          xlabel={$\gamma = \frac{p}{n}$},
          ylabel=Generalization Error,
          grid = both,
            minor tick num = 1,
            major grid style = {lightgray},
            minor grid style = {lightgray!25},
            ymode=log
          ]
        \addplot[ smooth, thin, red] table[ y=TestActivation, x=gamma]{Data/SyntheticElasticL21e-5_L1lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 1$}
        %
        \addplot[ smooth, thin, blue] table[ y=TestActivation, x=gamma]{Data/SyntheticElasticL21e-5_L1lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 10^{-2}$}
        \addplot[ smooth, thin, green] table[ y=TestActivation, x=gamma]{Data/SyntheticElasticL21e-5_L1lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addlegendentry{$\lambda = 10^{-5}$}
        \addplot[ dashed, thin, red] table[ y=TestActivation, x=gamma]{Data/SyntheticElastic1LayerL21e-5_L1lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[ dashed, thin, blue] table[ y=TestActivation, x=gamma]{Data/SyntheticElastic1LayerL21e-5_L1lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[ dashed, thin, green] table[ y=TestActivation, x=gamma]{Data/SyntheticElastic1LayerL21e-5_L1lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[red, mark = triangle, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]{Data/SyntheticElasticL21e-5_L1lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
       \addplot[blue, mark = triangle, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]
        {Data/SyntheticElasticL21e-5_L1lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[green, mark = triangle, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]
        {Data/SyntheticElasticL21e-5_L1lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[red, mark = square, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]{Data/SyntheticElastic1LayerL21e-5_L1lam=1Gamma02-4_70_Delta1.5_100averages_sum300.dat};
       \addplot[blue, mark = square, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]
        {Data/SyntheticElastic1LayerL21e-5_L1lam=1e-2Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \addplot[green, mark = square, mark size = 1pt, only marks] table[y=TestGaussian, x=gamma]
        {Data/SyntheticElastic1LayerL21e-5_L1lam=1e-5Gamma02-4_70_Delta1.5_100averages_sum300.dat};
        \end{axis}
    \end{tikzpicture}
        } 
    }
       \caption{Comparison of 1-Layer and 2-Layer RFs, with square loss function and $\ell_1 + \ell_2^2$ regularization with regularization strength $\lambda$ for the $\ell_1$ term and fixed $\ell_2$ regularization strength. Solid lines represent 2 layer and dashed lines are 1-Layer. Triangles are the CGMT result for 2-layers and squares for 1-layer}
    \label{fig:elacticNet}
\end{figure}



\subsection{Eigendistribution of the Covariance Matrix}
In the CGMT analysis performed above, where the input data is Gaussian, the Gaussian equivalent features $\gamma^{(L)}$ are distributed as $\mathcal{N}(\bm{0}, \mathbf{R}^{(L)})$ where $\mathbf{R}^{(L)}$ is a covariance matrix defined recursively as
\begin{eqnarray}\label{eq:RlyapunovDefinition}
    \mathbf{R}^{(0)} = \mathbf{I}\qquad \mathbf{R}^{(l)} = \rho_{1, l}^2\mathbf{W}^{(l)}\mathbf{R}^{(l-1)}\mathbf{W}^{(l)T} + \rho_{2, l}^2\mathbf{I}_{p_{l}},
\end{eqnarray}
where each $\mathbf{W}^{(l)}$ has rows $\mathbf{w}_j^{(l)}\sim\mathcal{N}(\bm{0}, \frac{1}{p_{l-1}}\mathbf{I}_{p_{l-1}})$. In the case of ridge regression of linear models, or any rotationally invariant setup, the optimal value is directly dependent upon the eigenvalues of the covariance matrix. As the covariance matrix is random we consider its eigendistribution, the marginal probability distribution over the eigenvalues.

We note that the type of recursion for $\mathbf{R}^{(l)}$ is a form of a Lyapanov recursion, which has been studied in the literature \citep{vakili2011RandomMatrix, emery2007lectures}.
We denote the eigendistribution of the matrix $\mathbf{R}^{(l)}$ as $f_{\mathbf{R}^{(l)}}(\lambda)$ for eigenvalues $\lambda$. In the case of $l = 1$, the matrix $\mathbf{R}^{(1)}$ is a scaled Wishart matrix plus an identity, whose eigendistribution is given by a shifted version of the Marchenko–Pastur distribution. In figure~\ref{fig:HistEigDistribution} we consider the empirical eigendistribution of $\mathbf{R}^{(2)}$, corresponding to the two layer case studied above. We choose $p_0 = 1000$ and $p_2= 1500$ fixing the input and output dimensions of the layers, and vary the size of the hidden layer $p_1$. We note as the size of the hidden layer grows the more concentrated the eigendistribution become around zero, while decreasing it results in a more flat structure. In the case of ridge regression, the decreased in the support of the eigenvalues could represent in an increase in model uncertainty at large sizes of the hidden layers.   

We also examine the eigendistribution analytically. We make use of the Stieltjes transform $S_{\mathbf{R}^{(l)}}(z)$ of the distribution $f_{\mathbf{R}^{(l)}}$. This transform and its inverse are give by 
\begin{eqnarray}
    S_{\mathbf{R}^{(l)}}(z) = \int \frac{f_{\mathbf{R}^{(l)}}(\lambda)}{\lambda - z}\mathrm{d}\lambda\qquad f_{\mathbf{R}^{(l)}}(\lambda) = \frac{1}{\pi}\lim_{\omega \rightarrow 0^+}\mathrm{Im}[S(\lambda + i \omega)]
\end{eqnarray}
where $i$ is the imaginary unit, and $z$ is complex. We can demonstrate that the Stieltjes transform of the matrices $\mathbf{R}^{(l)}$ follows the following recursion.
\begin{theorem}
    Let $\beta_l = \frac{p_l}{p_{l-1}}$, then the Stieltjes transform $S_l(z)$ of $\mathbf{R}^{(l)}$ in \eqref{eq:RlyapunovDefinition} is given recursively by
    \begin{eqnarray}
        S_{l+1}(z) = \frac{1}{\rho_{1, l+1}^2}\Omega_{l}\left(\frac{z - \rho_{2, l+1}^2}{\rho_{1,l+1}^2} \right)\\
        \Omega_{l}(z) = \frac{1}{1 - \beta - \beta z\Omega_l(z)}S_l\left(\frac{z}{1 - \beta - \beta z\Omega_l(z)}\right)
    \end{eqnarray}
    Where $\Omega_0$ is the Stieltjes transform of a Wishart matrix, given by
    \begin{eqnarray}
        \Omega_0 = \frac{1-\beta_1 - z + \sqrt{z^2 - 2(\beta_1 +1)z + (\beta_1 -1)^2}}{2\beta_1z}
    \end{eqnarray}
\end{theorem}
\begin{proof}
The proof is given in appendix~\ref{app:sec:LyapanovProofs}.
\end{proof}
The recursive definitions given are difficult to compute empirically, as such we will leave visualizing these results to future work. However the recursive structure suggests that there exists a limiting distribution over the eigenvalues in the limit of infinite depth characterized by the different ratio in size between the various layers. 


 \begin{figure}[!hb]
    \centering
    \label{fig:HistEigDistribution}
    \subfigure [$p_1 \leq p_0$] {
    \label{subfig:StieltjesTransforms1}
    \resizebox{0.45 \textwidth}{!}{%
        \begin{tikzpicture}
        \begin{axis}[
          xlabel={$\lambda$},
          ylabel=Probability,
          grid = both,
            minor tick num = 1,
            major grid style = {lightgray},
            minor grid style = {lightgray!25},
            bar width = 0.01320055239333573,
           
            ymax =2
          ]
        \addplot[ybar,fill, black, opacity=0.33] table [x, y, col sep=space] {StieltjesData/HistVarmiddle=100.dat};
                \addlegendentry{$p_1= 100$}
        \addplot[ybar,fill ,blue, opacity=0.33] table [x, y, col sep=space] {StieltjesData/HistVarmiddle=500.dat};
                \addlegendentry{$p_1= 500$}
        \addplot[ybar,fill, red, opacity=0.33] table [x, y, col sep=space] {StieltjesData/HistVarmiddle=1000.dat};
                \addlegendentry{$p_1= 1000$}
       
       
        \end{axis}
        \end{tikzpicture}
         } 
        }
    \subfigure [$p_1 > p_0$] {
    \label{subfig:StieltjesTransforms2}
    \resizebox{0.45 \textwidth}{!}{%
        \begin{tikzpicture}
        \begin{axis}[
          xlabel={$\lambda$},
          ylabel=Probability,
          grid = both,
            minor tick num = 1,
            major grid style = {lightgray},
            minor grid style = {lightgray!25},
            bar width = 0.01320055239333573,
           
            ymax =5
          ]
        \addplot[ybar,fill, black, opacity=0.33] table [x, y, col sep=space] {StieltjesData/HistVarmiddle=2000.dat};
        \addlegendentry{$p_1= 2000$}
        \addplot[ybar,fill ,blue, opacity=0.33] table [x, y, col sep=space] {StieltjesData/HistVarmiddle=5000.dat};
                \addlegendentry{$p_1= 5000$}
        \addplot[ybar,fill, red, opacity=0.33] table [x, y, col sep=space] {StieltjesData/HistVarmiddle=10000.dat};
                \addlegendentry{$p_1= 10000$}
       
       
        \end{axis}
    \end{tikzpicture}
        } 
    }
    \caption{Empirical Eigendistribution of $\mathbf{R}^{(l)}$ for various sizes $p_1$ of the 1st hidden layer}
\end{figure}







\section{Conclusion}\label{sec:Conclusion}
In this paper, we prove an asymptotic equivalence between deep random feature models and linear Gaussian models with respect to the training and generalization error. As a result of this universality, we can study a Gaussian equivalent model to the DRF model, in the asymptotic limit. We use this fact to provide an exact asymptotic analysis by means of the convex Gaussian min max theorem for an $L$-layer deep random feature model with Gaussian inputs. We further demonstrate that depth has an effect on training and generalization error both experimentally and by studying the eigendistribution of the Gaussian equivalent model's Covariance matrix. 



\bibliographystyle{abbrv}

