Random Features Outperform Linear Models: Effect of Strong Input-Label Correlation in Spiked Covariance Data

## Abstract

Abstract Random Feature Model (RFM) with a nonlinear activation function is instrumental in understanding training and generalization performance in high-dimensional learning. While existing research has established an asymptotic equivalence in performance between the RFM and noisy linear models under isotropic data assumptions, empirical observations indicate that the RFM frequently surpasses linear models in practical applications. To address this gap, we ask, “When and how does the RFM outperform linear models?” In practice, inputs often have additional structures that significantly influence learning. Therefore, we explore the RFM under anisotropic input data characterized by spiked covariance in the proportional asymptotic limit, where dimensions diverge jointly while maintaining finite ratios. Our analysis reveals that a high correlation between inputs and labels is a critical factor enabling the RFM to outperform linear models. Moreover, we show that the RFM performs equivalent to noisy polynomial models, where the polynomial degree depends on the strength of the correlation between inputs and labels. Our numerical simulations validate these theoretical insights, confirming the performance-wise superiority of RFM in scenarios characterized by strong input-label correlation.

## 1 Introduction

Random Feature Model (RFM) (Rahimi & Recht, ) has received significant attention in recent years due to its compelling theoretical properties and its connections to neural networks. In this work, we focus on a typical supervised learning framework utilizing the RFM, expressed as

$$ ${\bm{\omega}}^{T}\sigma({\mathbf{F}}{\mathbf{x}}),$ (1) $$

where ${\mathbf{x}}\in\mathbb{R}^{n}$ is an input vector, ${\mathbf{F}}\in\mathbb{R}^{k\times n}$ is called the random feature matrix, $\sigma:\mathbb{R}\to\mathbb{R}$ is an activation function, and ${\bm{\omega}}\in\mathbb{R}^{k}$ is a vector of learnable parameters. The feature matrix ${\mathbf{F}}$ is randomly sampled and remains fixed, allowing the RFM to be interpreted as a two-layer neural network with fixed weights in the first layer. The learning process involves optimizing the following objective function based on a set of training samples $\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{m}$

$$ $\hat{{\bm{\omega}}}_{\sigma}:=\underset{{\bm{\omega}}\in\mathbb{R}^{k}}{\mbox{ argmin}}\frac{1}{m}\sum_{i=1}^{m}(y_{i}-{\bm{\omega}}^{T}\sigma({\mathbf{F}}{ \mathbf{x}}_{i}))^{2}+\lambda||{\bm{\omega}}||_{2}^{2},$ (2) $$

where $\lambda\in\mathbb{R}^{+}$ is a regularization constant. While the RFM can be studied under various loss functions, we concentrate on squared loss for simplicity in our theoretical development. The training error of the RFM with activation function $\sigma$ is defined as

$$ ${\mathcal{T}}_{\sigma}:=\frac{1}{m}\sum_{i=1}^{m}(y_{i}-\hat{{\bm{\omega}}}_{ \sigma}^{T}\sigma({\mathbf{F}}{\mathbf{x}}_{i}))^{2}+\lambda||\hat{{\bm{\omega }}}_{\sigma}||_{2}^{2}.$ (3) $$

Similarly, one can measure the performance of a trained RFM via the generalization error defined as

$$ ${\mathcal{G}}_{\sigma}:=\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}(y-\hat{{\bm{ \omega}}}_{\sigma}^{T}\sigma({\mathbf{F}}{\mathbf{x}}))^{2}.$ (4) $$

We study the training and generalization errors of RFM (1) in the “proportional asymptotic limit,” where the number of training samples $m$, the input dimension $n$, and the number of intermediate features $k$ jointly diverge: $m,n,k\to\infty$ with $n/m,n/k\in(0,\infty)$. Under an isotropic data assumption, the RFM has been shown to perform equivalent to the following noisy linear model (Montanari et al., ; Gerace et al., ; Goldt et al., ; Mei & Montanari, ; Hu & Lu, ):

$$ $\displaystyle{\bm{\omega}}^{T}(\mu_{0}{\mathbf{1}}+\mu_{1}{\mathbf{F}}{\mathbf {x}}+\mu_{*}{\mathbf{z}}),$ (5) $$

where ${\mathbf{1}}$ is an all-one vector, ${\mathbf{z}}\sim{\mathcal{N}}(0,{\mathbf{I}}_{k})$, and $\mu_{0},\mu_{1},\mu_{*}$ are constants depending on the activation function $\sigma$ as $\mu_{0}=\mathop{{\mathbb{E}}}[\sigma(z)]$, $\mu_{1}=\mathop{{\mathbb{E}}}[z\sigma(z)]$, and $\mu_{*}=\sqrt{{\mathop{{\mathbb{E}}}[\sigma(z)^{2}]-\mu_{1}^{2}-\mu_{0}^{2}}}$ where $z\sim\mathcal{N}(0,1)$. The noisy linear model enables the performance characterization of the RFM (Dhifallah & Lu, ). It is also useful for studying double-descent phenomenon (Belkin et al., ; Geiger et al., ) with regards to loss functions (Mei & Montanari, ), activation functions (Wang & Bento, ; Demir & Doğan, ), and regularization constant (Nakkiran et al., ). However, such a noisy linear model falls short in explaining the empirical performance of the RFM on real-world data (Ghorbani et al., ). This motivates us to investigate the following question:

“Under which conditions does the random feature model outperform linear models?”

The equivalence of the RFM to the noisy linear model is closely tied to the isotropic data assumption, i.e., ${\mathbf{x}}\sim{\mathcal{N}}(0,{\mathbf{I}}_{n})$. However, in practice, inputs often exhibit additional structures. To address this limitation, we consider spiked covariance data model (Johnstone, ; Baik et al., )

$$ $\displaystyle{\mathbf{x}}$ $\displaystyle\sim{\mathcal{N}}(0,{\mathbf{I}}_{n}+\theta{\bm{\gamma}}{\bm{ \gamma}}^{T}),\quad y:=\sigma_{*}\left(\frac{{\bm{\xi}}^{T}{\mathbf{x}}}{\sqrt {1+\theta\alpha^{2}}}\right),$ (6) $$

which introduces anisotropic characteristics due to the structured covariance of ${\mathbf{x}}$. Here, ${\bm{\gamma}}\in\mathbb{R}^{n}$ serves as a fixed but unknown “spike signal” with $||{\bm{\gamma}}||_{2}=1$, ${\bm{\xi}}\in\mathbb{R}^{n}$ represents a fixed but unknown “label signal” also normalized to one, i.e., $||{\bm{\xi}}||_{2}=1$. The parameter $\theta$, termed “spike magnitude” scales with the dimension as $\theta\asymp n^{\beta}$ where $\beta\in[0,1/2)$, and $\sigma_{*}:\mathbb{R}\to\mathbb{R}$ is a nonlinear target function. Note that the input of the target function $\sigma_{*}$ is scaled to have unit variance. Additionally, we define the alignment parameter $\alpha:={\bm{\gamma}}^{T}{\bm{\xi}}$, which governs the input-label correlation.

Under the spiked data model, we first establish a “universality” theorem indicating that the RFM exhibits equivalent performance across two different activation functions if the first two statistical moments of $(\sigma({\mathbf{F}}{\mathbf{x}}),y)$ are consistent for both activations. Leveraging this theorem, we extend our analysis to demonstrate that high-order polynomial models perform equivalent to the RFM beyond traditional linear frameworks (5). Notably, we show that the equivalence to the noisy linear model still holds provided that at least one of the following is sufficiently small: (i) the cosine similarity between the rows of ${\mathbf{F}}$ and the ${\bm{\gamma}}$, (ii) the spike magnitude $\theta$, and (iii) the alignment parameter $\alpha$. A precise characterization of these conditions will be detailed in Section 4.
Moreover, we introduce a new equivalence between the RFM and noisy polynomial models, generalizing previous results on noisy linear models (5). Finally, this framework allows us to compare different activation functions within the context of spiked data, highlighting how nonlinearity can enhance learning performance beyond linear models.

Overall, our contributions are as follows:

- 1.
We show the conditions under which the RFM outperforms linear models under anisotropic data settings, which indicates a key role played by the input-label correlation in data.
- 2.
We establish a theorem that shows high-order polynomial models perform equivalent to the RFM under spiked covariance data assumption.
- 3.
We relax the isotropic data assumption previously utilized by Hu & Lu () for the “universality of random features,” which may be of particular interest on its own.

By addressing these aspects, our work not only clarifies when and how RFMs can excel in the spiked covariance model setting but also enriches the theoretical understanding of their performance relative to linear models.

## 2 Related Work

The RFM has emerged as a vital tool in the realm of machine learning, initially introduced as a randomized approximation to kernel methods (Rahimi & Recht, ). Its significance has grown in parallel with the development of neural networks, particularly in how both the RFM and the Neural Tangent Kernel (NTK) offer distinct linear approximations of two-layer neural networks (Ghorbani et al., ). While empirical evidence shows that neural networks often outperform both the RFM and NTK (Ghorbani et al., ), analyzing the RFM provides critical insights into the behavior of nonlinear models, serving as a bridge toward understanding more complex architectures.

A substantial body of research has focused on the equivalence between the RFM and noisy linear models with regards to their performances,(Montanari et al., ; Gerace et al., ; Goldt et al., ; ; Dhifallah & Lu, ; Mei & Montanari, ). For instance, Hu & Lu () demonstrated this equivalence in terms of training and generalization performance using Lindeberg’s method. Additionally, Montanari & Saeed () showed that nonlinear features could be substituted with equivalent Gaussian features in empirical risk minimization settings. However, these studies mainly examined the RFM under isotropic data assumptions, which limits their applicability to real-world scenarios characterized by more complex data structures. Our work extends these findings by investigating the RFM under anisotropic data conditions, thereby addressing a significant gap in the literature. Recent works have indicated that equivalent Gaussian features may not consistently provide equivalent training and generalization errors for nonlinear models (Gerace et al., ; Pesce et al., ). This discrepancy has led researchers to explore alternative feature representations, including Gaussian mixtures (Dandi et al., ). However, our study diverges from this line of inquiry by concentrating on identifying equivalent activation functions rather than seeking alternatives to Gaussian features.

In light of the RFM’s limitations in surpassing linear models, recent research has turned to neural networks trained with one gradient step as a potential solution to bridge this performance gap (Ba et al., ; Damian et al., ). Studies by Ba et al. () and Mousavi-Hosseini et al. () have explored one-step gradient descent techniques in relation to sample complexity under spiked covariance data assumption. Furthermore, Moniri et al. () explored neural networks trained with one gradient step in an isotropic data setting and showed that the update can be characterized by the appearance of a spike in the spectrum of the feature matrix. Dandi et al. (); Cui et al. () studied the one-step gradient descent with a higher learning rate in comparison to Moniri et al. (). Ba et al. (); Moniri et al. (); Dandi et al. (); Cui et al. () also stated Gaussian equivalence results in terms of training and generalization errors when studying the one gradient step technique under isotropic Gaussian data. While these works primarily focus on surpassing the performance of linear models through feature learning in neural networks, our research aims to clarify when and how the RFM can outperform linear models without relying on feature learning.

Additionally, there is an ongoing exploration into deep random features that integrate multiple layers of random projections interleaved with nonlinearities (Schröder et al., ; Bosch et al., ). Research by Zavatone-Veth & Pehlevan () has also examined deep random features with anisotropic weight matrices using statistical physics techniques. Moreover, studies have investigated the RFM under polynomial scaling limits where $k,n,m\to\infty$ while maintaining specific ratios between these parameters, i.e., $k/n^{c_{1}},m/k^{c_{2}}\in(0,\infty)$ for some constants $c_{1},c_{2}>0$ (Hu et al., ), further providing equivalence results for the RFM under the polynomial scaling limit.

## 3 Preliminaries

##### Notations

We denote vectors and matrices using bold letters, with lowercase letters representing vectors and uppercase letters representing matrices. Scalars are indicated by non-bold letters. The notation $||{\mathbf{x}}||_{2}$ refers to the Euclidean norm of the vector ${\mathbf{x}}$ while $||{\mathbf{F}}||$ denotes to spectral norm of the matrix ${\mathbf{F}}$. We use the term $\text{polylog }k$ to represent any polylogarithmic function of $k$. The symbol $\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}$ indicates convergence in probability. We employ big-O notation, denoted $\mathcal{O}(.)$, and small-o notation, represented as $o(.)$ with respect to the parameters $k,n,m$. These notations are

$$ $\displaystyle h(k)=\mathcal{O}(g(k))\quad\Leftrightarrow\quad\limsup_{k\to \infty}\frac{h(k)}{g(k)}<\infty,$ (7) $\displaystyle h(k)=o(g(k))\quad\Leftrightarrow\quad\lim_{k\to\infty}\frac{h(k) }{g(k)}=0.$ (8) $$

It is important to note that $h(k)=o(1/\text{polylog }k)$ implies that $\lim_{k\to\infty}h(k)\text{ polylog }k=0$ for all polylogarithmic functions. We also use the notation $\theta\asymp n^{\beta}$ to indicate that there exists some constants $C_{1},C_{2}>0$ such that $C_{1}n^{\beta}\leq\theta\leq C_{2}n^{\beta}$. Finally, we denote the $i$-th Hadamard power of a vector ${\mathbf{x}}$ as ${\mathbf{x}}^{\circ i}$.

##### Assumptions

We establish our results based on the following assumptions.

- (A.1)
The spike and label signal vectors ${\bm{\gamma}},{\bm{\xi}}\in\mathbb{R}^{n}$ are deterministic and $||{\bm{\gamma}}||_{2}=||{\bm{\xi}}||_{2}=1$.
- (A.2)
${\mathbf{x}}\sim{\mathcal{N}}(0,{\mathbf{I}}_{n}+\theta{\bm{\gamma}}{\bm{
\gamma}}^{T})$ for spike magnitude $\theta\asymp n^{\beta}$ where $\beta\in[0,1/2)$.
- (A.3)
$y:=\sigma_{*}\left(\frac{{\bm{\xi}}^{T}{\mathbf{x}}}{\sqrt{1+\theta\alpha^{2}}
}\right)$ where $\sigma_{*}:\mathbb{R}\to\mathbb{R}$ is a function satisfying $|\sigma_{*}(x)|<C(1+|x|^{K})$ for all $x\in\mathbb{R}$ and some constants $C>0$, $K\in\mathbb{Z}^{+}$.
- (A.4)
The number of training samples $m$, dimension of input vector $n$ and number of intermediate features $k$ jointly diverge: $m,n,k\to\infty$ with $n/m,n/k\in(0,\infty)$.
- (A.5)
${\mathbf{F}}:=[{\mathbf{f}}_{1},{\mathbf{f}}_{2},\dots,{\mathbf{f}}_{k}]^{T}$ where ${\mathbf{f}}_{i}\sim{\mathcal{N}}(0,\frac{1}{n+\theta}{\mathbf{I}}_{n})$. Note that the covariance of ${\mathbf{f}}_{i}$ is selected such that $\mathop{{\mathbb{E}}}_{{\mathbf{x}},{\mathbf{f}}_{i}}[({\mathbf{f}}_{i}^{T}{
\mathbf{x}})^{2}]=1$.
- (A.6)
The activation function $\sigma:\mathbb{R}\to\mathbb{R}$ is an odd function satisfying $|\sigma(x)|<C(1+|x|^{K})$ for all $x\in\mathbb{R}$ and some constants $C>0$, $K\in\mathbb{Z}^{+}$.

##### Discussion of Assumptions

The assumptions outlined in (A.1)-(A.3) establish the foundational framework for our data model and setup. Specifically, (A.1) describes the general characteristics of the spike and label signals, while (A.2) introduces the spike magnitude $\theta$. Although there may be potential to extend the range of $\theta$ in future work, our current proofs necessitate that $\beta<1/2$. This restriction is crucial for ensuring that the mathematical properties we derive remain valid within the specified limits. In (A.3), we scale the input of the target function $\sigma_{*}$ to get unit variance. This scaling not only simplifies our proofs but also aligns with standard practices in theoretical analysis. Assumption (A.4) specifies the proportional asymptotic limit, which is a critical aspect of our study. By requiring that the dimensions jointly diverge while maintaining finite ratios, we can effectively analyze the behavior of the RFM in high-dimensional settings. In (A.5), we detail the distribution of the feature matrix ${\mathbf{F}}$, emphasizing that the $1/(n+\theta)$ scaling is vital for our results to hold. This specific scaling ensures that the covariance structure of ${\mathbf{F}}$ is appropriately calibrated with respect to both the input data and the spike magnitude, facilitating accurate performance characterization. Assumption (A.6) pertains to the properties of the activation function $\sigma$, which warrants further elaboration. Previous research by Hu & Lu () established a performance-wise equivalence between the RFM and noisy linear models under isotropic data conditions, relying on the premise that $\sigma$ is an odd function with bounded first, second, and third derivatives. In our case, (A.6) implies the derivatives of $\sigma(x)$ bounded by $\text{polylog }k$ if $|x|\leq\text{polylog }k$, which happens with high probability for $x\sim{\mathcal{N}}(0,1)$. Finally, it is noteworthy that while ReLU (9) does not conform to the odd function assumption stipulated in (A.6), empirical evidence suggests that our findings remain valid even when using ReLU as an activation function.

##### Activation Functions for Numerical Results

In our numerical simulations, we employ three widely used activation functions: ReLU (Rectified Linear Unit), tanh (hyperbolic tangent), and Softplus. The mathematical definitions of these activation functions are as follows:

$$ $\displaystyle\sigma_{ReLU}(x):=\text{max}(0,x),\quad\sigma_{tanh}(x):=\frac{e^ {2x}-1}{e^{2x}+1},\quad\sigma_{Softplus}:=\text{log}(1+e^{x}).$ (9) $$

## 4 Main Results

In this section, we describe our main theoretical results, significantly enhancing the understanding of RFM under structured data conditions. In Section 4.1, we introduce a new “universality” theorem that establishes the performance-wise equivalence of the RFM when utilizing two distinct activation functions, $\sigma$ and $\hat{\sigma}$. This theorem asserts that if the first two statistical moments of the pairs $(\sigma({\mathbf{F}}{\mathbf{x}}),y)$ and $(\hat{\sigma}({\mathbf{F}}{\mathbf{x}}),y)$ are matched, the training and generalization performance of the RFM will remain equivalent across both functions. Building on this foundation, Section 4.2 explores the scenario where $\hat{\sigma}$ is represented as a finite-degree Hermite expansion of $\sigma$, demonstrating that the RFM retains equivalent performance regardless of the chosen activation function. Notably, we reveal that the degree of this expansion is intricately tied to the input-label correlation, offering insights into how data structure influences model performance. In Section 4.3, we establish conditions for the equivalence between the RFM and a noisy linear model, providing a crucial corollary to our earlier findings. When these conditions are not satisfied, we identify that a high-order polynomial model serves as an equivalent to the RFM, a topic further elaborated in Section 4.4.

### 4.1 Universality of Random Features under Spiked Data

Here, we extend the universality laws established by Hu & Lu () for RFM to encompass spiked covariance data. The concept of “universality of random features” implies that one can replace $\sigma({\mathbf{F}}{\mathbf{x}})$ in the problem with an equivalent Gaussian vector, as articulated in equation (5), without compromising the training and generalization performance of the model. While Hu & Lu () demonstrated this result under mild assumptions regarding the feature matrix ${\mathbf{F}}$, the nonlinearity $\sigma$, and the loss function, their findings were constrained by an isotropic data assumption. To address these limitations, we propose a new universality theorem tailored for scenarios involving spiked covariance data, as outlined in equation (6).

###### Theorem 1 .

Let $\sigma(x)$, $\hat{\sigma}(x)$ be two activation functions satisfying (A.6). Suppose that assumptions (A.1)-(A.6) given in Section 3 hold. If

$$ $\displaystyle\left\|\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}[\sigma({\mathbf{F }}{\mathbf{x}})\sigma({\mathbf{F}}{\mathbf{x}})^{T}]-\mathop{{\mathbb{E}}}_{({ \mathbf{x}},y)}[\hat{\sigma}({\mathbf{F}}{\mathbf{x}})\hat{\sigma}({\mathbf{F} }{\mathbf{x}})^{T}]\right\|=o\left(1/\text{polylog }k\right),$ (10) $\displaystyle\left\|\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left[\sigma({ \mathbf{F}}{\mathbf{x}})y\right]-\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left [\hat{\sigma}({\mathbf{F}}{\mathbf{x}})y\right]\right\|=o\left(1/\text{polylog }k\right),$ (11) $$

are satisfied, then

- (i)
the training error ${\mathcal{T}}_{\sigma}$ of the RFM with activation $\sigma$, and the training error ${\mathcal{T}}_{\hat{\sigma}}$ of the RFM with activation $\hat{\sigma}$, both converge in probability to the same value $e_{{\mathcal{T}}}\geq 0$,
- (ii)
the corresponding generalization errors ${\mathcal{G}}_{\sigma}$ and ${\mathcal{G}}_{\hat{\sigma}}$ also converge in probability to the same value $e_{{\mathcal{G}}}\geq 0$ under additional assumptions (A.7) and (A.8) provided in Appendix A.

###### Proof.

We provide a proof in Appendix A.

###### Theorem 1 .

###### Proof.

##### Proof Approach

We follow the proof technique used by Hu & Lu (). First, we define a perturbed training objective that includes additional terms related to the generalization error. Then, we bound the expected difference in the perturbed training objective resulting from replacing $\sigma({\mathbf{F}}{\mathbf{x}}_{i})$ with $\hat{\sigma}({\mathbf{F}}{\mathbf{x}}_{i})$ for each training sample ${\mathbf{x}}_{i}$, following the principles of Lindeberg’s method (Lindeberg, ; Korada & Montanari, ).
∎

A key advancement of our work is encapsulated in Theorem 4.1, which establishes the equivalence of the RFM across two different activation functions. This result is notably more general than previous findings that focused solely on the equivalence of specific models. Specifically, Theorem 4.1 asserts that substituting the activation function $\sigma$ with $\hat{\sigma}$ does not impact the training and generalization performance of the RFM, provided that conditions (10) and (11) are satisfied. In the subsequent sections, we will leverage this theorem to identify models that are equivalent to the RFM and outline the conditions necessary for these equivalences to hold, thereby broadening the applicability of our findings in various contexts.

### 4.2 Noisy Polynomial Models Equivalent to the Random Feature Model

Now, we consider generalizing the equivalent noisy linear model to equivalent high-order polynomial models by utilizing Theorem 4.1. To achieve this, we leverage Hermite polynomials (O’Donnell, , Chapter 11.2). For $i\in\mathbb{Z}^{+}$, $i$-th Hermite polynomial is defined as

$$ $\displaystyle H_{i}(x):=(-1)^{i}e^{x^{2}/2}\frac{d^{i}}{dx^{i}}e^{-x^{2}/2}$ (12) $$

for any $x\in\mathbb{R}$. There are two important properties of Hermite polynomials. First, they are orthogonal with respect to Gaussian measure, so $\mathop{{\mathbb{E}}}_{z\sim{\mathcal{N}}(0,1)}[H_{i}(z)H_{j}(z)]=0$ for $i\neq j$. Second, they form an orthogonal basis of the Hilbert space of functions $\sigma:\mathbb{R}\to\mathbb{R}$ satisfying $\mathop{{\mathbb{E}}}_{z\sim{\mathcal{N}}(0,1)}[\sigma(z)^{2}]<\infty$. This means that the function $\sigma(x)$ can be expressed as an infinite sum of Hermite polynomials as

$$ $\displaystyle\sigma(x)=\sum_{i=0}^{\infty}\frac{1}{i!}\mu_{i}H_{i}(x),\quad \text{ with }\quad\mu_{i}:=\mathop{{\mathbb{E}}}_{z\sim{\mathcal{N}}(0,1)}[H_{ i}(z)\sigma(z)].$ (13) $$

This sum is also known as “Hermite expansion” and $\mu_{i}$ is called $i$-th Hermite coefficient. Since we are working with Gaussian random variables, Hermite expansion simplifies the analysis. Indeed, we can construct a performance-wise equivalent activation function using the Hermite expansion. Specifically, the following theorem states that the RFM with activation function $\sigma(x)$ performs equivalent to the RFM with $\hat{\sigma}_{l}(x)$, which is a finite-degree Hermite expansion with additional noise.

###### Theorem 2 .

Let $\sigma$ be any activation function satisfying (A.6). Define another activation function

$$ $\displaystyle\hat{\sigma}_{l}(x):=\left(\sum_{j=0}^{l-1}\frac{1}{j!}\mu_{j}H_{ j}(x)\right)+\mu^{*}_{l}z\quad\text{with}\quad z\sim{\mathcal{N}}(0,1),$ (14) $$

where $\mu_{j}:=\mathop{{\mathbb{E}}}_{z\sim{\mathcal{N}}(0,1)}[H_{j}(z)\sigma(z)]$ and $\mu^{*}_{l}:=\sqrt{\mathop{{\mathbb{E}}}_{z\sim{\mathcal{N}}(0,1)}[\sigma(z)^{
2}]-\sum_{j=0}^{l-1}\mu_{j}^{2}/(j!)}$.

Suppose that assumptions (A.1)-(A.6) given in Section 3 hold. If

$$ $\eta:=\max_{1\leq i\leq k}\left|\frac{({\bm{\xi}}+\theta\alpha{\bm{\gamma}})^{ T}{\mathbf{f}}_{i}}{\sqrt{1+\theta\alpha^{2}}}\right|\leq\frac{C}{n^{1/l}}, \quad\text{for some $C>0$ and some $l\in\mathbb{Z}^{+}$,}$ (15) $$

is satisfied, then

- (i)
the training error ${\mathcal{T}}_{\sigma}$ of the RFM with activation $\sigma$, and the training error ${\mathcal{T}}_{\hat{\sigma}_{l}}$ of the RFM with activation $\hat{\sigma}_{l}$, both converge in probability to the same value $e_{{\mathcal{T}}}\geq 0$,
- (ii)
the corresponding generalization errors ${\mathcal{G}}_{\sigma}$ and ${\mathcal{G}}_{\hat{\sigma}_{l}}$ also converge in probability to the same value $e_{{\mathcal{G}}}\geq 0$ under additional assumptions (A.7) and (A.8) provided in Appendix A.

###### Proof.

In Appendix C, we show

$$ $\displaystyle\left\|\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}[\sigma({\mathbf{F }}{\mathbf{x}})\sigma({\mathbf{F}}{\mathbf{x}})^{T}]-(\mu_{1}^{2}{\mathbf{F}}( {\mathbf{I}}_{n}+\theta{\bm{\gamma}}{\bm{\gamma}}^{T}){\mathbf{F}}^{T}+\mu_{*} ^{2}{\mathbf{I}}_{k})\right\|=o\left(1/\text{polylog }k\right),$ (16) $$

for any $\sigma$ satisfying (A.6), $\mu_{1}=\mathop{{\mathbb{E}}}[z\sigma(z)]$ and $\mu_{*}=\sqrt{\mathop{{\mathbb{E}}}[\sigma(z)^{2}]-\mu_{1}^{2}}$ with $z\sim{\mathcal{N}}(0,1)$. By definition, we have $\mathop{{\mathbb{E}}}[z\sigma(z)]=\mathop{{\mathbb{E}}}[z\hat{\sigma}_{l}(z)]$. One can easily show $\mathop{{\mathbb{E}}}[\sigma(z)^{2}]=\mathop{{\mathbb{E}}}[\hat{\sigma}_{l}(z)
^{2}]$. Therefore, (16) holds for $\sigma$ and $\hat{\sigma}_{l}$ with same $\mu_{1}$ and $\mu_{*}$ values.
Using (16) and triangle inequality, we get

$$ $\displaystyle\left\|\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}[\sigma({\mathbf{F }}{\mathbf{x}})\sigma({\mathbf{F}}{\mathbf{x}})^{T}]-\mathop{{\mathbb{E}}}_{({ \mathbf{x}},y)}[\hat{\sigma}_{l}({\mathbf{F}}{\mathbf{x}})\hat{\sigma}_{l}({ \mathbf{F}}{\mathbf{x}})^{T}]\right\|=o\left(1/\text{polylog }k\right).$ (17) $$

Similarly, in Appendix D, we prove

$$ $\left\|\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left[\sigma({\mathbf{F}}{ \mathbf{x}})y\right]-\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left[\hat{\sigma }_{l}({\mathbf{F}}{\mathbf{x}})y\right]\right\|=o\left(1/\text{polylog }k \right),$ (18) $$

when $|\eta|=\mathcal{O}(n^{-1/l})$ for some $l\in\mathbb{Z}^{+}$ and all $i\in\{1,\dots,k\}$. Finally, using Theorem 4.1 with (17) and (18), we reach the statement of the theorem.
∎

Theorem 2 specifies a finite-degree polynomial function (14) that is equivalent to any given activation function $\sigma(x)$. Furthermore, the degree of the polynomial is related to $\eta$ defined in (15). It is important to note that $\eta$ is related to the input-label correlation, which can be decomposed as

$$ $\displaystyle\mathop{{\mathbb{E}}}[({\mathbf{x}}-\mathop{{\mathbb{E}}}[{ \mathbf{x}}])(y-\mathop{{\mathbb{E}}}[y])]=\tilde{\mu}_{1}\frac{{\bm{\xi}}+ \theta\alpha\gamma}{\sqrt{1+\theta\alpha^{2}}}+\sum_{j>1}\tilde{\mu}_{j} \mathop{{\mathbb{E}}}_{{\mathbf{x}}}[{\mathbf{x}}H_{j}({\bm{\xi}}^{T}{\mathbf{ x}})],$ (19) $$

where $\tilde{\mu}_{j}:=\mathop{{\mathbb{E}}}_{z\sim{\mathcal{N}}(0,1)}[\sigma_{*}(z)
H_{j}(z)]$ denotes $j$-th Hermite coefficient of $\sigma_{*}$. Here, $\eta$ depends on $({\bm{\xi}}+\theta\alpha\gamma)/\sqrt{1+\theta\alpha^{2}}$, which is the first term in the decomposition of the input-label correlation (19). Additionally, a higher value of $\|{\bm{\xi}}+\theta\alpha\gamma\|/\sqrt{1+\theta\alpha^{2}}$ implies a greater value of $\eta$ for any ${\mathbf{F}}$ sampled independent of ${\bm{\gamma}}$ and ${\bm{\xi}}$.

In the following sections, we will first examine the case where the equivalent activation function in equation (14) is linear, effectively reducing the model to the noisy linear model described in equation (5). Subsequently, we will focus on the more general case where the equivalent activation function is represented as a high-order polynomial.

###### Theorem 2 .

###### Proof.

### 4.3 Condition of the Equivalence to the Noisy Linear Model

In this section, we explore the “linear regime” of the RFM and delineate the conditions under which it aligns with the noisy linear model in the context of spiked data. Specifically, we examine the substitution of the activation function $\sigma(x)$ with $\hat{\sigma}_{1}(x):=(\mu_{0}+\mu_{1}x+\mu_{*}z)$ for $z\sim{\mathcal{N}}(0,1)$. This transformation effectively simplifies the model to that of a noisy linear model, as described in equation (5). The conditions necessary for establishing equivalence between the noisy linear model and the RFM are articulated in the following corollary, which emerges directly from Theorem 2.

Figure: (a) Percentage difference between the generalization errors of the RFM and the noisy linear model (5). The red line is predicted by Corollary 3.
Refer to caption: extracted/5889816/figures/alignment.png

###### Corollary 3 .

The noisy linear model (5) performs equivalent to the RFM with activation $\sigma$ in terms of training and generalization errors if $\eta=\mathcal{O}(n^{-1/2})$ in (15).

Figure 1(a) vividly illustrates the disparity in generalization errors between the RFM and the noisy linear model across various degrees of alignment $\alpha$ and different values of spike magnitude $\theta$. The red line delineates the boundary of equivalence, as established by Corollary 3. Given that ${\mathbf{f}}_{i}$ is sampled independently of ${\bm{\gamma}}$ and ${\bm{\xi}}$ , we observe that both $|{\bm{\gamma}}^{T}{\mathbf{f}}_{i}|$ and $|{\bm{\xi}}^{T}{\mathbf{f}}_{i}|$ scale as $\mathcal{O}(n^{-1/2})$. Consequently, the expression $|\theta\alpha/\sqrt{1+\theta\alpha^{2}}|=\mathcal{O}(1)$ leads to the conclusion that $\eta=\mathcal{O}(n^{-1/2})$, which fulfills the condition outlined in Corollary 3.

As an illustrative example, we examine a misaligned scenario where $\alpha:={\bm{\xi}}^{T}{\bm{\gamma}}=\mathcal{O}(n^{-1/2})$. This condition is likely to occur when both ${\bm{\gamma}}$ and ${\bm{\xi}}$ are independently and identically sampled from ${\mathcal{N}}(0,{\mathbf{I}}_{n})$ and are scaled such that $||{\bm{\gamma}}||_{2}=||{\bm{\xi}}||_{2}=1$. Under these circumstances, we find that $\eta=\mathcal{O}(n^{-1/2})$ with high probability for values of $\theta$ that scale as $\theta\asymp n^{\beta}$ , where $\beta\in[0,1/2)$. Therefore, the RFM demonstrates performance equivalence to the noisy linear model with high probability in this context. Figure 1(b) visually represents this relationship, showcasing the training and generalization errors of both the RFM and the noisy linear model. Notably, the results indicate that this performance-wise equivalence is maintained across different activation functions, including ReLU and tanh.

Figure: (a) $\sigma_{*}=\sigma_{ReLU}$
Refer to caption: extracted/5889816/figures/relu_polynomial_equivalence.png

###### Corollary 3 .

### 4.4 Beyond the Noisy Linear Model: High-Order Polynomial Models

In this section, we explore the performance-wise equivalence of the RFM to the noisy polynomial model beyond the linear regime. Specifically, we focus on the model represented by

$$ ${\bm{\omega}}^{T}\hat{\sigma}_{l}({\mathbf{F}}{\mathbf{x}}),$ (20) $$

for $l>2$, where $\hat{\sigma}_{l}$ is the equivalent polynomial function defined in (14).

If the condition in Corollary 3 does not hold, the equivalence between the RFM and the noisy linear model is no longer applicable. In this case, Theorem 2 provides high-order polynomial models (20) equivalent to the RFM. Thus, we delve into the “nonlinear regime” of the RFM under spiked data conditions. While much of the existing theoretical literature has concentrated on the linear regime of the RFM, Theorem 2 suggests that the nonlinear components of random features become increasingly significant when the spike signal ${\bm{\gamma}}$ and label signal ${\bm{\xi}}$ are aligned, indicating a strong input-label correlation. This relationship is illustrated in Figure 2, which demonstrates that the generalization errors of the RFM closely match those of the equivalent polynomial model in the aligned case ($\alpha:={\bm{\gamma}}^{T}{\bm{\xi}}=1$). Notably, although the equivalence between the RFM and the noisy linear model may not hold due to violated conditions for linear equivalence, it is essential to recognize that such equivalence can still be achieved depending on the choice of activation function $\sigma$ and label function $\sigma_{*}$. This nuanced understanding is further characterized in the following remark.

###### Remark 4 .

When we assume that $\eta=\mathcal{O}(n^{-1/4})$, this allows the equivalent activation function to be represented as a third-degree polynomial. Let $\mu_{0},\mu_{1},\mu_{2},\mu_{3}$ and $\tilde{\mu}_{0},\tilde{\mu}_{1},\tilde{\mu}_{2},\tilde{\mu}_{3}$ denote the first four Hermite coefficients (13) of the functions $\sigma$ and $\sigma_{*}$, respectively. In this case, we can express $\mathop{{\mathbb{E}}}[\hat{\sigma}_{l}({\mathbf{f}}_{i}^{T}{\mathbf{x}})y]=
\sum_{j=0}^{3}\frac{1}{j!}\mu_{j}\tilde{\mu}_{j}\eta_{i}^{j}$ due to (113). This formulation reveals that the noisy polynomial model represented in (20) simplifies to a second-degree model if the product $\mu_{3}\tilde{\mu}_{3}=0$. Furthermore, it reduces to the noisy linear model described in (5) under the conditions that both $\mu_{2}\tilde{\mu}_{2}=0$ and $\mu_{3}\tilde{\mu}_{3}=0$. These relationships highlight how specific properties of the Hermite coefficients dictate the behavior of the model, illustrating the nuanced interplay between activation functions and their corresponding polynomial representations.

The implications of Remark 4 are vividly illustrated in Figure 2. In 2(a), we observe that the noisy linear model suffices for the tanh activation function in terms of equivalence, while a high-order polynomial model is necessary for the ReLU activation, given that $\sigma_{*}$ is ReLU. Conversely, 2(b) showcases the same simulation with $\sigma_{*}$ as tanh, demonstrating that a noisy polynomial model is essential to achieve equivalence to the RFM for tanh activation, whereas the noisy linear model remains adequate for ReLU activation. This distinction arises from the specific properties of the Hermite coefficients: for tanh, we have $\mu_{2}=0$ and $\mu_{3}\neq 0$, while for ReLU, $\mu_{2}\neq 0$ and $\mu_{3}=0$. Together, Theorem 2, Remark 4, and Figure 2 underscore the critical influence of the label generation process on the study of equivalent models. This insight represents a novel contribution of our work, emphasizing how different activation functions necessitate distinct modeling approaches to maintain performance equivalence.

Figure: (a) impact of number of samples ($\alpha=1$ and $\theta=n^{1/2}$)
Refer to caption: extracted/5889816/figures/comparison.png

###### Remark 4 .

## 5 Nonlinearity in RFM Benefits High Input-Label Correlation

In this section, we conduct a comparative analysis of the RFM employing various activation functions, focusing specifically on their generalization performance. To establish a benchmark, we first consider the following “optimal” linear activation function for the RFM:

$$ $\displaystyle\sigma_{linear}(x)=a_{0}+a_{1}x,$ (21) $$

where $a_{0},a_{1}\in\mathbb{R}$ are coefficients determined numerically to minimize the generalization error. This “optimal” linear activation serves as a critical reference point for understanding the conditions under which the RFM surpasses traditional linear models in performance.

According to Theorem 2, the RFM with any activation function that meets our assumptions performs equivalently to a noisy polynomial model. Notably, we find that $\eta=\mathcal{O}(n^{-1/4})$ in (15) holds with high probability when $\beta<1/2$ for $\theta\asymp n^{\beta}$ and ${\mathbf{F}}$ is sampled independently of ${\bm{\gamma}}$ and ${\bm{\xi}}$. Therefore, in this context, we can define the following third-degree polynomial as an optimal activation function:

$$ $\displaystyle\sigma_{polynomial}(x)=b_{0}+b_{1}x+\frac{1}{2}b_{2}(x^{2}-1)+ \frac{1}{6}b_{3}(x^{3}-3x)+b_{4}z,$ (22) $$

where $z\sim{\mathcal{N}}(0,1)$ and the coefficients $b_{0},b_{1},b_{2},b_{3},b_{4}\in\mathbb{R}$ are determined numerically to minimize the generalization error, akin to the approach taken for the optimal linear activation in (21).

We proceed to compare the generalization performance of the RFM across various activation functions: linear (21), polynomial (22), ReLU (9), and Softplus (9). Figure 3(a) presents this comparison in relation to the number of training samples. Notably, the generalization error of the RFM aligns closely with that of the equivalent polynomial model, as predicted by Theorem 2. Additionally, we observe the double-descent phenomenon (Belkin et al., ; Geiger et al., ), characterized by an initial decrease in generalization error, followed by an increase, and then a second decrease as the number of samples (or features) increases, particularly for the ReLU and Softplus activation functions. In contrast, this phenomenon is absent for the linear (21) and polynomial (22) activation functions, as their coefficients are numerically optimized. Similar findings have been documented under isotropic data conditions (Wang & Bento, ; Demir & Doğan, ). Interestingly, the RFM with linear activation outperforms both ReLU and Softplus in the mid-range of $k/m$ due to the double-descent behavior exhibited by those two functions. However, throughout the entire range of $k/m$, the polynomial RFM consistently demonstrates lower generalization error than its counterparts, highlighting its superior performance and robustness across different scenarios.

Next, we investigate the impact of alignment, defined as $\alpha:={\bm{\gamma}}^{T}{\bm{\xi}}$ controlling the input-label correlation, on the generalization performance of the RFM across the aforementioned activation functions in Figure 3(b). In the misaligned scenario ($\alpha\leq 0.3$), we observe relatively high generalization errors, with the errors for both linear and polynomial activation functions closely aligned. As alignment increases, generalization errors decrease, reflecting the inherent relationship between inputs and labels. Specifically, for the range $0.3<\alpha\leq 0.6$, we note a gradual separation in the generalization errors of the linear and polynomial functions. In cases of high alignment ($\alpha>0.6$), the polynomial RFM distinctly outperforms the linear RFM, achieving significantly lower generalization errors. This observation corroborates our theoretical predictions that a strong input-label correlation—facilitated by high values of $\alpha$ —enables the RFM to outperform linear models. Consequently, this underscores the advantages of employing nonlinear activation functions in enhancing model performance when the input-label correlation is high.

In Figure 3(c), we examine the influence of the spike magnitude scale, denoted as $\theta$, on the generalization errors of the RFM across the previously discussed activation functions. Notably, we observe that generalization errors for all activation functions decline as the scale of the spike magnitude increases. In the low spike magnitude regime ($\beta<0.4$), both the RFM with linear and polynomial activation functions yield similar generalization errors, which remain lower than those of the RFM utilizing ReLU and Softplus activations. As we transition into the high spike magnitude regime ($\beta\geq 0.4$), the polynomial RFM consistently outperforms the linear RFM, aligning with our expectations. Although our analysis is confined to $\beta<0.5$, we find that the generalization errors for ReLU and Softplus continue to align with those of their equivalent noisy polynomial models, even at higher values of $\beta$. Importantly, our findings indicate that RFM with nonlinear activation functions surpasses the optimal linear model in scenarios characterized by high spike magnitudes. This reinforces our conclusion that nonlinearity is advantageous for the RFM when dealing with anisotropic data, particularly under conditions of strong alignment and elevated spike magnitudes.

## 6 Conclusion

In conclusion, this study explored the random feature model within the context of spiked covariance data, establishing a significant “universality” theorem that underpins our findings. We demonstrated that a strong correlation between inputs and labels is essential for the RFM to surpass linear models in generalization performance. Additionally, we established that the RFM achieves performance equivalent to noisy polynomial models, with the polynomial degree influenced by the magnitude of the input-label correlation. Our numerical simulations corroborated these theoretical findings, highlighting the RFM’s superior effectiveness in environments with strong input-label correlation. This work highlights the critical role of data modeling in understanding and enhancing model performance. Moving forward, we believe that extending the data model presents a promising avenue for future research, allowing for deeper exploration of complex relationships and further optimization of high-dimensional learning strategies.

#### Acknowledgments

We would like to thank Murat A. Erdogdu for beneficial discussions about this work.

This study was supported by the Scientific and Technological Research Council of Turkey (TÜBİTAK) under project number 124E063. We extend our gratitude to TÜBİTAK for its support. S.D. is supported by a KUIS AI Fellowship and a PhD Scholarship (BİDEB 2211) from TÜBİTAK.

## Appendix A Proof of Universality of Random Features for the Spiked Data

Here, we provide a proof of Theorem 4.1 by adapting the proof used by Hu & Lu (). For a given set of training samples $\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{m}$, let ${\mathbf{R}}:=[{\mathbf{r}}_{1},{\mathbf{r}}_{2},\dots,{\mathbf{r}}_{m}]^{T}$ where ${\mathbf{r}}_{i}:=\sigma({\mathbf{F}}{\mathbf{x}}_{i})$ for all $i\in\{1,2,\dots,m\}$. For a given ${\mathbf{R}}$, we then introduce the following perturbed optimization related to the original problem:

$$ $\displaystyle\Phi_{{\mathbf{R}}}(\tau_{1},\tau_{2}):=\inf_{{\bm{\omega}}\in \mathbb{R}^{k}}\left\{\sum_{i=1}^{m}\left({\bm{\omega}}^{T}{\mathbf{r}}_{i}-y_ {i}\right)^{2}+Q({\bm{\omega}})\right\},$ (23) $$

where $Q({\bm{\omega}}):=m\lambda||{\bm{\omega}}||_{2}^{2}+\tau_{1}k{\bm{\omega}}^{T}
\Sigma_{x}{\bm{\omega}}+\tau_{2}k{\bm{\omega}}^{T}\Sigma_{xy}$ with $\Sigma_{x}:=\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}[\sigma({\mathbf{F}}{
\mathbf{x}})\sigma({\mathbf{F}}{\mathbf{x}})^{T}]$ and $\Sigma_{xy}:=\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left[\sigma({\mathbf{F}}
{\mathbf{x}})y\right]$. Here, $Q({\bm{\omega}})$ is a function of $\tau_{1},\tau_{2}$ (in addition to ${\bm{\omega}}$) but we do not include them in the notation for simplicity. Note that $\Phi_{{\mathbf{R}}}(0,0)/m$ is equal to the training error (${\mathcal{T}}_{\sigma}$) of the original problem (3). The additional terms in $Q({\bm{\omega}})$ are related to the generalization error and we discuss them later. Since we are interested in two different functions $\sigma$ and $\hat{\sigma}$, we define ${\mathbf{A}}:=[{\mathbf{a}}_{1},{\mathbf{a}}_{2},\dots,{\mathbf{a}}_{m}]^{T}$ and ${\mathbf{B}}:=[{\mathbf{b}}_{1},{\mathbf{b}}_{2},\dots,{\mathbf{b}}_{m}]^{T}$ where ${\mathbf{a}}_{i}:=\sigma({\mathbf{F}}{\mathbf{x}}_{i})$ and ${\mathbf{b}}_{i}:=\hat{\sigma}({\mathbf{F}}{\mathbf{x}}_{i})$ for all $i\in\{1,2,\dots,m\}$. Then, the training error ${\mathcal{T}}_{\sigma}$ for the RFM with activation $\sigma$, and the training error ${\mathcal{T}}_{\hat{\sigma}}$ for the case of activation $\hat{\sigma}$ converge in probability to the same value $e_{{\mathcal{T}}}\geq 0$ if the following holds

$$ $\displaystyle\frac{\Phi_{{\mathbf{A}}}(0,0)}{m}\stackrel{{\scriptstyle{\mathbb {P}}}}{{\to}}e_{{\mathcal{T}}}\quad\text{and}\quad\frac{\Phi_{{\mathbf{B}}}(0, 0)}{m}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}e_{{\mathcal{T}}}.$ (24) $$

Since $\Phi_{{\mathbf{R}}}(0,0)/m=(k/m)(\Phi_{{\mathbf{R}}}(0,0)/k)$ and $(k/m)\in(0,\infty)$ by assumption (A.4), we focus on $\Phi_{{\mathbf{R}}}(0,0)/k$ for the rest of the proof in order to write the bounds in terms of $k$ only.

Now, we describe how we handle the generalization error. We introduce $\tau_{1}k{\bm{\omega}}^{T}\Sigma_{x}{\bm{\omega}}+\tau_{2}k{\bm{\omega}}^{T}
\Sigma_{xy}$ into $Q({\bm{\omega}})$ to be used in the calculation of the generalization error. Before describing the relationship between the additional terms and the generalization error, we need the following two additional assumptions for the study of generalization error.

##### Additional Assumptions for Generalization Error

- (A.7)
There exists $\tau_{*}=\mathcal{O}(1/\sqrt{k})$ and $s\in(0,\infty)$ such that $Q({\bm{\omega}})/m$ is $s$-strongly convex for all $\tau_{1},\tau_{2}\in(-\tau_{*},\tau_{*})$.
- (A.8)
There exists a limiting function $q^{*}(\tau_{1},\tau_{2})$ such that $\Phi_{{\mathbf{R}}}(\tau_{1},\tau_{2})/k\stackrel{{\scriptstyle{\mathbb{P}}}}{
{\to}}q^{*}(\tau_{1},\tau_{2})$ for all $\tau_{1},\tau_{2}\in(-\tau_{*},\tau_{*})$. Furthermore, the partial derivatives of $q^{*}(\tau_{1},\tau_{2})$ exist at $\tau_{1}=\tau_{2}=0$, denoted as $\frac{\partial}{\partial\tau_{1}}q^{*}(0,0)=\rho^{*}$ and $\frac{\partial}{\partial\tau_{2}}q^{*}(0,0)=\pi^{*}$.

Assumption (A.7) is required to keep the perturbed objective (23) convex and bound the norm of optimal ${\bm{\omega}}$ of (23). Based on the assumption (A.8), the generalization error of the RFM with $\sigma$ can be specified as follows ${\mathcal{G}}_{\sigma}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\rho^{*}-2
\pi^{*}+\mathop{{\mathbb{E}}}_{y}[y^{2}]$. Furthermore, Lemma 10 in Appendix E states that the generalization error ${\mathcal{G}}_{\sigma}$ for the RFM with activation $\sigma$ and the generalization error ${\mathcal{G}}_{\hat{\sigma}}$ for the case of activation $\hat{\sigma}$ converge in probability to the same value $e_{{\mathcal{G}}}\geq 0$ if the following holds

$$ $\displaystyle\frac{\Phi_{{\mathbf{A}}}(\tau_{1},\tau_{2})}{k}\stackrel{{ \scriptstyle{\mathbb{P}}}}{{\to}}\phi\quad\text{and}\quad\frac{\Phi_{{\mathbf{ B}}}(\tau_{1},\tau_{2})}{k}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\phi,$ (25) $$

for any $\phi\in\mathbb{R}$ and any $\tau_{1},\tau_{2}\in(-\tau_{*},\tau_{*})$. The rest of this appendix shows that (25) holds, which then implies the statement of Theorem 4.1. To show the convergence in probability (25), we take advantage of a test function $\varphi(x)$. Specifically, Lemma 11 in Appendix E states that $(\ref{eq:probabilistic_convergence})$ holds if

$$ $\displaystyle\left|\mathop{{\mathbb{E}}}\varphi\left(\frac{1}{k}\Phi_{\mathbf{ A}}\right)-\mathop{{\mathbb{E}}}\varphi\left(\frac{1}{k}\Phi_{\mathbf{B}} \right)\right|\leq\max\left\{||\varphi||_{\infty},||\varphi^{\prime}||_{\infty },||\varphi^{\prime\prime}||_{\infty}\right\}\kappa_{k},$ (26) $$

for some $\kappa_{k}=o(1)$ and every bounded test function $\varphi(x)$ which has bounded first two derivatives. From now on, we use $\Phi_{{\mathbf{R}}}$ to denote $\Phi_{{\mathbf{R}}}(\tau_{1},\tau_{2})$. We can continue as follows

$$ $\displaystyle\left|\mathbb{E}\varphi\left(\frac{1}{k}\Phi_{\bm{A}}\right)- \mathbb{E}\varphi\left(\frac{1}{k}\Phi_{\bm{B}}\right)\right|$ $\displaystyle\leq\mathbb{E}\left|\mathbb{E}_{\backslash\bm{F}}\left[\varphi \left(\frac{1}{k}\Phi_{\bm{A}}\right)\right]-\mathbb{E}_{\backslash\bm{F}} \left[\varphi\left(\frac{1}{k}\Phi_{\bm{B}}\right)\right]\right|,$ $\displaystyle=\mathbb{E}\left[\left|\mathbb{E}_{\backslash\bm{F}}\left[\varphi \left(\frac{1}{k}\Phi_{\bm{A}}\right)\right]-\mathbb{E}_{\backslash\bm{F}} \left[\varphi\left(\frac{1}{k}\Phi_{\bm{B}}\right)\right]\right|\left(\mathds{ 1}_{\mathcal{A}}(\bm{F})+\mathds{1}_{\mathcal{A}^{c}}(\bm{F})\right)\right],$ $\displaystyle\leq\sup_{\bm{F}\in\mathcal{A}}\left|\mathbb{E}_{\backslash\bm{F} }\left[\varphi\left(\frac{1}{k}\Phi_{\bm{A}}\right)\right]-\mathbb{E}_{ \backslash\bm{F}}\left[\varphi\left(\frac{1}{k}\Phi_{\bm{B}}\right)\right] \right|+2\|\varphi\|_{\infty}\mathbb{P}\left(\mathcal{A}^{c}\right),$ (27) $$

where $\mathbb{E}_{\backslash\bm{F}}[\cdot]$ denotes the conditional expectation for a fixed feature matrix $\bm{F}$ and $\mathds{1}_{\mathcal{A}}(\bm{F})$ is a indicator function (1 if ${\mathbf{F}}\in\mathcal{A}$). We refer to $\mathcal{A}$ as the admissible set of feature matrices defined in Appendix B. We also have $\mathbb{P}\left(\mathcal{A}^{c}\right)=o(1)$, proof of which is also given in Appendix B. Now, we need to show

$$ $\displaystyle\begin{aligned} &\sup_{\bm{F}\in\mathcal{A}}\left|\mathbb{E}_{ \backslash\bm{F}}\varphi\left(\frac{1}{k}\Phi_{\bm{A}}\right)-\mathbb{E}_{ \backslash\bm{F}}\varphi\left(\frac{1}{k}\Phi_{\bm{B}}\right)\right|\leq\max \left\{\left\|\varphi^{\prime}\right\|_{\infty},\left\|\varphi^{\prime\prime} \right\|_{\infty}\right\}\kappa_{k},\end{aligned}$ (28) $$

for some $\kappa_{k}=o(1)$. To find a useful bound for (28), we interpolate the path from $\Phi_{\mathbf{A}}$ to $\Phi_{\mathbf{B}}$ following Lindeberg’s method (Lindeberg, ; Korada & Montanari, )

$$ $\displaystyle\Phi_{t}$ $\displaystyle:=\min_{{\bm{\omega}}\in\mathbb{R}^{k}}\left\{\sum_{i=1}^{t}\left ({\bm{\omega}}^{T}{\mathbf{b}}_{i}-y_{i}\right)^{2}+\sum_{i=t+1}^{m}\left({\bm {\omega}}^{T}{\mathbf{a}}_{i}-y_{i}\right)^{2}+Q({\bm{\omega}})\right\},$ (29) $\displaystyle\hat{{\bm{\omega}}}_{t}$ $\displaystyle:=\operatorname*{arg\,min}_{{\bm{\omega}}\in\mathbb{R}^{k}}\left \{\sum_{i=1}^{t}\left({\bm{\omega}}^{T}{\mathbf{b}}_{i}-y_{i}\right)^{2}+\sum_ {i=t+1}^{m}\left({\bm{\omega}}^{T}{\mathbf{a}}_{i}-y_{i}\right)^{2}+Q({\bm{ \omega}})\right\},$ (30) $$

for $0\leq t\leq m$. Then,

$$ $\displaystyle\left|\mathbb{E}_{\backslash\bm{F}}\varphi\left(\frac{1}{k}\Phi_{ \mathbf{A}}\right)-\mathbb{E}_{\backslash\bm{F}}\varphi\left(\frac{1}{k}\Phi_{ \mathbf{B}}\right)\right|\leq\sum_{t=1}^{m}\left|\mathbb{E}_{\backslash\bm{F}} \varphi\left(\frac{1}{k}\Phi_{t}\right)-\mathbb{E}_{\backslash\bm{F}}\varphi \left(\frac{1}{k}\Phi_{t-1}\right)\right|,$ (31) $$

due to triangle inequality. Therefore, we focus on showing the following

$$ $\displaystyle\left|\mathbb{E}_{\backslash\bm{F}}\varphi\left(\frac{1}{k}\Phi_{ t}\right)-\mathbb{E}_{\backslash\bm{F}}\varphi\left(\frac{1}{k}\Phi_{t-1} \right)\right|\leq\max\left\{\left\|\varphi^{\prime}\right\|_{\infty},\left\| \varphi^{\prime\prime}\right\|_{\infty}\right\}\frac{\kappa_{k}}{k},$ (32) $$

since $k/m\in(0,\infty)$. Here, $\Phi_{t}$ and $\Phi_{t-1}$ can be seen as perturbations of a common “leave-one-out” problem defined as

$$ $\displaystyle\Phi_{\backslash t}$ $\displaystyle:=\min_{{\bm{\omega}}\in\mathbb{R}^{k}}\left\{\sum_{i=1}^{t-1} \left({\bm{\omega}}^{T}{\mathbf{b}}_{i}-y_{i}\right)^{2}+\sum_{i=t+1}^{m}\left ({\bm{\omega}}^{T}{\mathbf{a}}_{i}-y_{i}\right)^{2}+Q({\bm{\omega}})\right\},$ (33) $\displaystyle\hat{{\bm{\omega}}}_{\backslash t}$ $\displaystyle:=\operatorname*{arg\,min}_{{\bm{\omega}}\in\mathbb{R}^{k}}\left \{\sum_{i=1}^{t-1}\left({\bm{\omega}}^{T}{\mathbf{b}}_{i}-y_{i}\right)^{2}+ \sum_{i=t+1}^{m}\left({\bm{\omega}}^{T}{\mathbf{a}}_{i}-y_{i}\right)^{2}+Q({ \bm{\omega}})\right\}.$ (34) $$

Note that $t$-th sample is left out in the definition of $\Phi_{\backslash t}$ and ${\bm{\omega}}_{\backslash t}$. Since $\Phi_{t}\approx\Phi_{\backslash t}$, we apply Taylor’s expansion around $\Phi_{\backslash t}$

$$ $\displaystyle\begin{gathered}\varphi\left(\frac{1}{k}\Phi_{t}\right)=\varphi \left(\frac{1}{k}\Phi_{\backslash t}\right)+\frac{1}{k}\varphi^{\prime}\left( \frac{1}{k}\Phi_{\backslash t}\right)\left(\Phi_{t}-\Phi_{\backslash t}\right) +\frac{1}{2k^{2}}\varphi^{\prime\prime}(\zeta)\left(\Phi_{t}-\Phi_{\backslash t }\right)^{2},\end{gathered}$ (36) $$

where $\zeta$ is some value that lies between $\frac{1}{k}\Phi_{t}$ and $\frac{1}{k}\Phi_{\backslash t}$. Similarly, applying Taylor’s expansion to $\varphi\left(\Phi_{t-1}\right)$ around $\Phi_{\backslash t}$, and then subtracting it from the last equation, we get

$$ $\displaystyle\Bigg{|}\mathbb{E}_{\backslash\bm{F}}\left[\varphi\left(\frac{1}{ k}\Phi_{t}\right)\right]-\mathbb{E}_{\backslash\bm{F}}\left[\varphi\left(\frac {1}{k}\Phi_{t-1}\right)\right]\Bigg{|}$ $\displaystyle\leq\frac{\left\|\varphi^{\prime}(x)\right\|_{\infty}}{k}\mathbb{ E}_{\backslash\bm{F}}\left|\mathbb{E}_{t}\left(\Phi_{t}-\Phi_{t-1}\right) \right|+\frac{\left\|\varphi^{\prime\prime}(x)\right\|_{\infty}}{2k}\frac{ \left(\mathbb{E}_{\backslash\bm{F}}\left(\Phi_{t}-\Phi_{\backslash t}\right)^{ 2}+\mathbb{E}_{\backslash\bm{F}}\left(\Phi_{t-1}-\Phi_{\backslash t}\right)^{2 }\right)}{k},$ (37) $$

where $\mathbb{E}_{t}$ denotes the conditional expectation over the random vectors $\left\{\bm{a}_{t},\bm{b}_{t}\right\}$ associated with the $t$-th training sample, while $\left\{\bm{a}_{i},\bm{b}_{i}\right\}_{i\neq t}$ and $\bm{F}$ are fixed.

To bound the terms on the right side of (37), we first define a surrogate optimization problem:

$$ $\displaystyle\Psi_{t}({\mathbf{r}}):=\Phi_{\backslash t}+\min_{{\bm{\omega}} \in\mathbb{R}^{k}}\left\{\frac{1}{2}\left({\bm{\omega}}-\hat{{\bm{\omega}}}_{ \backslash t}\right)^{T}\bm{H}_{\backslash t}\left({\bm{\omega}}-\hat{{\bm{ \omega}}}_{\backslash t}\right)+\left({\bm{\omega}}^{T}{\mathbf{r}}-y_{t} \right)^{2}\right\},$ (38) $$

where

$$ $\displaystyle\bm{H}_{\backslash t}:=$ $\displaystyle 2\sum_{i=1}^{t-1}{\mathbf{b}}_{i}{\mathbf{b}}_{i}^{T}+2\sum_{i=t +1}^{m}{\mathbf{a}}_{i}{\mathbf{a}}_{i}^{T}+\nabla^{2}Q\left(\hat{{\bm{\omega} }}_{\backslash t}\right)$ (39) $$

is the the Hessian matrix of the objective function $\Phi_{\backslash t}$ evaluated at $\hat{{\bm{\omega}}}_{\backslash t}$.
Due to Lemma 12 in Appendix E, we have $\Phi_{t-1}=\Psi_{t}({\mathbf{a}}_{t})$ and $\Phi_{t}=\Psi_{t}({\mathbf{b}}_{t})$, which makes $\Psi_{t}({\mathbf{r}})$ particularly interesting. We continue by simplifying $\Psi_{t}({\mathbf{r}})$ as follows:

$$ $\displaystyle\Psi_{t}({\mathbf{r}}):$ $\displaystyle=\Phi_{\backslash t}+\min_{{\bm{\omega}}\in\mathbb{R}^{k}}\left\{ \frac{1}{2}\left({\bm{\omega}}-\hat{{\bm{\omega}}}_{\backslash t}\right)^{T} \bm{H}_{\backslash t}\left({\bm{\omega}}-\hat{{\bm{\omega}}}_{\backslash t} \right)+\left({\bm{\omega}}^{T}{\mathbf{r}}-y_{t}\right)^{2}\right\},$ (40) $\displaystyle=\Phi_{\backslash t}+\min_{\tau}\min_{{\mathbf{r}}^{T}({\bm{ \omega}}-\hat{{\bm{\omega}}}_{\backslash t})=\tau}\left\{\frac{1}{2}\left({\bm {\omega}}-\hat{{\bm{\omega}}}_{\backslash t}\right)^{T}\bm{H}_{\backslash t} \left({\bm{\omega}}-\hat{{\bm{\omega}}}_{\backslash t}\right)+\left(\hat{{\bm{ \omega}}}_{\backslash t}^{T}{\mathbf{r}}+\tau-y_{t}\right)^{2}\right\},$ (41) $\displaystyle=\Phi_{\backslash t}+\min_{\tau}\left\{\frac{\tau^{2}}{2\nu_{t}({ \mathbf{r}})}+\left(\hat{{\bm{\omega}}}_{\backslash t}^{T}{\mathbf{r}}+\tau-y_ {t}\right)^{2}\right\},\quad\text{ with }\quad\nu_{t}({\mathbf{r}}):={\mathbf{ r}}^{T}\bm{H}_{\backslash t}^{-1}{\mathbf{r}},$ (42) $\displaystyle=\Phi_{\backslash t}+\frac{\left(\hat{{\bm{\omega}}}_{\backslash t }^{T}{\mathbf{r}}-y_{t}\right)^{2}}{2\nu_{t}({\mathbf{r}})+1},$ (43) $\displaystyle\leq\Phi_{\backslash t}+\left(\hat{{\bm{\omega}}}_{\backslash t}^ {T}{\mathbf{r}}-y_{t}\right)^{2},$ (44) $$

where the inequality in the last line is achieved by using $\tau=0$ in (42) while the rest of the steps are trivial. Next, the following is due to Lemma 13 in Appendix E

$$ $\displaystyle\max\left\{\mathbb{E}_{\backslash\bm{F}}\left(\Psi_{t}\left(\bm{b }_{t}\right)-\Phi_{\backslash t}\right)^{2},\mathbb{E}_{\backslash\bm{F}}\left (\Psi_{t}\left(\bm{a}_{t}\right)-\Phi_{\backslash t}\right)^{2}\right\}$ $\displaystyle\leq k^{1-\epsilon}\text{ polylog }k,$ (45) $$

which holds uniformly over ${\mathbf{F}}\in\mathcal{A}$ and $t\in\{1,2,\dots,m\}$ for some $\epsilon>0$ satisfying $\theta^{2}\leq n^{1-\epsilon}$.
Then, we can reach

$$ $\displaystyle\frac{1}{k}\mathbb{E}_{\backslash\bm{F}}\left(\Phi_{t}-\Phi_{ \backslash t}\right)^{2}$ $\displaystyle=\frac{1}{k}\mathbb{E}_{\backslash\bm{F}}\left(\Psi_{t}\left(\bm{ b}_{t}\right)-\Phi_{\backslash t}\right)^{2}=o(1),$ (46) $$

and similarly,

$$ $\displaystyle\frac{1}{k}\mathbb{E}_{\backslash\bm{F}}\left(\Phi_{t-1}-\Phi_{ \backslash t}\right)^{2}$ $\displaystyle=\frac{1}{k}\mathbb{E}_{\backslash\bm{F}}\left(\Psi_{t}\left(\bm{ a}_{t}\right)-\Phi_{\backslash t}\right)^{2}=o(1).$ (47) $$

Next, by Lemma 12 and equation (42), we have

$$ $\displaystyle\mathbb{E}$ ${}_{\backslash\bm{F}}\left|\mathbb{E}_{k}\left(\Phi_{t}-\Phi_{t-1}\right) \right|=\mathbb{E}_{\backslash\bm{F}}\left|\mathbb{E}_{t}\left[\Psi_{t}\left( \bm{b}_{t}\right)-\Psi_{t}\left(\bm{a}_{t}\right)\right]\right|,$ (48) $\displaystyle=\mathbb{E}_{\backslash\bm{F}}\left|\mathbb{E}_{t}\left[\frac{ \left(\hat{{\bm{\omega}}}_{\backslash t}^{T}{\mathbf{b}}_{t}-y_{t}\right)^{2}} {2\nu_{t}({\mathbf{b}}_{t})+1}-\frac{\left(\hat{{\bm{\omega}}}_{\backslash t}^ {T}{\mathbf{a}}_{t}-y_{t}\right)^{2}}{2\nu_{t}({\mathbf{a}}_{t})+1}\right] \right|,$ (49) $\displaystyle\leq\frac{1}{2\bar{\nu}_{t}+1}\mathbb{E}_{\backslash\bm{F}}\left| \mathbb{E}_{t}\left[\left(\hat{{\bm{\omega}}}_{\backslash t}^{T}{\mathbf{b}}_{ t}-y_{t}\right)^{2}-\left(\hat{{\bm{\omega}}}_{\backslash t}^{T}{\mathbf{a}}_{ t}-y_{t}\right)^{2}\right]\right|+\Delta_{t}({\mathbf{a}}_{t})+\Delta_{t}({ \mathbf{b}}_{t}),$ (50) $$

where $\bar{\nu}_{t}:=\mathop{{\mathbb{E}}}_{t}\nu_{t}({\mathbf{a}}_{t})$ and $\Delta_{t}({\mathbf{r}})$ is defined as

$$ $\displaystyle\Delta_{t}({\mathbf{r}}):=\mathbb{E}_{\backslash\bm{F}}\left| \mathbb{E}_{t}\left(\frac{1}{2\nu_{t}({\mathbf{r}})+1}-\frac{1}{2\bar{\nu}_{t} +1}\right)\left(\hat{{\bm{\omega}}}_{\backslash t}^{T}{\mathbf{r}}-y_{t}\right )^{2}\right|,$ (51) $$

for which the following is expected to hold

$$ $\displaystyle\max\left\{\Delta_{t}({\mathbf{a}}_{t}),\Delta_{t}({\mathbf{b}}_{ t})\right\}=o(1),$ (52) $$

due to the concentration of $\nu_{t}({\mathbf{a}}_{t})$ and $\nu_{t}({\mathbf{b}}_{t})$ around $\bar{\nu}_{t}$ (Lemma 17) and the fact that $\frac{1}{2\bar{\nu}_{t}+1}\leq C$ for some $C>0$ by (44). Next, we bound the remaining term in (50) as

$$ $\displaystyle\mathbb{E}_{\backslash\bm{F}}$ $\displaystyle\left|\mathbb{E}_{t}\left[\left(\hat{{\bm{\omega}}}_{\backslash t }^{T}{\mathbf{b}}_{t}-y_{t}\right)^{2}\right]-\mathbb{E}_{t}\left[\left(\hat{{ \bm{\omega}}}_{\backslash t}^{T}{\mathbf{a}}_{t}-y_{t}\right)^{2}\right]\right|,$ (53) $\displaystyle=\mathbb{E}_{\backslash\bm{F}}\left|\hat{{\bm{\omega}}}_{ \backslash t}^{T}\left(\mathop{{\mathbb{E}}}_{t}[{\mathbf{b}}_{t}{\mathbf{b}}_ {t}^{T}]-\mathop{{\mathbb{E}}}_{t}[{\mathbf{a}}_{t}{\mathbf{a}}_{t}^{T}]\right )\hat{{\bm{\omega}}}_{\backslash t}-2\hat{{\bm{\omega}}}_{\backslash t}^{T} \left(\mathop{{\mathbb{E}}}_{t}[{\mathbf{b}}_{t}y_{t}]-\mathop{{\mathbb{E}}}_{ t}[{\mathbf{a}}_{t}y_{t}]\right)\right|$ (54) $\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\mathbb{E}_{\backslash\bm{F}} \left|\hat{{\bm{\omega}}}_{\backslash t}^{T}\left(\mathop{{\mathbb{E}}}_{t}[{ \mathbf{b}}_{t}{\mathbf{b}}_{t}^{T}]-\mathop{{\mathbb{E}}}_{t}[{\mathbf{a}}_{t }{\mathbf{a}}_{t}^{T}]\right)\hat{{\bm{\omega}}}_{\backslash t}\right|+2 \mathbb{E}_{\backslash\bm{F}}\left|\hat{{\bm{\omega}}}_{\backslash t}^{T}\left (\mathop{{\mathbb{E}}}_{t}[{\mathbf{b}}_{t}y_{t}]-\mathop{{\mathbb{E}}}_{t}[{ \mathbf{a}}_{t}y_{t}]\right)\right|$ (55) $\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\mathbb{E}_{\backslash\bm{F}} \left\|\hat{{\bm{\omega}}}_{\backslash t}\right\|^{2}\mathbb{E}_{\backslash\bm {F}}\left\|\mathop{{\mathbb{E}}}_{t}[{\mathbf{b}}_{t}{\mathbf{b}}_{t}^{T}]- \mathop{{\mathbb{E}}}_{t}[{\mathbf{a}}_{t}{\mathbf{a}}_{t}^{T}]\right\|+2 \mathbb{E}_{\backslash\bm{F}}\left\|\hat{{\bm{\omega}}}_{\backslash t}\right\| \mathbb{E}_{\backslash\bm{F}}\left\|\mathop{{\mathbb{E}}}_{t}[{\mathbf{b}}_{t} y_{t}]-\mathop{{\mathbb{E}}}_{t}[{\mathbf{a}}_{t}y_{t}]\right\|$ (56) $\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}o(1),$ (57) $$

where we use triangle inequality to reach (a). Then, we apply Cauchy-Schwarz inequality in order to reach (b). In the last step, we use a bound on $\mathbb{E}_{\backslash\bm{F}}\left\|\hat{{\bm{\omega}}}_{\backslash t}\right\|
^{l}$ (Lemma 18) with the equivalence of covariance matrices (10) and the equivalence of cross-covariance vectors (11) to reach (c). We can use (57) and (52) to bound the right-hand side of (50). Then, (46)-(50) together with (37) imply (32). Finally, (31)-(32) imply (28) and consequently (26), which completes our proof.

## Appendix B Admissible Set of Feature Matrices

A feature matrix ${\mathbf{F}}=[{\mathbf{f}}_{1},{\mathbf{f}}_{2},\dots,{\mathbf{f}}_{k}]^{T}$ is called admissible (${\mathbf{F}}\in\mathcal{A}$) if it satisfies

$$ $\displaystyle\max_{1\leq i,j\leq k}\left|{\mathbf{f}}_{i}^{T}{\mathbf{f}}_{j}+ \theta{\bm{\gamma}}^{T}{\mathbf{f}}_{j}{\mathbf{f}}_{i}^{T}{\bm{\gamma}}- \delta_{ij}\right|$ $\displaystyle\leq\frac{\text{polylog }k}{\sqrt{k}},$ (58) $\displaystyle\|{\mathbf{F}}({\mathbf{I}}_{n}+(\sqrt{1+\theta}-1){\bm{\gamma}}{ \bm{\gamma}}^{T})\|$ $\displaystyle\leq\sqrt{\theta}\text{ polylog }k,$ (59) $$

where $\delta_{ij}$ denote the Kronecker delta. Note that (58) (see Lemma 5) and (59) (see Lemma 6) hold with high probability for ${\mathbf{f}}_{i}\sim{\mathcal{N}}\left(0,\frac{1}{n+\theta}{\mathbf{I}}_{n}\right)$ which is sampled independent of ${\bm{\gamma}}$, ${\bm{\xi}}$ and ${\mathbf{f}}_{j}$ for $i\neq j$.

###### Lemma 5 .

Under assumptions (A.1)-(A.6), the following holds with high probability

$$ $\displaystyle\max_{1\leq i<j\leq k}\left|{\mathbf{f}}_{i}^{T}{\mathbf{f}}_{j}+ \theta{\bm{\gamma}}^{T}{\mathbf{f}}_{j}{\mathbf{f}}_{i}^{T}{\bm{\gamma}}\right|$ $\displaystyle\leq\frac{\text{polylog }k}{\sqrt{k}},$ (60) $\displaystyle\max_{1\leq i\leq k}\left|\|{\mathbf{f}}_{i}\|^{2}+\theta({\bm{ \gamma}}^{T}{\mathbf{f}}_{i})^{2}-1\right|$ $\displaystyle\leq\frac{\text{polylog }k}{\sqrt{k}}.$ (61) $$

###### Proof.

We first apply triangle equality to reach

$$ $\displaystyle\left|{\mathbf{f}}_{i}^{T}{\mathbf{f}}_{j}+\theta{\bm{\gamma}}^{T }{\mathbf{f}}_{j}{\mathbf{f}}_{i}^{T}{\bm{\gamma}}\right|$ $\displaystyle\leq|{\mathbf{f}}_{i}^{T}{\mathbf{f}}_{j}|+\theta|{\bm{\gamma}}^{ T}{\mathbf{f}}_{j}{\mathbf{f}}_{i}^{T}{\bm{\gamma}}|,$ (62) $\displaystyle\left|\|{\mathbf{f}}_{i}\|^{2}+\theta({\bm{\gamma}}^{T}{\mathbf{f }}_{i})^{2}-1\right|$ $\displaystyle\leq\left|\|{\mathbf{f}}_{i}\|^{2}-\frac{n}{n+\theta}\right|+ \theta\left|({\bm{\gamma}}^{T}{\mathbf{f}}_{i})^{2}-\frac{1}{n+\theta}\right|.$ (63) $$

Next, we find probability bounds for each of the terms on the right-hand side of (62) and (63). Let’s start with $|{\mathbf{f}}_{i}^{T}{\mathbf{f}}_{j}|$. Remember that ${\mathbf{f}}_{i}\sim{\mathcal{N}}\left(0,\frac{1}{n+\theta}{\mathbf{I}}_{n}\right)$. Then, for any $l\in\{1,\dots,n\}$, $f_{i,l}$ and $f_{j,l}$ are sub-Gaussian random variables with sub-Gaussian norm bounded by $\frac{C}{\sqrt{n+\theta}}$ for some $C>0$ (Vershynin, , Example 2.5.8). Therefore, $f_{i,l}f_{j,l}$ is a sub-exponential random variable with sub-exponential norm bounded by $\frac{C^{2}}{n+\theta}$ (Vershynin, , Lemma 2.7.7). Since ${\mathbf{f}}_{i}^{T}{\mathbf{f}}_{j}=\sum_{l=1}^{n}f_{i,l}f_{j,l}$, we can apply Bernstein’s inequality (Vershynin, , Theorem 2.8.2)

$$ $\displaystyle{\mathbb{P}}\left(|{\mathbf{f}}_{i}^{T}{\mathbf{f}}_{j}|\geq \epsilon\right)\leq 2e^{-c\min\left(\frac{\epsilon^{2}}{K^{2}n},\frac{\epsilon }{K}\right)},$ (64) $$

where $K=C^{2}/(n+\theta)$ is the sub-exponential norm bound and $c>0$ is a constant. We set $\epsilon=(\log k)^{2}/\sqrt{k}$ and apply union bound so that we get

$$ $\displaystyle{\mathbb{P}}\left(\max_{1\leq i<j\leq k}|{\mathbf{f}}_{i}^{T}{ \mathbf{f}}_{j}|\geq\frac{(logk)^{2}}{\sqrt{k}}\right)\leq c_{1}e^{-(\log k)^{ 2}/c_{2}},$ (65) $$

for some $c_{1},c_{2}>0$. Note that we use $n/k\in(0,\infty)$ and $\theta\leq\sqrt{n}$ to write $n$ and $\theta$ in terms of $k$. Also, we may reach the following inequality using the analogous steps:

$$ $\displaystyle{\mathbb{P}}\left(\max_{1\leq i\leq k}\left|\|{\mathbf{f}}_{i}\|^ {2}-\frac{n}{n+\theta}\right|\geq\frac{(logk)^{2}}{\sqrt{k}}\right)\leq c_{1}e ^{-(\log k)^{2}/c_{2}},$ (66) $$

for some $c_{1},c_{2}>0$. Note that $n/(n+\theta)$ appears in the equation since $\mathop{{\mathbb{E}}}[\|{\mathbf{f}}_{i}\|^{2}]=n/(n+\theta)$.
Now, we focus on $|{\bm{\gamma}}^{T}{\mathbf{f}}_{j}{\mathbf{f}}_{i}^{T}{\bm{\gamma}}|$. First, note that $({\mathbf{f}}_{i}^{T}{\bm{\gamma}})\sim{\mathcal{N}}(0,1/(n+\theta))$ since ${\bm{\gamma}}$ is a deterministic vector with $\|{\bm{\gamma}}\|=1$. Therefore, $({\mathbf{f}}_{i}^{T}{\bm{\gamma}})$ and $({\mathbf{f}}_{j}^{T}{\bm{\gamma}})$ are sub-Gaussian random variables with sub-Gaussian norm bounded by $\frac{C}{\sqrt{n+\theta}}$ for some $C>0$ (Vershynin, , Example 2.5.8) so $({\bm{\gamma}}^{T}{\mathbf{f}}_{j}{\mathbf{f}}_{i}^{T}{\bm{\gamma}})$ is a sub-exponential random variable with sub-exponential norm bounded by $\frac{C^{2}}{n+\theta}$ (Vershynin, , Lemma 2.7.7). Thus, we have sub-exponential tail bound (Vershynin, , Proposition 2.7.1) as follows

$$ $\displaystyle{\mathbb{P}}\left(|{\bm{\gamma}}^{T}{\mathbf{f}}_{j}{\mathbf{f}}_ {i}^{T}{\bm{\gamma}}|\geq\epsilon\right)\leq 2e^{-C\epsilon(n+\theta)},$ (67) $$

for some $C>0$. Setting $\epsilon=(\log k)^{2}/k$, we reach

$$ $\displaystyle{\mathbb{P}}\left(\max_{1\leq i<j\leq k}|{\bm{\gamma}}^{T}{ \mathbf{f}}_{j}{\mathbf{f}}_{i}^{T}{\bm{\gamma}}|\geq\frac{(logk)^{2}}{k} \right)\leq 2e^{-C(\log k)^{2}},$ (68) $$

for some $C>0$. Following the analogous steps, we also get

$$ $\displaystyle{\mathbb{P}}\left(\max_{1\leq i\leq k}\left|({\mathbf{f}}_{i}^{T} {\bm{\gamma}})^{2}-\frac{1}{n+\theta}\right|\geq\frac{(logk)^{2}}{k}\right) \leq 2e^{-C(\log k)^{2}},$ (69) $$

for some $C>0$ since $\mathop{{\mathbb{E}}}[({\mathbf{f}}_{i}^{T}{\bm{\gamma}})^{2}]=1/(n+\theta)$. Using the found probability bounds for each of the terms on the right-hand side of (62) and (63), we complete the proof.
∎

###### Lemma 6 .

Under assumptions (A.1)-(A.6), the following hold with high probability

$$ $\|{\mathbf{F}}({\mathbf{I}}_{n}+\theta{\bm{\gamma}}{\bm{\gamma}}^{T}){\mathbf{ F}}^{T}\|\leq\theta\text{ polylog }k.$ (70) $$

###### Proof.

Using the triangle inequality, we get

$$ $\|{\mathbf{F}}({\mathbf{I}}_{n}+\theta{\bm{\gamma}}{\bm{\gamma}}^{T}){\mathbf{ F}}^{T}\|\leq\|{\mathbf{F}}{\mathbf{F}}^{T}\|+\theta\|{\mathbf{F}}{\bm{\gamma} }{\bm{\gamma}}^{T}{\mathbf{F}}^{T}\|.$ (71) $$

First, we focus on $\|{\mathbf{F}}{\mathbf{F}}^{T}\|=\|{\mathbf{F}}\|^{2}$. Here, we can take advantage of a well-known result on the spectral norm of Gaussian random matrices (Vershynin, , Corollary 7.3.3):

$$ $\displaystyle{\mathbb{P}}\left(\|{\mathbf{F}}\|\geq\sqrt{\frac{n}{n+\theta}}+ \sqrt{\frac{k}{n+\theta}}+\epsilon\right)$ $\displaystyle\leq 2e^{-c(n+\theta)\epsilon^{2}},$ (72) $\displaystyle{\mathbb{P}}\left(\|{\mathbf{F}}\|\geq\sqrt{\frac{n}{n+\theta}}+2 \sqrt{\frac{k}{n+\theta}}\right)$ $\displaystyle\leq 2e^{-ck}$ (73) $$

for any $\epsilon>0$ and some $c>0$. This result indicates $\|{\mathbf{F}}{\mathbf{F}}^{T}\|\leq\text{polylog }k$ with high probability.
Next, we work on $\|{\mathbf{F}}{\bm{\gamma}}{\bm{\gamma}}^{T}{\mathbf{F}}^{T}\|$. Since $({\mathbf{F}}{\bm{\gamma}}{\bm{\gamma}}^{T}{\mathbf{F}}^{T})_{i,j}={\mathbf{f}
}_{i}^{T}{\bm{\gamma}}{\bm{\gamma}}^{T}{\mathbf{f}}_{j}$, we may use the related results (68) - (69) from the proof of Lemma 5. Specifically, we have the following inequality with high probability

$$ $\displaystyle\max_{1\leq i,j\leq k}\left|{\bm{\gamma}}^{T}{\mathbf{f}}_{j}{ \mathbf{f}}_{i}^{T}{\bm{\gamma}}-\frac{\delta_{ij}}{n+\theta}\right|$ $\displaystyle\leq\frac{\text{polylog }k}{k},$ (74) $\displaystyle\max_{1\leq i,j\leq k}\left|{\bm{\gamma}}^{T}{\mathbf{f}}_{j}{ \mathbf{f}}_{i}^{T}{\bm{\gamma}}\right|$ $\displaystyle\leq\frac{\text{polylog }k}{k},$ (75) $$

where $\delta_{ij}$ denote the Kronecker delta function and. we use $1/(n+\theta)<C/k$ for some $C>0$ in reaching the last line. It follows that the following holds with high probability

$$ $\|{\mathbf{F}}{\bm{\gamma}}{\bm{\gamma}}^{T}{\mathbf{F}}^{T}\|\leq\|{\mathbf{F }}{\bm{\gamma}}{\bm{\gamma}}^{T}{\mathbf{F}}^{T}\|_{F}=\sqrt{\sum_{i=1}^{k} \sum_{j=1}^{k}({\mathbf{f}}_{i}^{T}{\bm{\gamma}}{\bm{\gamma}}^{T}{\mathbf{f}}_ {j})^{2}}\leq\text{polylog }k.$ (76) $$

Using (71), we combine the results and reach the statement of the lemma.
∎

###### Lemma 5 .

###### Proof.

###### Lemma 6 .

###### Proof.

## Appendix C Equivalence of Covariance Matrices

Here, we show that (10) holds for any activation function $\sigma$ satisfying assumption (A.6) and a second activation function defined as $\hat{\sigma}(x):=\mu_{1}x+\mu_{*}z$ where $z\sim{\mathcal{N}}(0,1)$. Note that trivially we have $\mathop{{\mathbb{E}}}_{\mathbf{x}}[\hat{\sigma}({\mathbf{F}}{\mathbf{x}})\hat{
\sigma}({\mathbf{F}}{\mathbf{x}})^{T}]=\mu_{1}^{2}{\mathbf{F}}({\mathbf{I}}_{n
}+\theta{\bm{\gamma}}{\bm{\gamma}}^{T}{\mathbf{F}}^{T})+\mu_{*}^{2}{\mathbf{I}
}_{k}$. The following lemma specifies the equivalence of covariance matrices for the activation functions $\sigma$ and $\hat{\sigma}$.

###### Lemma 7 .

Under assumptions (A.1)-(A.6) given in Section 3, the following holds

$$ $\displaystyle\left\|\mathbb{E}_{{\mathbf{x}}}[\sigma({\mathbf{F}}{\mathbf{x}}) \sigma({\mathbf{F}}{\mathbf{x}})^{T}]-(\mu_{1}^{2}\hat{{\mathbf{F}}}\hat{{ \mathbf{F}}}^{T}+\mu_{*}^{2}{\mathbf{I}}_{k})\right\|=o\left(1/\text{polylog } k\right),$ (77) $$

where $\hat{{\mathbf{F}}}:={\mathbf{F}}({\mathbf{I}}_{n}+(\sqrt{1+\theta}-1){\bm{
\gamma}}{\bm{\gamma}}^{T})$.

###### Proof.

(Hu & Lu, , Lemma 5) showed
$\|\mathop{{\mathbb{E}}}_{\mathbf{g}}[\sigma({\mathbf{F}}{\mathbf{g}})\sigma({
\mathbf{F}}{\mathbf{g}})^{T}]-(\mu_{1}^{2}{\mathbf{F}}{\mathbf{F}}^{T}+\mu_{*}
^{2}{\mathbf{I}}_{k})\|=o\left(1/\text{polylog }k\right)$ under ${\mathbf{g}}\sim{\mathcal{N}}(0,{\mathbf{I}}_{n})$ and some other assumptions on ${\mathbf{F}}$ matrix and activation function $\sigma$. Here, we adapt their proof to our case. First, note that we can consider $\sigma(\hat{{\mathbf{F}}}{\mathbf{g}})$ instead of $\sigma({\mathbf{F}}{\mathbf{x}})$ since $\sigma(\hat{{\mathbf{F}}}{\mathbf{g}})$ is equal (in distribution) to $\sigma({\mathbf{F}}{\mathbf{x}})$. Now, our goal is to show $\|\mathop{{\mathbb{E}}}_{\mathbf{g}}[\sigma(\hat{{\mathbf{F}}}{\mathbf{g}})
\sigma(\hat{{\mathbf{F}}}{\mathbf{g}})^{T}]-(\mu_{1}^{2}\hat{{\mathbf{F}}}\hat
{{\mathbf{F}}}^{T}+\mu_{*}^{2}{\mathbf{I}}_{k})\|=o\left(1/\text{polylog }k\right)$. The $(i,j)$-th entry of $\mathop{{\mathbb{E}}}_{\mathbf{g}}[\sigma(\hat{{\mathbf{F}}}{\mathbf{g}})
\sigma(\hat{{\mathbf{F}}}{\mathbf{x}})^{T}]$ is $\mathbb{E}\left[\sigma\left(\hat{{\mathbf{f}}}_{i}^{T}{\mathbf{g}}\right)
\sigma\left(\hat{{\mathbf{f}}}_{j}^{T}{\mathbf{g}}\right)\right]$. Since $(\hat{{\mathbf{f}}}_{i}^{T}{\mathbf{g}},\hat{{\mathbf{f}}}_{j}^{T}{\mathbf{g}})$ are jointly Gaussian, we can rewrite their joint distribution as that of $\left(z_{i},\rho_{ij}z_{i}+\sqrt{1-\rho_{ij}\rho_{ji}}z_{j}\right)$, where $z_{i}\sim\mathcal{N}(0,\|\hat{{\mathbf{f}}}_{i}\|^{2})$, $z_{j}\sim\mathcal{N}(0,\|\hat{{\mathbf{f}}}_{j}\|^{2})$ are two independent Gaussian variables and $\rho_{ij}:=\hat{{\mathbf{f}}}_{i}^{T}\hat{{\mathbf{f}}}_{j}/\|\hat{{\mathbf{f}
}}_{i}\|^{2}$. Then, for $i\neq j$, we have,

$$ $\displaystyle\mathbb{E}$ $\displaystyle\left[\sigma\left(\hat{{\mathbf{f}}}_{i}^{T}{\mathbf{g}}\right) \sigma\left(\hat{{\mathbf{f}}}_{j}^{T}{\mathbf{g}}\right)\right]=\mathbb{E} \left[\sigma\left(z_{i}\right)\sigma\left(\rho_{ij}z_{i}+\sqrt{1-\rho_{ij}\rho _{ji}}z_{j}\right)\right],$ (78) $\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbb{E}[\sigma\left(z_{i} \right)\sigma\left(\sqrt{1-\rho_{ij}\rho_{ji}}z_{j}\right)]+\rho_{ij}\mathbb{E }\left[\sigma\left(z_{i}\right)z_{i}\sigma^{\prime}\left(\sqrt{1-\rho_{ij}\rho _{ji}}z_{j}\right)\right]$ $\displaystyle\qquad+\frac{1}{2}\rho_{ij}^{2}\mathbb{E}\left[\sigma\left(z_{i} \right)z_{i}^{2}\sigma^{\prime\prime}\left(\sqrt{1-\rho_{ij}\rho_{ji}}z_{j} \right)\right]+\frac{1}{6}\rho_{ij}^{3}\mathbb{E}\left[\sigma\left(z_{i}\right )z_{i}^{3}\sigma^{\prime\prime\prime}\left(\zeta_{ij}\right)\right],$ (79) $\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\hat{{\mathbf{f}}}_{i}^{T}\hat{{ \mathbf{f}}}_{j}\mathbb{E}\sigma^{\prime}\left(z_{i}\right)\mathbb{E}\sigma^{ \prime}\left(\sqrt{1-\rho_{ij}\rho_{ji}}z_{j}\right)+\frac{1}{6}\rho_{ij}^{3} \mathbb{E}\left[\sigma\left(z_{i}\right)z_{i}^{3}\sigma^{\prime\prime\prime} \left(\zeta_{ij}\right)\right],$ (80) $\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\hat{{\mathbf{f}}}_{i}^{T}\hat{{ \mathbf{f}}}_{j}\mathbb{E}\sigma^{\prime}\left(z_{i}\right)\mathbb{E}\sigma^{ \prime}\left(z_{j}\right)+R_{ij},$ (81) $$

where $\zeta_{ij}$ is some point between $\sqrt{1-\rho_{ij}\rho_{ji}}z_{j}$ and $\rho_{ij}z_{i}+\sqrt{1-\rho_{ij}\rho_{ji}}z_{j}$. Here, we apply Taylor’s series expansion of $\sigma(\rho_{ij}z_{i}+\sqrt{1-\rho_{ij}\rho_{ji}}z_{j})$ around $\sqrt{1-\rho_{ij}\rho_{ji}}z_{j}$ to reach (a). Then, we use the independence between $z_{i}$ and $z_{j}$, and the following identities: $\mathbb{E}\sigma(z_{i})=\mathbb{E}[\sigma(z_{i})z_{i}^{2}]=0$ (due to $\sigma$ being an odd function) and $\mathbb{E}[\sigma(z_{i})z_{i}]=\|\hat{{\mathbf{f}}}_{i}\|^{2}\mathbb{E}[\sigma
^{\prime}(z_{i})]$ (by Stein’s lemma). To reach (c), we define the remainder term $R_{ij}$ as

$$ $\displaystyle R_{ij}=\hat{{\mathbf{f}}}_{i}^{T}\hat{{\mathbf{f}}}_{j}\mathbb{E }\sigma^{\prime}\left(z_{i}\right)\left(\mathbb{E}\sigma^{\prime}\left(\sqrt{1 -\rho_{ij}\rho_{ji}}z_{j}\right)-\mathbb{E}\sigma^{\prime}\left(z_{j}\right) \right)+\frac{1}{6}\rho_{ij}^{3}\mathbb{E}\left[\sigma\left(z_{i}\right)z_{i}^ {3}\sigma^{\prime\prime\prime}\left(\zeta_{ij}\right)\right].$ (82) $$

For $i=j$, we define $R_{ii}=0$. Then, we can verify the following decomposition using (81):

$$ $\displaystyle\mathop{{\mathbb{E}}}[\sigma(\hat{{\mathbf{F}}}{\mathbf{g}}) \sigma(\hat{{\mathbf{F}}}{\mathbf{g}})^{T}]=\left(\mu_{1}{\mathbf{I}}_{k}+{ \mathbf{D}}_{1}\right)\hat{{\mathbf{F}}}\hat{{\mathbf{F}}}^{T}\left(\mu_{1}{ \mathbf{I}}_{k}+{\mathbf{D}}_{1}\right)+\mu_{*}^{2}{\mathbf{I}}_{k}+{\mathbf{D }}_{2}+{\mathbf{D}}_{3}+\mathbf{R}$ (83) $$

where we define ${\mathbf{D}}_{1},{\mathbf{D}}_{2},{\mathbf{D}}_{3}$ to be diagonal matrices as follows: ${\mathbf{D}}_{1}:=\text{diag }(\mathbb{E}\sigma^{\prime}\left(z_{i}\right)-\mu
_{1})$, ${\mathbf{D}}_{2}=\text{diag }(\mu_{1}^{2}-\|\hat{{\mathbf{f}}}_{i}\|^{2}(
\mathbb{E}\sigma^{\prime}(z_{i}))^{2})$, and ${\mathbf{D}}_{3}=\text{diag }(\mathbb{E}\sigma^{2}(z_{i})-\mu_{1}^{2}-\mu_{*}^
{2})$.

Then, we have

$$ $\displaystyle\left\|\mathop{{\mathbb{E}}}[\sigma({\mathbf{F}}{\mathbf{x}}) \sigma({\mathbf{F}}{\mathbf{x}})^{T}]-(\mu_{1}^{2}\hat{{\mathbf{F}}}\hat{{ \mathbf{F}}}^{T}+\mu_{*}^{2}{\mathbf{I}}_{k})\right\|$ $\displaystyle\leq\left(2\mu_{1}+\right.\left.\left\|\mathbf{D}_{1}\right\| \right)\|\hat{{\mathbf{F}}}\|^{2}\left\|\mathbf{D}_{1}\right\|+\left\|\mathbf{ D}_{2}\right\|+\left\|\mathbf{D}_{3}\right\|+\left\|\mathbf{R}\right\|,$ (84) $\displaystyle\leq\frac{\|\hat{\mathbf{F}}\|^{2}\text{ polylog}k}{\sqrt{k}},$ (85) $\displaystyle=o\left(1/\text{polylog }k\right),$ (86) $$

where Lemma 8 is used for the bounds of $\left\|\mathbf{D}_{1}\right\|,\left\|\mathbf{D}_{2}\right\|,\left\|\mathbf{D}_
{3}\right\|$ and $\|\mathbf{R}\|$. To reach the conclusion, we then take advantage of $\|\hat{{\mathbf{F}}}\|^{2}\leq k^{1/2-\epsilon/2}\text{ polylog }k$ (for some $\epsilon>0$ satisfying $\theta\leq n^{1/2-\epsilon/2}$), which is due to (59).
∎

###### Lemma 8 .

Suppose the setting in Lemma 7 and the definitions in its proof. Then,

$$ $\displaystyle\max\left\{\|{\mathbf{D}}_{1}\|,\|{\mathbf{D}}_{2}\|,\|{\mathbf{D }}_{3}\|,\|\mathbf{R}\|\right\}\leq\frac{\text{polylog }k}{\sqrt{k}}.$ (87) $$

###### Proof.

Before bounding the terms, we define a truncated version of the activation function $\sigma$ as

$$ $\sigma_{\epsilon}(x):=\begin{cases}\sigma(x),&\text{if }|x|<\epsilon,\\ 0,&\text{otherwise,}\end{cases}$ (88) $$

for $\epsilon=2logk$, which is useful due to Lemma 19. Furthermore, $\sigma_{\epsilon}(x)$, its finite powers and its derivatives are bounded due to assumption (A.6).
First, we focus on ${\mathbf{D}}_{1}$

$$ $\displaystyle\|{\mathbf{D}}_{1}\|$ $\displaystyle\leq\max_{i}\mathbb{E}|\sigma^{\prime}\left(z_{i}\right)-\mu_{1}|$ (89) $\displaystyle\leq\max_{i}\mathbb{E}_{z\sim{\mathcal{N}}(0,1)}|\sigma^{\prime}( \|\hat{{\mathbf{f}}}_{i}\|z)-\sigma^{\prime}(z)|$ (90) $\displaystyle\leq\max_{i}\mathbb{E}_{z\sim{\mathcal{N}}(0,1)}|\sigma_{\epsilon }^{\prime}(\|\hat{{\mathbf{f}}}_{i}\|z)-\sigma_{\epsilon}^{\prime}(z)|+\frac{ \text{ polylog }k}{k^{logk}}$ (91) $\displaystyle\leq\left(\|\sigma^{\prime\prime}_{\epsilon}(x)\|_{\infty}\mathop {{\mathbb{E}}}[|z|]\right)\max_{i}\left|\|\hat{{\mathbf{f}}}_{i}\|-1\right|+ \frac{\text{ polylog }k}{k^{logk}}$ (92) $\displaystyle\leq\left(\|\sigma^{\prime\prime}_{\epsilon}(x)\|_{\infty}\mathop {{\mathbb{E}}}[|z|]\right)\max_{i}\left|\|\hat{{\mathbf{f}}}_{i}\|^{2}-1\right |+\frac{\text{ polylog }k}{k^{logk}}$ (93) $\displaystyle\leq\frac{\text{ polylog }k}{\sqrt{k}},$ (94) $$

where we use Lemma 19 to get (91) and while the last line is due to (58) (Lemma 5) and bounded second derivative of $\sigma_{\epsilon}$. Similarly, one can show $\|{\mathbf{D}}_{2}\|\leq\text{ polylog }k/\sqrt{k}$. Next, we study $\|{\mathbf{D}}_{3}\|$ as

$$ $\displaystyle\|{\mathbf{D}}_{3}\|$ $\displaystyle\leq\max_{i}|\mathbb{E}\sigma^{2}\left(z_{i}\right)-\mu_{1}^{2}- \mu_{*}^{2}|$ (95) $\displaystyle=\max_{i}\mathbb{E}_{z\sim{\mathcal{N}}(0,1)}|\sigma^{2}(\|\hat{{ \mathbf{f}}}_{i}\|z)-\sigma^{2}(z)|$ (96) $\displaystyle=\max_{i}\mathbb{E}_{z\sim{\mathcal{N}}(0,1)}|\sigma_{\epsilon}^{ 2}(\|\hat{{\mathbf{f}}}_{i}\|z)-\sigma_{\epsilon}^{2}(z)|+\frac{\text{ polylog }k}{k^{logk}}$ (97) $\displaystyle\leq\text{ polylog }k\max_{i}\left|\|\hat{{\mathbf{f}}}_{i}\|-1 \right|+\frac{\text{ polylog }k}{k^{logk}}$ (98) $\displaystyle\leq\text{ polylog }k\max_{i}\left|\|\hat{{\mathbf{f}}}_{i}\|^{2} -1\right|+\frac{\text{ polylog }k}{k^{logk}}$ (99) $\displaystyle\leq\frac{\text{ polylog }k}{\sqrt{k}},$ (100) $$

where we use the bounded derivative of $\sigma_{\epsilon}^{2}$ to reach (98) while the rest is similar to the bounding of $\|{\mathbf{D}}_{1}\|$. Lastly, to bound $\|\mathbf{R}\|$, one can first easily show that

$$ $\max_{1\leq i,j\leq k}|R_{ij}|\leq\frac{\text{polylog }k}{k\sqrt{k}}.$ (101) $$

using (58) (Lemma 5) and the definition of $R_{ij}$ in (82). Then, we have

$$ $\displaystyle\|\mathbf{R}\|\leq\|\mathbf{R}\|_{F}=\left(\sum_{i=1}^{k}\sum_{j= 1}^{k}R_{ij}^{2}\right)^{1/2}\leq\frac{\text{ polylog }k}{\sqrt{k}},$ (102) $$

which concludes our proof.
∎

###### Lemma 7 .

###### Proof.

###### Lemma 8 .

###### Proof.

## Appendix D Equivalence of Cross-Covariance Vectors

Here, we are interested in the cross-covariance $\mathop{{\mathbb{E}}}[\sigma({\mathbf{F}}{\mathbf{x}})y]$. The following lemma specifies the equivalence of the cross-covariance vectors for activation function $\sigma$ and its finite-degree Hermite expansion $\hat{\sigma}_{l}$.

###### Lemma 9 .

Suppose assumptions (A.1)-(A.6). Let $\eta_{i}:=\mathop{{\mathbb{E}}}\left[\frac{({\mathbf{f}}_{i}{\mathbf{x}})({\bm
{\xi}}^{T}{\mathbf{x}})}{\sqrt{1+\theta\alpha^{2}}}\right]=\frac{{\mathbf{f}}_
{i}^{T}{\bm{\xi}}+\theta\alpha{\mathbf{f}}_{i}^{T}{\bm{\gamma}}}{\sqrt{1+
\theta\alpha^{2}}}$. Assume that

$$ $\eta:=\max_{1\leq i\leq k}\left|\eta_{i}\right|\leq\frac{C}{n^{1/l}},\quad \text{for some $C>0$ and some $l\in\mathbb{Z}^{+}$.}$ (103) $$

Then, the following holds

$$ $\left\|\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left[\sigma({\mathbf{F}}{ \mathbf{x}})y\right]-\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left[\hat{\sigma }_{l}({\mathbf{F}}{\mathbf{x}})y\right]\right\|=o\left(1/\text{polylog }k \right),$ (104) $$

for $\hat{\sigma}_{l}(x)$ defined in (14).

###### Proof.

Let ${\mathbf{x}}\sim{\mathcal{N}}(0,{\mathbf{I}}_{k}+\theta{\bm{\gamma}}{\bm{
\gamma}}^{T})$. By assumption (A.3), we have

$$ $\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left[\sigma({\mathbf{F}}{\mathbf{x}}) y\right]=\mathop{{\mathbb{E}}}_{{\mathbf{x}}}\left[\sigma({\mathbf{F}}{\mathbf {x}})\sigma_{*}\left(\frac{{\bm{\xi}}^{T}{\mathbf{x}}}{\sqrt{1+\theta\alpha^{2 }}}\right)\right]=\mathop{{\mathbb{E}}}_{{\mathbf{g}}}\left[\sigma\left(\hat{{ \mathbf{F}}}{\mathbf{g}}\right)\sigma_{*}\left(\hat{{\bm{\xi}}}^{T}{\mathbf{g} }\right)\right],$ (105) $$

where ${\mathbf{g}}\sim{\mathcal{N}}(0,{\mathbf{I}}_{n})$, $\hat{{\mathbf{F}}}:={\mathbf{F}}({\mathbf{I}}_{n}+(\sqrt{1+\theta}-1){\bm{
\gamma}}{\bm{\gamma}}^{T})$ and $\hat{{\bm{\xi}}}:=\frac{({\mathbf{I}}_{n}+(\sqrt{1+\theta}-1){\bm{\gamma}}{\bm
{\gamma}}^{T}){\bm{\xi}}}{\sqrt{1+\theta\alpha^{2}}}$. Before dealing with the functions $\sigma$ and $\sigma_{*}$, we consider the joint distribution of $(\hat{{\mathbf{f}}}_{i}^{T}{\mathbf{g}},\hat{{\bm{\xi}}}{}^{T}{\mathbf{g}})$, which is jointly Gaussian

$$ $\displaystyle(\hat{{\mathbf{f}}}_{i}^{T}{\mathbf{g}},\hat{{\bm{\xi}}}{}^{T}{ \mathbf{g}})\sim{\mathcal{N}}\left(0,\begin{bmatrix}\|\hat{{\mathbf{f}}}_{i}\| ^{2}&\eta_{i}\\ \eta_{i}&1\end{bmatrix}\right),$ (106) $$

Indeed, we can write an equivalent of $(\hat{{\mathbf{f}}}_{i}^{T}{\mathbf{g}},\hat{{\bm{\xi}}}{}^{T}{\mathbf{g}})$ as follows

$$ $\displaystyle(\eta_{i}z+\sqrt{\|\hat{{\mathbf{f}}}_{i}\|^{2}-\eta_{i}^{2}}z_{i },z),$ (107) $$

where $z_{i},z\stackrel{{\scriptstyle i.i.d}}{{\sim}}{\mathcal{N}}(0,1)$. This decomposition will be useful since we split correlated part $\eta_{i}z$ and uncorrelated part $\sqrt{\|\hat{{\mathbf{f}}}_{i}\|^{2}-\eta_{i}^{2}}z_{i}$. Now, we use Hermite expansions of $\sigma$ and $\sigma_{*}$

$$ $\displaystyle\sigma(x)=\sum_{j=0}^{\infty}\frac{1}{j!}\mu_{j}H_{j}(x),\quad \text{ and }\quad\sigma_{*}(x)=\sum_{j=0}^{\infty}\frac{1}{j!}\tilde{\mu}_{j}H _{j}(x),$ (108) $$

where we define $H_{j}(z)$ to be the probabilist’s $j$-th Hermite polynomial, $\mu_{j}:=\mathop{{\mathbb{E}}}_{z\sim{\mathcal{N}}(0,1)}[H_{j}(z)\sigma(z)]$ and $\tilde{\mu}_{j}:=\mathop{{\mathbb{E}}}_{z\sim{\mathcal{N}}(0,1)}[H_{j}(z)
\sigma_{*}(z)]$.
Then, we have

$$ $\displaystyle\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}$ $\displaystyle\left[\sigma({\mathbf{f}}_{i}^{T}{\mathbf{x}})y\right]=\mathop{{ \mathbb{E}}}_{{\mathbf{g}}}\left[\sigma(\hat{{\mathbf{f}}}_{i}^{T}{\mathbf{g}} )\sigma_{*}(\hat{{\bm{\xi}}}{}^{T}{\mathbf{g}})\right],$ (109) $\displaystyle=\mathop{{\mathbb{E}}}_{z_{i},z}\left[\sigma\left(\eta_{i}z+\sqrt {\|\hat{{\mathbf{f}}}_{i}\|^{2}-\eta_{i}^{2}}z_{i}\right)\sigma_{*}(z)\right]$ (110) $\displaystyle=\mathop{{\mathbb{E}}}_{z_{i},z}\left[\left(\sum_{j=0}^{\infty} \frac{1}{j!}\mu_{j}H_{j}\left(\eta_{i}z+\sqrt{\|\hat{{\mathbf{f}}}_{i}\|^{2}- \eta_{i}^{2}}z_{i}\right)\right)\left(\sum_{j=0}^{\infty}\frac{1}{j!}\tilde{ \mu}_{j}H_{j}(z)\right)\right],$ (111) $\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\sum_{j=0}^{\infty}\frac{1}{(j!) ^{2}}\mu_{j}\tilde{\mu}_{j}\mathop{{\mathbb{E}}}_{z_{i},z}\left[H_{j}\left( \eta_{i}z+\sqrt{\|\hat{{\mathbf{f}}}_{i}\|^{2}-\eta_{i}^{2}}z_{i}\right)H_{j}( z)\right]$ (112) $\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\sum_{j=0}^{\infty}\frac{1}{j!} \mu_{j}\tilde{\mu}_{j}\eta_{i}^{j}$ (113) $\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\left(\sum_{j=0}^{l-1}\frac{1}{j !}\mu_{j}\tilde{\mu}_{j}\eta_{i}^{j}\right)+R_{i},$ (114) $\displaystyle\stackrel{{\scriptstyle(d)}}{{=}}\mathop{{\mathbb{E}}}_{{\mathbf{ x}},y}\left[\hat{\sigma}_{l}({\mathbf{f}}_{i}^{T}{\mathbf{x}})y\right]+R_{i},$ (115) $$

where $R_{i}:=\sum_{j=l}^{\infty}\frac{1}{j!}\mu_{j}\tilde{\mu}_{j}\eta_{i}^{j}$ to reach (c). To reach (a)-(b), we use the orthogonality of Hermite polynomials (Lemma 20). Finally, we get (d) by the definition of $\hat{\sigma}_{l}$. For $|R_{i}|$, we derive an upper bound as follows

$$ $\displaystyle|R_{i}|$ $\displaystyle=\left|\sum_{j=l}^{\infty}\frac{1}{j!}\mu_{j}\tilde{\mu}_{j}\eta_ {i}^{j}\right|,$ (116) $\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\sum_{j=l}^{\infty}\frac{1}{j !}|\mu_{j}\tilde{\mu}_{j}||\eta_{i}^{j}|,$ (117) $\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\frac{C}{k},$ (118) $$

for some $C>0$. The triangle inequality is used to reach (a) while we get (b) using $|\eta_{i}^{l}|=\mathcal{O}(1/k)$ (by the assumption of the lemma). Finally, we have

$$ $\displaystyle\left\|\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left[\sigma({ \mathbf{F}}{\mathbf{x}})y\right]-\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}\left [\hat{\sigma}_{l}({\mathbf{F}}{\mathbf{x}})y\right]\right\|$ $\displaystyle\leq\sqrt{\sum_{i=1}^{k}\left|\mathop{{\mathbb{E}}}[\sigma({ \mathbf{f}}_{i}^{T}{\mathbf{x}})y]-\mathop{{\mathbb{E}}}[\hat{\sigma}_{l}({ \mathbf{f}}_{i}^{T}{\mathbf{x}})y]\right|^{2}}=\sqrt{\sum_{i=1}^{k}R_{i}^{2}} \leq\frac{C}{\sqrt{k}},$ (119) $$

which completes the proof.
∎

###### Lemma 9 .

###### Proof.

## Appendix E Auxiliary Results for the Proofs

In this section, we provide auxiliary results that are used in the proofs given in the previous sections.

###### Lemma 10 .

Suppose assumptions (A.1)-(A.6) in Section 3 and additional assumptions (A.7)-(A.8) in Appendix A. Then,

$$ $\displaystyle{\mathcal{G}}_{\sigma}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to} }e_{{\mathcal{G}}}\quad\textit{ if and only if }\quad{\mathcal{G}}_{\hat{ \sigma}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}e_{{\mathcal{G}}},$ (120) $$

if the following holds,

$$ $\displaystyle\frac{\Phi_{{\mathbf{A}}}(\tau_{1},\tau_{2})}{k}\stackrel{{ \scriptstyle{\mathbb{P}}}}{{\to}}\phi\quad\text{if and only if}\quad\frac{\Phi _{{\mathbf{B}}}(\tau_{1},\tau_{2})}{k}\stackrel{{\scriptstyle{\mathbb{P}}}}{{ \to}}\phi,$ (121) $$

where $\Phi_{{\mathbf{R}}}(\tau_{1},\tau_{2})$ is the perturbed problem defined in (23). Furthermore, ${\mathbf{A}}:=[{\mathbf{a}}_{1},{\mathbf{a}}_{2},\dots,{\mathbf{a}}_{m}]^{T}$ and ${\mathbf{B}}:=[{\mathbf{b}}_{1},{\mathbf{b}}_{2},\dots,{\mathbf{b}}_{m}]^{T}$ where ${\mathbf{a}}_{i}:=\sigma({\mathbf{F}}{\mathbf{x}}_{i})$ and ${\mathbf{b}}_{i}:=\hat{\sigma}({\mathbf{F}}{\mathbf{x}}_{i})$ for all $i\in\{1,2,\dots,m\}$.

###### Proof.

We start by writing ${\mathcal{G}}_{\sigma}$ as follows

$$ $\displaystyle{\mathcal{G}}_{\sigma}:$ $\displaystyle=\mathop{{\mathbb{E}}}_{({\mathbf{x}},y)}(y-\hat{{\bm{\omega}}}_{ \sigma}^{T}\sigma({\mathbf{F}}{\mathbf{x}}))^{2},$ (122) $\displaystyle=\mathop{{\mathbb{E}}}[y^{2}]-2\hat{{\bm{\omega}}}_{\sigma}^{T} \mathop{{\mathbb{E}}}[\sigma({\mathbf{F}}{\mathbf{x}})y]+\hat{{\bm{\omega}}}_{ \sigma}^{T}\mathop{{\mathbb{E}}}[\sigma({\mathbf{F}}{\mathbf{x}})\sigma({ \mathbf{F}}{\mathbf{x}})^{T}]\hat{{\bm{\omega}}}_{\sigma},$ (123) $\displaystyle=\mathop{{\mathbb{E}}}[y^{2}]-2\hat{{\bm{\omega}}}_{\sigma}^{T} \Sigma_{xy}+\hat{{\bm{\omega}}}_{\sigma}^{T}\Sigma_{x}\hat{{\bm{\omega}}}_{ \sigma},$ (124) $\displaystyle=\mathop{{\mathbb{E}}}[y^{2}]-2\pi_{{\mathbf{A}}}+\rho_{{\mathbf{ A}}},$ (125) $$

where we define $\pi_{{\mathbf{A}}}:=\hat{{\bm{\omega}}}_{0}^{T}\Sigma_{xy}$ and $\rho_{{\mathbf{A}}}:=\hat{{\bm{\omega}}}_{0}^{T}\Sigma_{x}\hat{{\bm{\omega}}}_
{0}$ in the last step. Note that $\hat{{\bm{\omega}}}_{\sigma}=\hat{{\bm{\omega}}}_{0}$ and $\hat{{\bm{\omega}}}_{\hat{\sigma}}=\hat{{\bm{\omega}}}_{m}$ by definition (30). Similarly, we can arrive at ${\mathcal{G}}_{\hat{\sigma}}=\mathop{{\mathbb{E}}}[y^{2}]-2\pi_{{\mathbf{B}}}+
\rho_{{\mathbf{B}}}+R$ for activation function $\hat{\sigma}$ where $\pi_{{\mathbf{B}}}:=\hat{{\bm{\omega}}}_{m}^{T}\Sigma_{xy}$ and $\rho_{{\mathbf{B}}}:=\hat{{\bm{\omega}}}_{m}^{T}\Sigma_{x}\hat{{\bm{\omega}}}_
{m}$. Furthermore, we define $R$ as the remainder term as follows

$$ $\displaystyle R$ $\displaystyle:=\hat{{\bm{\omega}}}_{m}^{T}\left(\mathop{{\mathbb{E}}}[\hat{ \sigma}({\mathbf{F}}{\mathbf{x}})\hat{\sigma}({\mathbf{F}}{\mathbf{x}})^{T}]- \Sigma_{x}\right)\hat{{\bm{\omega}}}_{m}-2\hat{{\bm{\omega}}}_{m}^{T}\left( \mathop{{\mathbb{E}}}[\hat{\sigma}({\mathbf{F}}{\mathbf{x}})y]-\Sigma_{xy} \right),$ (126) $$

which can be bounded as

$$ $\displaystyle\mathbb{E}_{\backslash\bm{F}}|R|$ $\displaystyle\leq\mathbb{E}_{\backslash\bm{F}}\left(\left\|\hat{{\bm{\omega}}} _{m}\right\|^{2}\right)\left\|\mathop{{\mathbb{E}}}[\hat{\sigma}({\mathbf{F}}{ \mathbf{x}})\hat{\sigma}({\mathbf{F}}{\mathbf{x}})^{T}]-\Sigma_{x}\right\|+2 \mathbb{E}_{\backslash\bm{F}}\left(\left\|\hat{{\bm{\omega}}}_{m}\right\| \right)\left\|\mathop{{\mathbb{E}}}[\hat{\sigma}({\mathbf{F}}{\mathbf{x}})y]- \Sigma_{xy}\right\|,$ (127) $\displaystyle=o(1),$ (128) $$

where we use a bound on $\mathbb{E}_{\backslash\bm{F}}[\left\|\hat{{\bm{\omega}}}_{m}\right\|^{l}]$ (Lemma 18) with the equivalence of covariance matrices (10) and the equivalence of cross-covariance vectors (11) to reach the last line. Now, to prove (120), we need to show $\rho_{{\mathbf{A}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\rho^{*}=\frac{
\partial}{\partial\tau_{1}}q^{*}(0,0)$ while $\rho_{{\mathbf{B}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\rho^{*}$ and $\pi_{{\mathbf{A}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\pi^{*}=\frac{
\partial}{\partial\tau_{2}}q^{*}(0,0)$ while $\pi_{{\mathbf{B}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\pi^{*}$. We show $\rho_{{\mathbf{A}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\rho^{*}$ in the following while the rest can be shown analogously. First, observe that, for any $\tau_{1},\tau_{2}\in(-\tau_{*},\tau_{*})$,

$$ $\Phi_{{\mathbf{A}}}(\tau_{1},\tau_{2})\leq\Phi_{{\mathbf{A}}}(0,0)+\tau_{1}k \hat{{\bm{\omega}}}_{0}^{T}\Sigma_{x}\hat{{\bm{\omega}}}_{0}+\tau_{2}k\hat{{ \bm{\omega}}}_{0}^{T}\Sigma_{xy},$ (129) $$

which we arrive at by using the definition of the perturbed optimization (23). Then, for any $\tau\in(0,\tau_{*})$,

$$ $\frac{\Phi_{{\mathbf{A}}}(\tau,0)-\Phi_{{\mathbf{A}}}(0,0)}{k\tau}\leq\rho_{{ \mathbf{A}}}\leq\frac{\Phi_{{\mathbf{A}}}(-\tau,0)-\Phi_{{\mathbf{A}}}(0,0)}{- k\tau}.$ (130) $$

For a fixed $\epsilon>0$, there exists $\delta>0$ such that

$$ $\left|\frac{q^{*}(\delta,0)-q^{*}(0,0)}{\delta}-\rho^{*}\right|\leq\epsilon/3,$ (131) $$

due to assumption (A.8) in Appendix A. Focusing on the inequality on the left in (130), we get:

$$ $\displaystyle{\mathbb{P}}(\rho_{{\mathbf{A}}}-\rho^{*}<-\epsilon)$ $\displaystyle\leq{\mathbb{P}}\left(\frac{\Phi_{{\mathbf{A}}}(\delta,0)-\Phi_{{ \mathbf{A}}}(0,0)}{k\delta}-\rho^{*}<-\epsilon\right),$ (132) $\displaystyle\leq{\mathbb{P}}(|\Phi_{{\mathbf{A}}}(\delta,0)/k-q^{*}(\delta,0) |>\delta\epsilon/3)+{\mathbb{P}}(|\Phi_{{\mathbf{A}}}(0,0)/k-q^{*}(0,0)|> \delta\epsilon/3),$ (133) $\displaystyle=0,$ (134) $$

where we use (131) in the middle line and $\Phi_{{\mathbf{A}}}(\tau_{1},\tau_{2})/k\stackrel{{\scriptstyle{\mathbb{P}}}}{
{\to}}q^{*}(\tau_{1},\tau_{2})$ is used in the last line.
Applying the same steps to the inequality on the right in (130), we arrive at ${\mathbb{P}}(\rho_{{\mathbf{A}}}-\rho^{*}>\epsilon)=0$. Therefore, we conclude that $\rho_{{\mathbf{A}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\rho^{*}$. Using the same reasoning that we used to show $\rho_{{\mathbf{A}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\rho^{*}$, we may also show $\rho_{{\mathbf{B}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\rho^{*}$, $\pi_{{\mathbf{A}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\pi^{*}$, and $\pi_{{\mathbf{B}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\pi^{*}$ but we omit them here. Since we have $\rho_{{\mathbf{A}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\rho^{*}$, $\rho_{{\mathbf{B}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\rho^{*}$, $\pi_{{\mathbf{A}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\pi^{*}$, and $\pi_{{\mathbf{B}}}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\to}}\pi^{*}$, we reach (120) by using (125).
∎

###### Lemma 11 .

Suppose the setting in Appendix A. Then, for any $\phi\in\mathbb{R}$,

$$ $\displaystyle\frac{\Phi_{{\mathbf{A}}}(\tau_{1},\tau_{2})}{k}\stackrel{{ \scriptstyle{\mathbb{P}}}}{{\to}}\phi\quad\text{if and only if}\quad\frac{\Phi _{{\mathbf{B}}}(\tau_{1},\tau_{2})}{k}\stackrel{{\scriptstyle{\mathbb{P}}}}{{ \to}}\phi,$ (135) $$

if the following holds:

$$ $\displaystyle\left|\mathop{{\mathbb{E}}}\varphi\left(\frac{1}{k}\Phi_{\mathbf{ A}}\right)-\mathop{{\mathbb{E}}}\varphi\left(\frac{1}{k}\Phi_{\mathbf{B}} \right)\right|\leq\max\left\{||\varphi||_{\infty},||\varphi^{\prime}||_{\infty },||\varphi^{\prime\prime}||_{\infty}\right\}\kappa_{k},$ (136) $$

for some $\kappa_{k}=o(1)$ and every bounded test function $\varphi(x)$ which has bounded first two derivatives.

###### Proof.

The proof is adapted from Section II-D in (Hu & Lu, ). We need to select an appropriate test function to start the proof. For any fixed $\epsilon>0$ and $\phi$, let

$$ $\varphi_{\epsilon}(x):=\mathds{1}_{|x|\geq 3\epsilon/2}*\zeta_{\epsilon/2}(x- \phi),$ (137) $$

where $*$ denotes the convolution operation and $\zeta_{\epsilon}(x):=\epsilon^{-1}\zeta(x/\epsilon)$ is a scaled mollifier. We define a standard mollifier as:

$$ $\zeta(x):=\begin{cases}ce^{-1/(1-x^{2})},&\text{if }|x|<1,\\ 0,&\text{otherwise,}\end{cases}$ (138) $$

for some constant $c$ ensuring $\int_{\mathbb{R}}\zeta(x)dx=1$. We may easily verify that $||\varphi^{\prime}_{\epsilon}(x)||_{\infty}<C/\epsilon$ and $||\varphi^{\prime\prime}_{\epsilon}(x)||_{\infty}<C/\epsilon^{2}$ for some constant $C$. Furthermore,

$$ $\mathds{1}_{|x-\phi|\geq 2\epsilon}\leq\varphi_{\epsilon}(x)\leq\mathds{1}_{|x -\phi|\geq\epsilon}.$ (139) $$

Let $x=\Phi_{\mathbf{A}}/k$ in (139) and take the expectation:

$$ ${\mathbb{P}}(|\Phi_{\mathbf{A}}/k-\phi|\geq 2\epsilon)\leq\mathop{{\mathbb{E}} }[\varphi_{\epsilon}(\Phi_{\mathbf{A}}/k)].$ (140) $$

Similarly, let $x=\Phi_{\mathbf{B}}/k$ in (139) and take the expectation:

$$ $\mathop{{\mathbb{E}}}[\varphi_{\epsilon}(\Phi_{\mathbf{B}}/k)]\leq{\mathbb{P}} (|\Phi_{\mathbf{B}}/k-\phi|\geq\epsilon).$ (141) $$

Using (136) together with the last two equations, we arrive at

$$ ${\mathbb{P}}(|\Phi_{\mathbf{A}}/k-\phi|\geq 2\epsilon)\leq{\mathbb{P}}(|\Phi_{ \mathbf{B}}/k-\phi|\geq\epsilon)+\text{max}\left\{C,C/\epsilon,C/\epsilon^{2} \right\}\kappa_{k},$ (142) $$

which leads to

$$ ${\mathbb{P}}(|\Phi_{\mathbf{A}}/k-\phi|\geq 2\epsilon)\leq{\mathbb{P}}(|\Phi_{ \mathbf{B}}/k-\phi|\geq\epsilon)+\frac{C\kappa_{k}}{\epsilon^{2}},$ (143) $$

for $\epsilon\in(0,1)$ satisfying $\kappa_{k}/\epsilon^{2}=o(1)$. Applying the same steps after switching ${\mathbf{A}}$ with ${\mathbf{B}}$, we get

$$ ${\mathbb{P}}(|\Phi_{\mathbf{B}}/k-\phi|\geq 2\epsilon)\leq{\mathbb{P}}(|\Phi_{ \mathbf{A}}/k-\phi|\geq\epsilon)+\frac{C\kappa_{k}}{\epsilon^{2}}.$ (144) $$

Combining the last two results and letting $\epsilon\to 0^{+}$ with $\kappa_{k}/\epsilon^{2}=o(1)$, we reach (135).
∎

###### Lemma 12 .

Suppose the definitions in Appendix A. Then, the following holds:

$$ $\Phi_{t-1}=\Psi_{t}({\mathbf{a}}_{t}),\qquad\textit{and}\qquad\Phi_{t}=\Psi_{t }({\mathbf{b}}_{t}).$ (145) $$

###### Proof.

We apply Taylor expansion to $\Phi_{t-1}$ around $\Phi_{\backslash t}$:

$$ $\displaystyle\Phi_{t-1}=\Phi_{\backslash t}+\min_{{\bm{\omega}}\in\mathbb{R}^{ k}}\left\{({\bm{\omega}}-\hat{{\bm{\omega}}}_{\backslash t})^{T}\bm{d}_{ \backslash t}+\frac{1}{2}\left({\bm{\omega}}-\hat{{\bm{\omega}}}_{\backslash t }\right)^{T}\bm{H}_{\backslash t}\left({\bm{\omega}}-\hat{{\bm{\omega}}}_{ \backslash t}\right)+\left({\bm{\omega}}^{T}{\mathbf{a}}_{t}-y_{t}\right)^{2} \right\},$ (146) $$

where $\bm{d}_{\backslash t}$ and $\bm{H}_{\backslash t}$ are respectively gradient vector and Hession matrix of $\Phi_{\backslash t}$ evaluated at $\hat{{\bm{\omega}}}_{\backslash t}$. Note that we do not need higher order terms since $\Phi_{\backslash t}$ only involves quadratic terms, which means higher order terms are equal to 0 0. Furthermore, the gradient vector $\bm{d}_{\backslash t}$ is equal to 0 0 due to the first-order optimality condition in (34). This leads us to:

$$ $\displaystyle\Phi_{t-1}=\Phi_{\backslash t}+\min_{{\bm{\omega}}\in\mathbb{R}^{ k}}\left\{\frac{1}{2}\left({\bm{\omega}}-\hat{{\bm{\omega}}}_{\backslash t} \right)^{T}\bm{H}_{\backslash t}\left({\bm{\omega}}-\hat{{\bm{\omega}}}_{ \backslash t}\right)+\left({\bm{\omega}}^{T}{\mathbf{a}}_{t}-y_{t}\right)^{2} \right\}=\Psi_{t}({\mathbf{a}}_{t}).$ (147) $$

Similarly, one can show $\Phi_{t}=\Psi_{t}({\mathbf{b}}_{t})$ using the same reasoning, which is omitted here.
∎

###### Lemma 13 .

Suppose the definitions in Appendix A. Then, the following holds:

$$ $\max\left\{\mathbb{E}_{\backslash\bm{F}}\left(\Psi_{t}\left(\bm{b}_{t}\right)- \Phi_{\backslash t}\right)^{2},\mathbb{E}_{\backslash\bm{F}}\left(\Psi_{t} \left(\bm{a}_{t}\right)-\Phi_{\backslash t}\right)^{2}\right\}\leq\text{ polylog }k.$ (148) $$

###### Proof.

Let ${\mathbf{r}}_{t}={\mathbf{a}}_{t}$ or ${\mathbf{r}}_{t}={\mathbf{b}}_{t}$. Then,

$$ $\displaystyle\mathbb{E}_{\backslash\bm{F}}\left(\Psi_{t}\left({\mathbf{r}}_{t} \right)-\Phi_{\backslash t}\right)^{2}$ $\displaystyle\leq\mathbb{E}_{\backslash\bm{F}}\left(\hat{{\bm{\omega}}}_{ \backslash t}^{T}{\mathbf{r}}_{t}-y_{t}\right)^{4},$ (149) $\displaystyle\leq\sqrt{\mathbb{E}_{\backslash\bm{F}}\left(\left|\hat{{\bm{ \omega}}}_{\backslash t}^{T}{\mathbf{r}}_{t}\right|+1\right)^{8}}\text{ polylog }k,$ (150) $\displaystyle\leq Q\left(\|\hat{{\bm{\omega}}}_{\backslash t}\|\right)\|\hat{{ \mathbf{F}}}\|^{4}\text{ polylog }k,$ (151) $\displaystyle\leq k^{1-\epsilon}\text{ polylog }k$ (152) $$

where $Q(x)$ is some finite-degree polynomial, $\hat{{\mathbf{F}}}:={\mathbf{F}}({\mathbf{I}}_{n}+(\sqrt{1+\theta}-1){\bm{
\gamma}}{\bm{\gamma}}^{T})$ and some $\epsilon>0$ satisfying $\theta^{2}\leq n^{1-\epsilon}$. We use (44) in the first step. Then, we reach the second line using Lemma 14. In the third step, we use Lemma 15. Finally, we use (59) and Lemma 18 to conclude the proof.
∎

###### Lemma 14 .

Suppose the definitions in Appendix A. Then, the following holds:

$$ $\displaystyle\mathbb{E}_{\backslash\bm{F}}\left(\hat{{\bm{\omega}}}_{ \backslash t}^{T}{\mathbf{r}}-y_{t}\right)^{2l}\leq\sqrt{\mathbb{E}_{ \backslash\bm{F}}\left(\left|\hat{{\bm{\omega}}}_{\backslash t}^{T}{\mathbf{r} }\right|+1\right)^{4l}}\text{ polylog }k,$ (153) $$

for $l\in\mathbb{Z}^{+}$.

###### Proof.

Let $y_{t}:=\sigma_{*}(s_{t})$ where $s_{t}:=\frac{{\bm{\xi}}^{T}{\mathbf{x}}_{t}}{\sqrt{1+\theta\alpha^{2}}}$. Then, for any $x\in\mathbb{R}$,

$$ $\displaystyle(x-y_{t})^{2}$ $\displaystyle\leq|x|^{2}+2|x||y_{t}|+|y_{t}|^{2},$ (154) $\displaystyle\leq|x|^{2}+2|x|C\left(|s_{t}|^{K}+1\right)+C^{2}\left(|s_{t}|^{K }+1\right)^{2},$ (155) $\displaystyle\leq\hat{C}(|x|+1)^{2}\left(|s_{t}|^{K}+1\right)^{2},$ (156) $$

for some $C,\hat{C}>0$. To reach the second line, we use $|y_{t}|\leq C\left(|s_{t}|^{K}+1\right)$ due to assumption (A.3). Using this result, we get:

$$ $\displaystyle\mathbb{E}_{\backslash\bm{F}}\left(\hat{{\bm{\omega}}}_{ \backslash t}^{T}{\mathbf{r}}-y_{t}\right)^{2l}$ $\displaystyle\leq\hat{C}^{l}\mathbb{E}_{\backslash\bm{F}}\left[\left(\left| \hat{{\bm{\omega}}}_{\backslash t}^{T}{\mathbf{r}}\right|+1\right)^{2l}\left( \left|s_{t}\right|^{K}+1\right)^{2l}\right],$ (157) $\displaystyle\leq\hat{C}^{l}\sqrt{\mathbb{E}_{\backslash\bm{F}}\left(\left| \hat{{\bm{\omega}}}_{\backslash t}^{T}{\mathbf{r}}\right|+1\right)^{4l}}\sqrt{ \mathbb{E}_{\backslash\bm{F}}\left(\left|s_{t}\right|^{K}+1\right)^{4l}}.$ (158) $$

$s_{t}\sim{\mathcal{N}}(0,1)$ so ${\mathbb{P}}(|s_{t}|>\epsilon)\leq 2e^{-\epsilon^{2}/2}$ (Vershynin, , Eq. (2.10)). Therefore, we can reach (153) by using moment bounds for $|s_{t}|$ from Lemma 16.
∎

###### Lemma 15 .

Suppose assumptions (A.1)-(A.6). Let ${\mathbf{a}}:=\sigma({\mathbf{F}}{\mathbf{x}})$. Furthermore, let $\hat{{\mathbf{F}}}:={\mathbf{F}}({\mathbf{I}}_{n}+(\sqrt{1+\theta}-1){\bm{
\gamma}}{\bm{\gamma}}^{T})$. Then, there exist some $c,C>0$ such that the following holds

$$ $\displaystyle{\mathbb{P}}_{\backslash{\mathbf{F}}}\left(\left|{\bm{\omega}}^{T }{\mathbf{a}}\right|\geq\epsilon\right)$ $\displaystyle\leq 2e^{\frac{-\epsilon^{2}}{c\|{\bm{\omega}}\|^{2}\|\hat{{ \mathbf{F}}}\|^{2}\text{ polylog }k}},$ (159) $\displaystyle\mathbb{E}_{\backslash{\mathbf{F}}}\left(\left|{\bm{\omega}}^{T}{ \mathbf{a}}\right|^{l}\right)$ $\displaystyle\leq(l!)\left(C\|{\bm{\omega}}\|^{2}\|\hat{{\mathbf{F}}}\|^{2} \text{ polylog }k\right)^{l/2},$ (160) $$

for any $l\in\mathbb{Z}^{+}$, any fixed ${\bm{\omega}}\in\mathbb{R}^{k}$ and $\epsilon\geq 0$.

###### Proof.

Let ${\mathbf{g}}\sim{\mathcal{N}}(0,{\mathbf{I}}_{k})$. One can verify that the function $f({\mathbf{g}}):={\bm{\omega}}^{T}{\mathbf{a}}={\bm{\omega}}^{T}\sigma(\hat{{
\mathbf{F}}}{\mathbf{g}})$ is $(\|\sigma^{\prime}\|_{\infty}\|{\bm{\omega}}\|\|\hat{{\mathbf{F}}}\|)$-Lipschitz continuous. Furthermore, $\mathop{{\mathbb{E}}}_{\backslash{\mathbf{F}}}\left({\bm{\omega}}^{T}{\mathbf{
a}}\right)=0$ due to $\sigma(x)$ being an odd function. Furthermore, $\sigma^{\prime}(z)\leq\text{ polylog }k$ with high probability for $z\sim{\mathcal{N}}(0,\hat{C})$ with all $\hat{C}>0$, which allows us to continue with $\|\sigma^{\prime}\|_{\infty}\leq\text{ polylog }k$ by using a truncation argument (similar to Lemma 19). Therefore, we can reach (159)-(160) using Lemma 16.
∎

###### Lemma 16 (Concentration of Lipschitz Functions - (Hu & Lu, 2023 , Theorem 3) ) .

Consider ${\mathbf{g}}\sim{\mathcal{N}}(0,{\mathbf{I}}_{k})$. Then, for any $\kappa$-Lipschitz function $f:\mathbb{R}^{k}\to\mathbb{R}$ and $\epsilon\geq 0$,

$$ ${\mathbb{P}}(|f({\mathbf{g}})-\mathop{{\mathbb{E}}}[f({\mathbf{g}})]|\geq \epsilon)\leq 2e^{-\epsilon^{2}/(4\kappa^{2})}.$ (161) $$

Furthermore, let $x$ be a random variable. If $x$ satisfying ${\mathbb{P}}(|x|>v)\leq ce^{-Cv}$ for some $C,c>0$,

$$ $\mathop{{\mathbb{E}}}\left[|x|^{l}\right]\leq clC^{-l}\int_{0}^{\infty}e^{-v}v ^{l-1}dv=c(l!)C^{-l}\quad\text{for any $l\in\mathbb{Z}^{+}$}.$ (162) $$

Similarly, if $x$ satisfying ${\mathbb{P}}(|x|>v)\leq ce^{-Cv^{2}}$ for some $C,c>0$,

$$ $\mathop{{\mathbb{E}}}\left[|x|^{l}\right]\leq 2c(l!)C^{-l/2}\quad\text{for any $l\in\mathbb{Z}^{+}$}.$ (163) $$

###### Proof.

See (Talagrand, , Theorem 1.3.4).
∎

###### Lemma 17 .

Suppose the definitions in Appendix A. Let $\epsilon>0$ be a constant satisfying $\theta\leq n^{1/2-\epsilon/2}$. For ${\mathbf{r}}_{t}={\mathbf{a}}_{t}$ or ${\mathbf{r}}_{t}={\mathbf{b}}_{t}$, we have

$$ $\displaystyle{\mathbb{P}}(|\nu_{t}({\mathbf{r}}_{t})-\mathop{{\mathbb{E}}}_{t} [\nu_{t}({\mathbf{r}}_{t})]|>\delta)\leq ce^{-Ck^{\epsilon}\delta^{2}/\text{ polylog }k},$ (164) $$

for some $c,C>0$. Furthermore,

$$ $\displaystyle\mathop{{\mathbb{E}}}_{t}[|\nu_{t}({\mathbf{r}}_{t})-\mathop{{ \mathbb{E}}}_{t}[\nu_{t}({\mathbf{r}}_{t})]|^{l}]\leq 2c(l!)(Ck^{\epsilon}/ \text{polylog }k)^{-l/2},$ (165) $$

for some $c,C>0$ and any $l\in\mathbb{Z}^{+}$. Finally,

$$ $\displaystyle\left|\mathop{{\mathbb{E}}}_{t}\nu_{t}({\mathbf{a}}_{t})-\mathop{ {\mathbb{E}}}_{t}[\nu_{t}({\mathbf{b}}_{t})]\right|=o(1/\text{ polylog }k).$ (166) $$

###### Proof.

Consider

$$ $\displaystyle\nu_{t}({\mathbf{r}}_{t})$ $\displaystyle={\mathbf{r}}_{t}^{T}\bm{H}_{\backslash t}^{-1}{\mathbf{r}}_{t}= \frac{{\mathbf{r}}_{t}^{T}\bar{\bm{H}}_{\backslash t}^{-1}{\mathbf{r}}_{t}}{k},$ (167) $$

where $\bar{\bm{H}}_{\backslash t}$ is defined as follows to take the $1/k$ scaling out

$$ $\displaystyle\bar{\bm{H}}_{\backslash t}:=$ $\displaystyle\frac{2}{k}\sum_{i=1}^{t-1}{\mathbf{b}}_{i}{\mathbf{b}}_{i}^{T}+ \frac{2}{k}\sum_{i=t+1}^{m}{\mathbf{a}}_{i}{\mathbf{a}}_{i}^{T}+\frac{1}{k} \nabla^{2}Q\left(\hat{{\bm{\omega}}}_{\backslash t}\right).$ (168) $$

Then, using $s$-strong convexity of $Q({\bm{\omega}})/m$ due to assumption (A.7), we have $\bar{\bm{H}}_{\backslash t}\succeq(m/k)s{\mathbf{I}}_{k}$, which lead to $\|\bar{\bm{H}}_{\backslash t}^{-1}\|\leq(k/m)(1/s)$. Furthermore, we can rewrite ${\mathbf{a}}_{i}=\sigma({\mathbf{F}}{\mathbf{x}}_{i})$ equivalently as $\sigma(\hat{{\mathbf{F}}}{\mathbf{g}})$ for $\hat{{\mathbf{F}}}:={\mathbf{F}}({\mathbf{I}}_{n}+(\sqrt{1+\theta}-1){\bm{
\gamma}}{\bm{\gamma}}^{T})$ and ${\mathbf{g}}\sim{\mathcal{N}}(0,{\mathbf{I}}_{n})$. We have $\|\hat{{\mathbf{F}}}\|\leq k^{1/4-\epsilon/4}\text{ polylog }k$ (for some $\epsilon>0$ satisfying $\theta\leq n^{1/2-\epsilon/2}$), which is due to (59). We also have $\|\sigma^{\prime}\|_{\infty}\leq\text{ polylog }k$ by using a truncation argument (similar to Lemma 19) as mentioned in the proof of Lemma 15. Then, we use (Louart et al., , Lemma 1) with the bounds for $\|\hat{{\mathbf{F}}}\|$ and $\|\bar{\bm{H}}_{\backslash t}^{-1}\|$ to reach (164). Finally, (165) can be obtained using Lemma 16. To show (166), we proceed as

$$ $\displaystyle\left|\mathop{{\mathbb{E}}}_{t}\nu_{t}({\mathbf{a}}_{t})-\mathop{ {\mathbb{E}}}_{t}[\nu_{t}({\mathbf{b}}_{t})]\right|$ $\displaystyle=\frac{1}{k}\left|\mathop{{\mathbb{E}}}_{t}\left[{\mathbf{a}}_{t} ^{T}\bar{\bm{H}}_{\backslash t}^{-1}{\mathbf{a}}_{t}\right]-\mathop{{\mathbb{E }}}_{t}\left[{\mathbf{b}}_{t}^{T}\bar{\bm{H}}_{\backslash t}^{-1}{\mathbf{b}}_ {t}\right]\right|,$ (169) $\displaystyle=\frac{1}{k}\left|\text{Tr}\left(\bar{\bm{H}}_{\backslash t}^{-1} \left(\mathop{{\mathbb{E}}}_{t}\left[{\mathbf{a}}_{t}{\mathbf{a}}_{t}^{T} \right]-\mathop{{\mathbb{E}}}_{t}\left[{\mathbf{b}}_{t}{\mathbf{b}}_{t}^{T} \right]\right)\right)\right|,$ (170) $\displaystyle=o(1/\text{ polylog }k),$ (171) $$

where Tr denotes the trace while we use $\|\bar{\bm{H}}_{\backslash t}^{-1}\|\leq(k/m)(1/s)$ and (10) to reach the last line.
∎

###### Lemma 18 .

Suppose the definitions in Appendix A. Then, the following holds:

$$ $\displaystyle\mathbb{E}_{\backslash{\mathbf{F}}}\|\hat{{\bm{\omega}}}_{t}\|^{l}$ $\displaystyle\leq\text{polylog }k,$ (172) $\displaystyle\mathbb{E}_{\backslash{\mathbf{F}}}\|\hat{{\bm{\omega}}}_{ \backslash t}\|^{l}$ $\displaystyle\leq\text{polylog }k,$ (173) $$

for any $l\in\mathbb{Z}^{+}$ and any $t\in\{0,\dots,m\}$.

###### Proof.

We show (172) here while (173) can be proved similarly. Let ${\mathbf{r}}_{i}:={\mathbf{b}}_{i}$ for $i\in\{1,\dots,k\}$ and ${\mathbf{r}}_{i}:={\mathbf{a}}_{i}$ for $t\in\{k+1,\dots,m\}$. Furthermore, let $Q({\bm{\omega}}):=m\lambda||{\bm{\omega}}||_{2}^{2}+\tau_{1}k{\bm{\omega}}^{T}
\Sigma_{x}{\bm{\omega}}+\tau_{2}k{\bm{\omega}}^{T}\Sigma_{xy}$. Then,

$$ $\displaystyle\frac{1}{m}Q(\hat{{\bm{\omega}}}_{t})$ $\displaystyle\leq\frac{1}{m}\sum_{i=1}^{m}\left(\hat{{\bm{\omega}}}_{t}^{T}{ \mathbf{r}}_{i})-y_{i}\right)^{2}+\frac{1}{m}Q(\hat{{\bm{\omega}}}_{t})\leq \frac{1}{m}\sum_{i=1}^{m}y_{i}^{2}+\frac{1}{m}Q(\bm{0}),$ (174) $$

since $\hat{{\bm{\omega}}}_{t}$ is the optimal solution defined in (30). Furthermore, $Q({\bm{\omega}})/m$ is $s$-strongly convex by assumption (A.7) in Appendix A. Using this fact, we reach

$$ $\displaystyle\frac{1}{m}{\mathbf{Q}}(\hat{{\bm{\omega}}}_{t})$ $\displaystyle\geq\frac{1}{m}Q(\bm{0})+\frac{1}{m}\nabla Q(\bm{0})^{T}\hat{{\bm {\omega}}}_{t}+(s/2)\|\hat{{\bm{\omega}}}_{t}\|^{2},$ (175) $\displaystyle\geq\frac{1}{m}Q(\bm{0})-\frac{1}{m}\|\nabla Q(\bm{0})\|\|\hat{{ \bm{\omega}}}_{t}\|+(s/2)\|\hat{{\bm{\omega}}}_{t}\|^{2}.$ (176) $$

We may then combine the found lower and upper bounds for $Q(\hat{{\bm{\omega}}}_{t})$

$$ $(s/2)\|\hat{{\bm{\omega}}}_{t}\|^{2}-\frac{1}{m}\|\nabla Q(\bm{0})\|\|\hat{{ \bm{\omega}}}_{t}\|\leq\frac{1}{m}\sum_{i=1}^{m}y_{i}^{2},$ (177) $$

which leads to

$$ $\displaystyle\|\hat{{\bm{\omega}}}_{t}\|$ $\displaystyle\leq\frac{\|\nabla Q(\bm{0})\|/m+\sqrt{\|\nabla Q(\bm{0})\|^{2}/m ^{2}+2(s/m)\sum_{i=1}^{m}y_{i}^{2}}}{s},$ (178) $\displaystyle\leq\frac{2\|\nabla Q(\bm{0})\|/m+\sqrt{2(s/m)\sum_{i=1}^{m}y_{i} ^{2}}}{s},$ (179) $\displaystyle\leq\left(\frac{1}{m}\sum_{i=1}^{m}y_{i}^{2}\right)^{1/2}\text{ polylog }k,$ (180) $$

where we use $\|\nabla Q(\bm{0})\|/m=|\tau_{2}|(k/m)\|\Sigma_{xy}\|\leq\text{polylog }k$ since we have $(k/m)\in(0,\infty)$ by assumption (A.4) in Section 3, $\tau_{2}=\mathcal{O}(1/\sqrt{k})$ by assumption (A.7) in Appendix A and one can easily show $\|\Sigma_{xy}\|\leq\sqrt{k}\text{ polylog }k$ using (113). Then, we consider $l$-th power of $\|\hat{{\bm{\omega}}}_{t}\|$ as

$$ $\|\hat{{\bm{\omega}}}_{t}\|^{l}\leq\left(\frac{1}{m}\sum_{i=1}^{m}y_{i}^{2} \right)^{l/2}\text{ polylog }k.$ (181) $$

Let $y_{i}:=\sigma_{*}(s_{i})$ where $s_{i}:=\frac{{\bm{\xi}}^{T}{\mathbf{x}}_{i}}{\sqrt{1+\theta\alpha^{2}}}$. Using $|y_{i}|\leq C\left(|s_{i}|^{K}+1\right)$ due to assumption (A.3), we get

$$ $\displaystyle\mathbb{E}_{\backslash{\mathbf{F}}}\|\hat{{\bm{\omega}}}_{t}\|^{l}$ $\displaystyle\leq\mathbb{E}_{\backslash{\mathbf{F}}}\left(\frac{1}{m}\sum_{i=1 }^{m}\left(|s_{i}|^{K}+1\right)^{2}\right)^{l/2}\text{ polylog }k.$ (182) $$

Finally, $s_{i}\sim{\mathcal{N}}(0,1)$ so we reach (172) by using moment bounds for $|s_{i}|$ from Lemma 16. Similarly, (173) can shown using analogous steps.
∎

###### Lemma 19 .

Let $f:\mathbb{R}\to\mathbb{R}$ be a function satisfying $|f(x)|\leq C(1+|x|^{K})$ for all $x\in\mathbb{R}$ and some constants $C>0,K\in\mathbb{Z}^{+}$. Furthermore, define a truncated version of $f$ as

$$ $\displaystyle f_{\epsilon}(x):=\begin{cases}f(x),&\text{if }|x|<\epsilon,\\ 0,&\text{otherwise,}\end{cases}$ (183) $$

for $\epsilon>0$. If $\epsilon=2logk$, the following holds for any constant $a>0$ and $z\sim{\mathcal{N}}(0,1)$,

$$ $\mathbb{E}_{z}|f(az)-f_{\epsilon}(az)|=\frac{\text{ polylog }k}{k^{logk}}.$ (184) $$

###### Proof.

We start by using $1=\mathds{1}_{|z|<\epsilon}+\mathds{1}_{|z|\geq\epsilon}$ as follows

$$ $\displaystyle\mathbb{E}_{z}|f(az)-f_{\epsilon}(az)|$ $\displaystyle\leq\mathbb{E}_{z}[|f(az)-f_{\epsilon}(az)|(\mathds{1}_{|z|< \epsilon}+\mathds{1}_{|z|\geq\epsilon})],$ (185) $\displaystyle=\mathbb{E}_{z}[|f(az)-f_{\epsilon}(az)|\mathds{1}_{|z|\geq \epsilon}],$ (186) $\displaystyle=\mathbb{E}_{z}[|f(az)|\mathds{1}_{|z|\geq\epsilon}]$ (187) $\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\sqrt{\mathbb{E}_{z}[f(az)^{2 }]}\sqrt{\mathbb{E}_{z}[(\mathds{1}_{|z|\geq\epsilon})^{2}]},$ (188) $\displaystyle=\sqrt{\mathbb{E}_{z}[f(az)^{2}]}\sqrt{{\mathbb{P}}(|z|\geq \epsilon)},$ (189) $\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}C\sqrt{\mathbb{E}_{z}[(a^{2K} z^{2K}+2a^{K}|z|^{K}+1)]}\sqrt{2e^{-\epsilon^{2}/2}},$ (190) $\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}\frac{\text{ polylog }k}{e^{ log^{2}(k)}}=\frac{\text{ polylog }k}{k^{logk}},$ (191) $$

where we use Cauchy–Schwarz inequality to reach (a) while (b) is due to $|f(x)|\leq C(1+|x|^{K})$ and ${\mathbb{P}}(|z|>\epsilon)\leq 2e^{-\epsilon^{2}/2}$ (Vershynin, , Eq. (2.10)). Finally, we reach (c) by using moment bounds for $|z|$ from Lemma 16.
∎

###### Lemma 20 .

Let $H_{i}(x)$ be the probabilist’s $i$-th Hermite polynomial and $z_{1},z_{2}\stackrel{{\scriptstyle i.i.d}}{{\sim}}{\mathcal{N}}(0,1)$. Then, for any $\rho\in[0,1]$ and any $c\geq|\rho|$,

$$ $\displaystyle\mathbb{E}_{z_{1},z_{2}}[H_{i}(\rho z_{1}+\sqrt{c^{2}-\rho^{2}}z_ {2})H_{j}(z_{1})]=(i!)\rho^{i}\delta_{i,j},$ (192) $$

where $\delta_{ij}$ is Kronecker delta function.

###### Proof.

Let $a=(\rho/c)z_{1}+\sqrt{1-\rho^{2}/c^{2}}z_{2}$ and $b=z_{1}$. Observe that $(a,b)$ is jointly Gaussian with $\mathop{{\mathbb{E}}}[a^{2}]=\mathop{{\mathbb{E}}}[b^{2}]=1$ and $\mathop{{\mathbb{E}}}[ab]=(\rho/c)$. Furthermore, let $p(a,b)$ be the joint probability density function of $(a,b)$ while $p(a)$ and $p(b)$ denoting the individual probability density functions. Using Mehler’s formula (Kibble, ; Mehler, ), we get

$$ $\displaystyle p(a,b)=p(a)p(b)\sum_{l=0}^{\infty}\frac{(\rho/c)^{l}}{l!}H_{l}(a )H_{l}(b).$ (193) $$

Then, it follows

$$ $\displaystyle\mathbb{E}_{z_{1},z_{2}}$ $\displaystyle[H_{i}(\rho z_{1}+\sqrt{c^{2}-\rho^{2}}z_{2})H_{j}(z_{1})]$ (194) $\displaystyle=\mathbb{E}_{a,b}[H_{i}(ca)H_{j}(b)],$ (195) $\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}H_{i}(ca)H_{j}(b)p (a,b)dadb,$ (196) $\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}H_{i}(ca)H_{j}(b)p (a)p(b)\sum_{l=0}^{\infty}\frac{(\rho/c)^{l}}{l!}H_{l}(a)H_{l}(b)dadb,$ (197) $\displaystyle=\sum_{l=0}^{\infty}\frac{(\rho/c)^{l}}{l!}\int_{-\infty}^{\infty }H_{i}(ca)H_{l}(a)p(a)da\int_{-\infty}^{\infty}H_{j}(b)H_{l}(b)p(b)db,$ (198) $\displaystyle=\sum_{l=0}^{\infty}\frac{(\rho/c)^{l}}{l!}\mathop{{\mathbb{E}}}_ {a}[H_{i}(ca)H_{l}(a)]\mathop{{\mathbb{E}}}_{b}[H_{j}(b)H_{l}(b)],$ (199) $\displaystyle=\sum_{l=0}^{\infty}\frac{(\rho/c)^{l}}{l!}(i!)(j!)(c^{i})\delta_ {i,l}\delta_{j,l},$ (200) $\displaystyle=(i!)\rho^{i}\delta_{ij},$ (201) $$

where we use Lemma 21 to reach the line before the last one.
∎

###### Lemma 21 .

Let $H_{i}(x)$ be the probabilist’s $i$-th Hermite polynomial. Also, let $z\sim{\mathcal{N}}(0,1)$. Then, for any $c\in\mathbb{R}$,

$$ $\mathop{{\mathbb{E}}}[H_{i}(cz)H_{j}(z)]=(i!)(c^{i})\delta_{ij}.$ (202) $$

###### Proof.

It is known that $\mathop{{\mathbb{E}}}[H_{i}(z)H_{j}(z)]=(i!)\delta_{ij}$ (O’Donnell, , Chapter 11.2). Here, we extend it to our case. For $i\neq j$, $\mathop{{\mathbb{E}}}[H_{i}(cz)H_{j}(z)]=0$ by orthogonality. Let $a_{l}$ be the coefficient of $x^{l}$ in $H_{i}(x)$. Then, we have

$$ $\displaystyle\mathop{{\mathbb{E}}}[H_{i}(cz)H_{i}(z)]$ $\displaystyle=\sum_{l=0}^{i}c^{l}\mathop{{\mathbb{E}}}[a_{l}z^{l}H_{i}(z)],$ (203) $\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}c^{i}\mathop{{\mathbb{E}}}[a_{l} z^{i}H_{i}(z)],$ (204) $\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}c^{i}\mathop{{\mathbb{E}}}[H_{i} (z)H_{i}(z)],$ (205) $\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}(i!)(c^{i}),$ (206) $$

where we use $\mathop{{\mathbb{E}}}[z^{l}H_{i}(z)]=\delta_{lj}$ due to orthogonality to reach (a) and (b) while (c) is because of $\mathop{{\mathbb{E}}}[H_{i}(z)H_{i}(z)]=(i!)$. Combining $i=j$ and $i\neq j$ cases, we get (202).
∎

###### Lemma 10 .

###### Proof.

###### Lemma 11 .

###### Proof.

###### Lemma 12 .

###### Proof.

###### Lemma 13 .

###### Proof.

###### Lemma 14 .

###### Proof.

###### Lemma 15 .

###### Proof.

###### Lemma 16 (Concentration of Lipschitz Functions - (Hu & Lu, 2023 , Theorem 3) ) .

###### Proof.

###### Lemma 17 .

###### Proof.

###### Lemma 18 .

###### Proof.

###### Lemma 19 .

###### Proof.

###### Lemma 20 .

###### Proof.

###### Lemma 21 .

###### Proof.

## Appendix F Plots for Training Errors

For the sake of completeness, here, we provide the supplementary training error plots corresponding to the setting of Figure 1(b) and Figure 2 in Section 4.

Figure: Figure 4: Misaligned Case (the setting in Figure 1(b))
Refer to caption: extracted/5889816/figures/linear_equivalence_training.png

Figure: (a) $\sigma_{*}=\sigma_{ReLU}$
Refer to caption: extracted/5889816/figures/relu_polynomial_equivalence_training.png

## Peer Reviews

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: My main concern for the present submission is the lack of clear elements of novelty that do not meet the high ICLR standards. 

As the auhtors correctly report, (Moniri et al. 2024) have already provided the rigorous Random Matrix Theory characterization when the spike appears in the weights of the Random Feature map. These results have been extended up to the maximal learning rate scaling regime by (Cui et al. 2024) which describes the emergence of a fully non-polynomial equivalent feature map in this regime. 

I fail to see notable differences between the setting of the present submission and the one in the above-mentioned works. 

In different parts of the manuscript, the authors significantly overclaim their contribution. For example, on page 5 "new universality theorem", "this result is notably more general than previous findings". Could the authors clarify what are the novel aspect in their contribution and distinguish them clearly from the results in (Moniri et al. 2024)?
soundness: 1
strengths: The paper is nicely written and the authors provide a nice introduction to the related literature.
summary: This work examines Random Feature Models (RFMs) under the assumption of spiked covariance for the input data. The authors establish a Gaussian Equivalence (GE) principle for the errors achieved by RFMs in this setting by revealing an equivalence with noisy polynomial equivalent models. When the alignment between the spike and the target is small, the classical GE principle holds, but for larger alignments, an extended version is required. The analysis is supported by numerical simulations that validate the theoretical findings.
weaknesses: The main weakness for this work is the strong relationship with already publsihed works, namely [Moniri et al. 2024, Cui et al. 2024]. The authors overclaim the depth of their contribution in different parts of the manuscript.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 4
questions: Please see the weakness section.
soundness: 4
strengths: - The paper is well written and well organized.

- Theorem 1 is something new and interesting. Although the proof technique is not novel, and is similar to the approach of Hu and Lu, the such universality across models is not something explicitly studied before (to the best of my knowledge). It can be of independent interest.
summary: The paper studies a random features (RF) model trained using ridge regression on a data with a spiked covariance matrix where the covariance is aligned to the target function. Conditions under which such alignment break the gaussian university  (and hence, RF outperforms linear models) is studied.
weaknesses: 1- Although the analysis of RF models with a spiked covariance assumption for the covariates is new, unlike what is stated in line 143-144, it can still be seen as "feature-learning": One can think of spiked covariance as a model for feature learning in the first layers of a deep neural network where the second-to-last layer is random (i.e., the matrix F) and the last layer is trained with ridge regression (i.e., the vector \omega). In particular, the connections to the following papers should be discussed. These papers study layer-wise updates for a three-layer neural network.

- [R1] Eshaan Nichani, Alex Damian, Jason D. Lee, Provable Guarantees for Nonlinear Feature Learning in Three-Layer Neural Networks.
 
- [R2] Zihao Wang, Eshaan Nichani, Jason D Lee, Learning hierarchical polynomials with three-layer neural networks

2- Ba et al, (2023); Mousavi-Hosseini et al. (2023) study a problem similar to the one studied here, but for kernel ridge regression instead of random features regression. There needs to be a detailed comparison of the results presented here to the results of these two papers

3- What happens when we set \beta = 1/2 in Assumption A.2? Also, as the paper studies squared losses, there is probably no need to make the odd activation function assumption (A.6); there can be a easier direct proof of universality using Lindeberg exchange without the need to use the results of Hu & Lu.

4- A precise characterization of the training/test errors will improve the quality of this work significantly.

---

confidence: 3
contribution: 4
initial_rating: 5
presentation: 3
questions: I am a bit confused about why the RFM outperforms the linear model in mathematical expression. Could the authors provide a clear comparison of their learning rates in the presence of strong input-label correlation?
soundness: 2
strengths: This paper solidly verified an interesting phenomenon within the framework of the random feature model that the performance of a learning algorithm has essential dependence on the input or the input-label correlation.  I find it interesting that in the case of a strong correlation between inputs and labels, the RFM outperforms the linear model.
summary: This paper studied whether, and if so under what conditions, the performance of the random feature method (RFM) is better than that of a 
 simpler linear model. The authors considered the spiked covariance data model, where the anisotropic characteristics of $x$ relax the previous work based on a more restrictive isotropic data assumption, thereby broadening the applicability of their established results. One important discovery in this paper is that in the case of a strong correlation between inputs and labels, the RFM outperforms the linear model in the sense that the latter has worse generalization performance.
weaknesses: **Writing:**

 I strongly suggest that the authors include some definitions, symbols, and notations (e.g. the training and generalization errors) in the main text so that the paper can be easily followed. I think that there is significant room for improvement in the paper's arrangement.

**Typographical remarks:**

The first expectation in Eq.(10) should not be taken over $(x,y)$ since the quantity involved does not contain $y$.  Some similar concerns also appear elsewhere.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - How realistic is the spike covariance model? Are there any datasets that (approximately) follow such distribution?
- For the assumptions (equation 10) and (equation 11), do existing datasets satisfy them? If so, is there any references or empirical plots? Also, how does these assumptions compared to those of Hu and Lu (2023)'s?
- In (10), is the equation independent of y?
- It is stated in the literature review that RFMs are used in explaining the double descent phenomenon. Can this study on generalization performance provide information about the double descent phenomenon for data drawn from the spike covariance model?
soundness: 3
strengths: - The paper is clearly written and motivated
- The topic of study is important, since RFMs are closely related to other machine learning models like neural networks. 
- The work extents existing studies by performing analysis on a more general data distribution, which is a step toward more general analysis.
- The theoretical analysis is derived in detail.
- The theoretical analysis is verified with numerical experiments
summary: This paper targets to address the question of "When and how does the RFM outperform linear models?". To this end, the paper considers the setting where the data are distributed according to a spike covariance model, which is more general than the isotropic setting existing analyses based on. The author(s) show that RFM can outperform linear models when there is a high correlation between inputs and labels. They also show that RFM is equivalent to noisy polynomial models, where the polynomial degree depends on the input-output correlations.
weaknesses: - Is it unclear how realistic the assumptions in p.4 are. It would be helpful to verify them using real data and compare the assumptions with the ones used in related analyses
- It seems that every single equation in the paper is numbered. It would be more readable to remove the numbers of the unreferenced equations.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: 1. While the authors say *"it is noteworthy that while ReLU (9) does not conform to the odd function assumptions stipulated in(A.6)"* in line 208, the simulation uses ReLU as an activation function and gets good results, can you explain it more intuitively?
2. The alignment parameter $\alpha$ seems to be a simple multiplication of two parameters in the structure of x and y, what is the real meaning of the input-label correlation in practice?
3. The main theorem in this paper is only an asymptotic results, i.e., equation (56), can you provide some non-asymptotic results such as the converge speed between the corresponding generalization errors $G_{\sigma}$ and $G_{\hat{\sigma}}$? 
4. While this paper is the extension to [1], and [1] considers a general loss function, can the framework from square loss extend to a general loss?

[1]  *Universality Laws for High-Dimensional Learning with Random Features*. (2023), Hu and Lu, TIT.
soundness: 2
strengths: 1. This paper looks good and well-written. Random features model and its connection to neural networks is of interest in machine learning theory.
2. The main contribution is its theoretical aspect, which extends the previous work in [1] to the anisotropic data structure. The universality theorem extends the understanding of RFMs, demonstrating their performance across different activation functions.
3.  Authors suggest that exploiting the non-linearities of RFMs could lead to significant performance gains in high correlation scenarios, which is interesting and confirms practitioners’ intuition.
summary: This paper explores the conditions when and how random feature models (RFM) outperform linear models, focusing on the anisotropic conditions with spiked covariance data with strong input-label correlations. Specifically, they prove that RFM performs equivalent to noisy polynomial models with polynomial degrees influenced by input-label correlation. Numerical simulations are conducted to validate the theoretical results.
weaknesses: 1. **Assumptions**: The author only discusses these assumptions from a technical or the math proof aspect, the practical aspects and their limitations should be clarified. And it seems hard to verify these assumptions in practice.
2. **Comparisons with previous work**:  The comparison with previous work is not enough, the proof technique in this paper is based on [1], and more comparisons and detailed differences should be remarked in the paper.
3.  **Limited experiments**:  While the focus of this paper is theory, the experiments in this paper are not adequate, only evaluated in some simulated data, while the real data application with unknown data structure is missing. The detailed setup of the experiment is not clearly presented in this paper, which should be put in a separate section.
3.  **Scope of Application**: While detailed, the focus is specifically on spiked covariance data, which may not generalize to all data structures and might not hold in all real-world scenarios.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. To what extent would these findings generalise to other anisotropic data settings, and to what extent are they expected to be particular to the spiked covariance model?
2. Are there any practical implications from this work for model selection, given access to a labelled dataset? 
3. Given the odd activation function assumption (A.6), do the authors see any empirical difference in results depending on the parity of $\sigma$? As mentioned above, including $\tanh$ in Fig. 3 might help answer this question.

I sincerely thank the authors for their time and efforts. It’s great to see such detailed, careful mathematical work, which is sometimes lacking in the field! If they can convince me of the broader applicability of their findings (especially, relating to real-world datasets) and clarify assumption A.7 I will be happy to raise my score.
soundness: 3
strengths: The paper makes solid technical contributions, making a careful detailed case that, for spiked covariance data (Eq. 6), RFMs will perform better except when the spike magnitude, alignment parameter, or cosine similarity between the rows of $\mathbf{F}$ and $\gamma$ are sufficiently small (how should I interpret the last condition?). It provides a nice extension to the work of Hu (2023) beyond the linear regime, and it seems possible that this type of data structure/anisotropy might be responsible for the superior performance of RFMs in the wild. The experiments (especially, measuring generalisation error vs. number of samples, alignment and spike magnitude for different models in Fig. 3) appear to back up the authors’ central claims.
summary: This paper examines the performance of random feature models (RFMs) with anisotropic input data, characterised by spiked covariance. Whilst RFMs and linear models are asymptotically equivalent for isotropic data, RFMs tend to perform better in practical applications – an observation that the authors put down to the structured, anisotropic nature of real world data. They establish a ‘universality theorem’, extending the work of Hu and Lu (2023) to prove that there exists two different activation functions with equivalent performances if the first two statistical moments of $\sigma(\mathbf{F} \mathbf{x},y)$ match. They use this result to show that higher-order polynomial models are equivalent to RFMs for spiked covariance data, and show that in particular limits for this data linear models become sufficient (roughly, when input correlations and covariance spiking becomes small). Input-label correlation is found to be crucial in determining whether RFMs outperform linear models. The authors provide detailed theoretical guarantees and experimental validation on synthetic data.
weaknesses: 1. *Applicability to real-world data.* I appreciate that the specific parameterisation of Eq. 6 is needed for analytic tractability, but can the authors provide evidence that this kind of anisotropy in real-world datasets is responsible for the performance gap between RFMs and linear models? I think that the detailed analysis on this toy model is valuable in its own right, and would be of especial interest for readers of e.g. IEEE transactions on information theory (where Hu 2023, which this work builds on, was published). But to be suitable for ICLR, I think it would be helpful to provide some evidence that similar behaviour holds on real-world data. As a straightforward example, the authors could show performance on a real dataset with the polynomial kernel, ablating over different $l$ and inspecting when performance approaches RFM. Do any properties of real-world distributions roughly corresponds to $\theta$ and $\alpha$, spike magnitude and alignment? Can we see a difference in the gap for datasets where these properties differ? Another way to better orient the work in existing ML literature would be to go further explaining the relationship to double descent, beyond the brief empirical observations around line 490, or the relationship (if any) to feature learning in neural networks. 
2. *Other points.*
- *Thm 1*. An activation function $\sigma$ and its Hermite expansion $\widehat{\sigma}$ will perform identically if we include every term since the Hermite basis is complete (line 280). Under suitable asymptotic conditions (needed to bound the contribution from higher order terms), it’s not too surprising that it becomes sufficient to consider just the first two terms. This appears to be the central conclusion of Thm. 1. I don’t doubt that the authors’ derivations are correct and can see that this is an essential corollary for the arguments that follow, but I wonder whether emphasising this as a `key advancement of this work’ (line 263) is overstated/misplaces emphasis.
- *Assumption 7*. A.7, which is required to keep the perturbed objective convex (line 742), seems to be essential to Thm. 1. Is there any evidence that this will hold in practice? This seems like an important and nontrivial assumption, needed to make claims about the generalisation error, so I wonder why it is deferred to the appendix. Perhaps it would be better to include it in the main text, and probe whether it holds in practice in the experiments.  
- *Odd activation functions*. Assumption A.6. requires that activation functions are odd. The authors empirically show that their arguments hold even when this is relaxed, e.g. for ReLU. Can the authors give some idea of where their arguments break down for non-odd activation functions (I guess around line 1006)? And why doesn’t Fig. 3 include any odd activations, e.g. tanh, given that the proof only holds in this case?

As a stylistic comment, Fig. 2a is difficult to read since the noisy linear and random feature models for $\tanh$ are overlaid. Perhaps this could be presented a little differently to make things clearer.

Meta-review: (a) Summary:

The paper explores the conditions under which Random Feature Models (RFMs) outperform linear models in scenarios with spiked covariance data and strong input-label correlation. The authors demonstrate that RFMs can exceed linear models by leveraging this correlation and establish an equivalence between RFMs and noisy polynomial models, with polynomial degree determined by the input-label correlation. They present a universality theorem to support their claims and validate their results with numerical experiments.

(b) Strengths:

The paper addresses an important and timely topic in machine learning theory, providing insights into the advantages of RFMs under structured data conditions.
It introduces a new perspective on how input-label correlation influences model performance, which is a valuable addition to the existing literature.
The use of spiked covariance as a model enriches the theoretical analysis and broadens applicability to realistic datasets.
Numerical simulations are consistent with the theoretical findings.

(c) Weaknesses:

Novelty: Despite the efforts to differentiate the work, the contribution feels incremental relative to prior studies, such as those by Ba et al. (2023) and Mousavi-Hosseini et al. (2023), which addressed similar phenomena in kernel regression contexts. The parallels to existing studies raise questions about the novelty of the results.
Practical Relevance: While theoretical results are well-grounded, the assumptions about spiked covariance data and strong alignment between input and target seem restrictive and may limit real-world applicability.
Empirical Evaluation: The experimental validation lacks breadth, relying heavily on synthetic data and a single dataset (CIFAR-10 binary classification). This undermines the claim of practical relevance.
Clarity and Scope: Reviewers noted room for improvement in presentation, especially in defining key assumptions and notations, and explaining connections to related works.

(d) Decision:

The decision is to reject the paper. While the paper contributes to understanding RFMs, the incremental nature of the contribution and the limited empirical validation weigh against its acceptance. I encourage the authors to address these issues in a future submission by expanding empirical validation to real-world datasets, providing stronger distinctions from related works, and improving the clarity of assumptions and their practical implications.

Decision: Reject