<figure><img src="image_1.png" /></figure>

HAL
open science

# Generalization in Representation Models via Random Matrix Theory: Application to Recurrent Networks

Yessin Moakher, Malik Tiomoko, Cosme Louart, Zhenyu Liao

► To cite this version:

Yessin Moakher, Malik Tiomoko, Cosme Louart, Zhenyu Liao. Generalization in Representation Models via Random Matrix Theory: Application to Recurrent Networks. 2025. hal-05340623

HAL Id: hal-05340623

https://inria.hal.science/hal-05340623v1

Preprint submitted on 3 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Generalization in Representation Models via Random Matrix Theory: Application to Recurrent Networks

Yessin Moakher<fnref n="1" />,<fnref n="2" />, Malik Tiomoko<fnref n="1" />, Cosme Louart<fnref n="3" />, and Zhenyu Liao<fnref n="4" />

## Abstract

We first study the generalization error of models that use a fixed feature representation (frozen intermediate layers) followed by a trainable readout layer. This setting encompasses a range of architectures, from deep random-feature models to echo-state networks (ESNs) with recurrent dynamics. Working in the high-dimensional regime, we apply Random Matrix Theory to derive a closed-form expression for the asymptotic generalization error. We then apply this analysis to recurrent representations and obtain concise formula that characterize their performance. *Surprisingly*, we show that a linear ESN is equivalent to ridge regression with an exponentially time-weighted (“memory”) input covariance, revealing a clear inductive bias toward recent inputs. Experiments match predictions: ESNs win in low-sample, short-memory regimes, while ridge prevails with more data or long-range dependencies. Our methodology provides a general framework for analyzing overparameterized models and offers insights into the behavior of deep learning networks.

## 1 Introduction

Deep learning has achieved remarkable empirical success across a wide range of applications. Despite their extreme overparameterization, modern neural networks often generalize well, a phenomenon that challenges classical statistical learning theory. In particular, recent observations of double descent behavior reveal that conventional machine learning arguments fail to capture the complexities of high-dimensional learning dynamics. This has motivated the development of theoretical frameworks aimed at understanding when and why overparameterized models generalize.

A promising avenue in this direction is provided by random features models, which were originally introduced as scalable approximations to kernel methods [20]. Over time, these models have also been studied as surrogates for neural networks [18], offering a simplified yet powerful framework for understanding generalization. In these models, input vectors are mapped through a random nonlinear transformation, after which only a linear readout layer is trained. This decoupling of representation and readout retains much of the expressive power of neural networks while allowing for deeper theoretical analysis.

Building on this line of work, we focus on the broader question of learning under fixed representations. In this setting, the representation function is fixed a priori, and learning occurs exclusively through the readout. This setting arises in various contexts, including random-feature models [20], random intermediate layers in deep networks [23], and reservoir computing architectures such as

<footnote n="1">¹Huawei Noah’s Ark Lab, Huawei Technologies, Paris, France</footnote>

<footnote n="2">²École Polytechnique, France</footnote>

<footnote n="3">³Chinese University of Hong Kong, Shenzhen, China</footnote>

<footnote n="4">⁴Huazhong University of Science and Technology, Wuhan, China</footnote>
*Echo State Networks* (ESNs) [11]. Its popularity stems from the practical advantages it offers for analysis.

However, despite its prevalence, the generalization properties of such models particularly when the feature map is structured rather than i.i.d. are still not fully understood. In the recurrent setting, [6] analyzed ESNs trained and later tested on a single time series, assuming both the input and target to be independent and deterministic, and derived closed-form expressions for the asymptotic mean-square error. While this provided valuable insights, the restrictive assumptions of independence and determinism limit its applicability.

In recent years, Random Matrix Theory (RMT) has emerged as a central tool for analyzing learning algorithms in the high-dimensional regime, where the number of samples and the feature dimension grow proportionally. RMT provides precise asymptotic predictions for performance metrics such as training and test error [5] and has been also applied to estimation problems, including for multi-task regression [10].

**Our approach.** We address these gaps by developing a unified RMT framework for models with arbitrary fixed feature representations, including structured, recurrent, linear and non linear maps, under the assumption that the transformed features form a concentrated random vector [13] and does not require i.i.d. or full-rank projections. Unlike most prior work on static models, we tackle the recurrent setting, where the representation arises from a neural network with temporal dynamics.

**Contributions.** Our main contributions are:

1. **General risk characterization (Theorem 1).** We derive a closed-form asymptotic limit for the out-of-sample risk of ridge regression on any fixed representation under the concentrated random vector assumption, including second-order deterministic equivalents obtained via concentration-of-measure arguments.

2. **Specialization to ESNs (Theorem 2).** We derive a closed-form performance characterization for Linear Echo State Recurrent Networks.

3. **Insights.** We show that the test risk of a linear ESN is equivalent to ridge regression on an exponentially time-weighted ("memory") input covariance and derive intuition on why we don't observe double descent with Linear ESN.

## 2 Related Work

Our study is connected to the expanding literature that uses random matrix theory (RMT) to analyze machine learning models in high-dimensional settings. Variations across studies mainly stem from differences in the assumed distribution of the features X and from the relationship between X and the target Y (often referred to as the teacher-student framework).

Early analyses often assumed that Y was deterministic. For instance, [6] examined the dynamics of a zero-shot linear echo state network (ESN) under this setting. Similarly, [17] investigated a Gram random matrix student model with predictions of the form $\hat{Y} = \sigma(WX)$ to study random neural networks, also assuming deterministic targets Y. In the same spirit, [19] studied the Gram matrix $\sigma(WX)$ when both the data X and targets Y were independent Gaussian variables.

Subsequent work shifted towards settings where Y is linearly dependent on X, such as

$$Y = \theta^* X + \epsilon.$$
In the case of linear student predictors $\hat{Y} = \hat{\theta}X$ corresponding to ridge regression models a variety of theoretical analyses have been carried out, highlighting phenomena such as the double descent of the test error curve. Early results often assumed isotropic feature distributions. For example, [8] studied ridge regression when the inputs $x_i$ were sampled from an isotropic Gaussian distribution $x_i \sim N(0, I_d)$. Later works, such as [27] and [21], extended these results to more general covariance structures. [9] derived insights on the ridgeless least squares interpolation.

Another closely related line of work considers random projections prior to learning. For instance, [4] analyzed models where the training data are transformed via random matrices, considering settings of the form $WSX$ where $S$ is a random projection matrix. Such formulations naturally connect to the case of linear ESNs, in which the reservoir acts as a fixed, structured projection of the input. However, most existing analyses assume i.i.d. entries in the projection matrix and often full-rank transformations.

[18] adds a non linear component the teacher model and perform ridge regression on random features. More recently, [3] study the problem of learning a polynomial target function when data is provided with a spiked covariance structure $x_i \sim N(0, I_d + \theta\mu\mu^\top)$.

**Notation** Throughout the paper, we use capital letters to denote matrices, lowercase letters for scalars, and lowercase bold letters for vectors. The Frobenius norm is denoted by $\|\cdot\|_F$, the operator norm by $\|\cdot\|$, and the Euclidean norm by $\|\cdot\|_2$. We say that $u = O(v)$ if the ratio $u/v$ remains bounded in this limit.

**Organization** The remainder of the paper is organized as follows. Section 3 introduces the problem setting and defines the teacher and student models. Section 4 states our main assumptions and theoretical risk characterization, including comparisons between models. Section 5 presents experiments supporting our theory, and Section 6 concludes.

# 3 Problem Setting

We consider a supervised learning task in which the goal is to learn a predictor that maps an input sequence to an output vector, based on a finite set of training samples.

**Training and test setup.** Let $\{(u_i, y_i)\}_{i=1}^N$ be $N$ independent and identically distributed (i.i.d.) training pairs, where

$$u_i \in \mathbb{R}^T, \quad y_i \in \mathbb{R}^q.$$

Throughout, both the training and test pairs are assumed to be generated from the following model.

**Definition 1** (Noisy linear model). An input-output pair $(u, y) \in \mathbb{R}^T \times \mathbb{R}^q$ is said to follow a noisy linear model if

$$y = \Theta_{*}^{\top} u + \epsilon, \qquad (1)$$

where:

* $\Theta_* \in \mathbb{R}^{T \times q}$ is the (unknown) ground-truth parameter matrix,

* $\epsilon \in \mathbb{R}^q$ is a noise vector with i.i.d. entries of zero mean and variance $\sigma^2$, and $\epsilon$ is independent of $u$.

This model can be viewed as a linearization of more general nonlinear models in high dimensions (see, e.g., [19], [7]).
**Feature representation.** Rather than using the raw input **u** directly, we first transform it via a fixed (possibly nonlinear) representation map:

$$F: \mathbb{R}^T \to \mathbb{R}^n, \quad \mathbf{z} = F(\mathbf{u}). \qquad (2)$$

Examples of such $F$ include:

* the reservoir state of a (linear or nonlinear) Echo State Network (ESN),

* random feature maps,

* intermediate layers of a pretrained network.

In this work, $F$ is fixed and known; only the final linear readout is learned from data.

**Linear readout with ridge regularization.** Let $Z = [\mathbf{z}_1, \dots, \mathbf{z}_N] \in \mathbb{R}^{n \times N}$ be the matrix of feature vectors and $Y = [\mathbf{y}_1, \dots, \mathbf{y}_N] \in \mathbb{R}^{q \times N}$ the matrix of corresponding targets. We estimate the output weights via ridge-regularized regression:

$$\hat{W}_{\text{out}} := \arg \min_{W \in \mathbb{R}^{q \times n}} \frac{1}{N} \sum_{i=1}^{N} \|\mathbf{y}_i - W\mathbf{z}_i\|_2^2 + \lambda \|W\|_F^2 = \frac{1}{N} YZ^{\top} \left( \frac{1}{N}ZZ^{\top} + \lambda I_n \right)^{-1}, \quad (3)$$

where $\lambda > 0$ is the regularization parameter.

**Prediction.** Given a new test input $\mathbf{u}' \in \mathbb{R}^T$, we compute its feature vector $\mathbf{z}' := F(\mathbf{u}')$ and output

$$\hat{\mathbf{y}}' = \hat{W}_{\text{out}} \mathbf{z}'. \qquad (4)$$

**Definition 2 (Out-of-sample risk).** The out-of-sample risk of the predictor is the mean squared prediction error on an independent test sample $(\mathbf{u}', \mathbf{y}')$:

$$\mathcal{R} := \frac{1}{q} \mathbb{E} \left[ \|\mathbf{y}' - \hat{\mathbf{y}}'\|_2^2 \right], \qquad (5)$$

where the expectation is taken over both training and test data.

# 4 Main Technical Results

## 4.1 Asymptotic characterization of out-of-sample risk

We define the following quantities: $\Sigma_u := \mathbb{E}[\mathbf{u}\mathbf{u}^\top] = \text{Cov}(\mathbf{u}) + \mathbb{E}(\mathbf{u})\mathbb{E}(\mathbf{u})^\top \in \mathbb{R}^{T \times T}$, $\Sigma_z := \mathbb{E}[\mathbf{z}\mathbf{z}^\top] \in \mathbb{R}^{n \times n}$, $\Sigma_{uz} := \mathbb{E}[\mathbf{u}\mathbf{z}^\top] \in \mathbb{R}^{T \times n}$ and the resolvent $Q := (\frac{ZZ^\top}{N} + \lambda I_n)^{-1}$.

**Assumptions.** In order to use Random Matrix Theory (RMT) tools, we make assumptions on the data distribution and the asymptotic regime.

**Definition 3 (Concentrated random vector).** A random vector $\mathbf{x} \in \mathbb{R}^n$ is said to be concentrated if there exist constants $C, c > 0$, independent of $n$, such that: for every 1-Lipschitz function $f: \mathbb{R}^n \to \mathbb{R}$ and all $t \ge 0$,

$$\mathbb{P}(|f(\mathbf{x}) - \mathbb{E}[f(\mathbf{x})]| \ge t) \le C e^{-ct^2}.$$
**Assumption 1.** The representation vector $\mathbf{z} \in \mathbb{R}^T$ is a concentrated random vector in the sense of the definition 3 and satisfies $\|\mathbb{E}[\mathbf{z}]\|_2 = O(1)$.

This class includes Gaussian vectors with covariance matrices bounded in operator norm, uniform vectors on the sphere, and any Lipschitz transformation thereof (e.g., features from GANs [24]). This assumption allows us to apply the Hanson-Wright inequality in the context of random matrix theory with non-isotropic vectors [1]. Other works, (e.g. [4]) make the assumption $z = \Sigma^{1/2}S$, where $S$ is i.i.d. subgaussian and $\Sigma$ bounded, in order to apply the Hanson-Wright inequality in the isotropic case [22]. The two setups overlap (e.g., Gaussians), but neither contains the other in full generality, and in most cases the same results could be derived under either assumption.

**Assumption 2.** We work in the classical random matrix theory proportional asymptotics regime, where the number of reservoirs $n$ and the number of samples $N$ diverge proportionally. That is, $\frac{n}{N} \to \gamma \in (0, \infty)$, as $n, N \to \infty$.

**Asymptotic Expressions.** We will make use of the following quantities:

$$ \bar{Q} := \left( \frac{\Sigma_z}{1 + \delta} + \lambda I_n \right)^{-1}, \quad \delta := \frac{1}{N} \mathrm{Tr}(\Sigma_z \bar{Q}), \quad \alpha := \frac{1}{N} \mathrm{Tr} \left( \frac{\Sigma_z}{1 + \delta} \bar{Q} \frac{\Sigma_z}{1 + \delta} \bar{Q} \right) $$

**Theorem 1 (Fixed Representation Generalization).** Let $\mathbf{u} \in \mathbb{R}^T$ be an input vector, and let $\mathbf{z} = F(\mathbf{u}) \in \mathbb{R}^n$ be a representation vector obtained through a transformation of $\mathbf{u}$.

Under the linear model of Definition 1, where predictions are defined as in Equation 4, and under Assumptions 1 and 2, the following expressions hold: $\mathcal{R} = \mathcal{B}^2 + \mathcal{V} + \sigma^2$ where

$$ \mathcal{B}^2 \to \frac{1}{1-\alpha} \left( \mathrm{Tr}(\Theta_*^\top \Sigma_u \Theta_*) - \frac{2}{1+\delta} \mathrm{Tr} \left( \Theta_*^\top \Sigma_{uz} \bar{Q} \Sigma_{uz}^\top \Theta_* \right) + \frac{1}{(1+\delta)^2} \mathrm{Tr} \left( \Theta_*^\top \Sigma_{uz} \bar{Q} \Sigma_z \bar{Q} \Sigma_{uz}^\top \Theta_* \right) \right) \\ \mathcal{V} \to \sigma^2 \frac{\alpha}{1-\alpha} $$

*Proof.* See Appendix A. The argument combines Sherman-Morrison identities to disentangle the resolvent $Q$ from the other random variables, concentration of quadratic forms to replace terms by their expectations, and deterministic equivalents for the limit. □

**Remark 1 (Special Cases).** Under the assumption that $\mathbf{u}$ is concentrated, note that for $F$ equal to the identity we recover ridge regression. For $F$ given by an ESN or even a feedforward neural network, under the assumptions that the weight matrices are normalized, and activation functions $f$ are Lipschitz with $f(0) = 0$, we have that $\mathbf{z}$ is concentrated with $\|\mathbb{E}[\mathbf{z}]\|_2 = O(1)$. The same holds for random projections with bounded operator norm, not necessarily i.i.d. sub-Gaussian.

## 4.2 Application to Recurrent Models Representations

The goal of this section is to apply Theorem 1 to the case of recurrently generated representations. To this end, we consider an Echo State Network (ESN), a class of recurrent neural networks (RNNs) designed for sequential data processing, particularly in time series forecasting, speech recognition, and dynamical system modeling.

**Definition 4 (Echo State Network).** An Echo State Network consists of:

1. A fixed input layer $w_{in} \in \mathbb{R}^n$ that maps the input $u_i(t) \in \mathbb{R}$ into an $n$-dimensional space.
2. A fixed recurrent reservoir layer $W \in \mathbb{R}^{n \times n}$ that captures the temporal dependency of data.

3. A trainable output layer $W_{\text{out}} \in \mathbb{R}^{q \times n}$ that maps reservoir states to predictions.

For an input $u_i(t)$, the reservoir state of the ESN is denoted $\mathbf{x}_i(t) \in \mathbb{R}^n$ and evolves according to:

$$ \mathbf{x}_i(t) = f(u_i(t)\mathbf{w}_{\text{in}} + W\mathbf{x}_i(t-1)), \qquad (6) $$

where $\mathbf{x}_i(0) = 0$, and $f(\cdot)$ is the activation function such as $f(t) = \tanh(t)$ or $f(t) = \text{ReLU}(t) = \max(t, 0)$ that applies entry-wise. The representation vector $\mathbf{z}_i \in \mathbb{R}^n$ is defined as the reservoir state at the (final) time step $T$: $\mathbf{z}_i = \mathbf{x}_i(T)$.

Compared to standard RNNs, ESNs fix the input and recurrent weights (typically drawn randomly) and train only the output layer.

**Linear ESN** For linear ESN (i.e $f$ is the identity), we have a closed form for the relationship between $U$ and $Z$, that is:

$$ Z = SU, \quad \text{where } S = [W^{T-1}\mathbf{w}_{\text{in}}, W^{T-2}\mathbf{w}_{\text{in}}, \dots, W^0\mathbf{w}_{\text{in}}] \in \mathbb{R}^{n \times T} \qquad (7) $$

The matrix $U$ is being projected by a sort of a Kalman controllability matrix $S$. In that case we have a closed form expression of $\Sigma_{uz} = \Sigma S^\top$ and $\Sigma_z = S\Sigma S^\top$.

**Asymptotic Expressions.** To ensure dynamical stability and guarantee that $z$ remains a concentrated vector under the assumption that $u$ is concentrated, we impose the following assumption.

**Assumption 3 (Linear ESN setting).** We consider a linear ESN with recurrent matrix $W = \frac{W_0}{\varphi \rho(W_0)}$, where $W_0 \in \mathbb{R}^{n \times n}$ has i.i.d. entries $W_{0,ij} \sim \mathcal{N}(0, 1)$, $\rho(W_0)$ is its spectral radius, and $\varphi > 0$ controls the effective spectral radius. The input weight vector $\mathbf{w}_{\text{in}}$ has i.i.d. entries normalized $\mathcal{N}(0, \frac{1}{n})$.

The parameter $\varphi < 1$ plays a key role in the dynamical stability of the ESN. Larger $\varphi$ can capture longer temporal dependencies but may also lead to instability if the system enters a regime of diverging activations.

**Theorem 2 (Linear ESN Generalization).** Let $(\mu_i, v_i)$ be the eigenvalue-eigenvector pairs of $\Sigma_u^{1/2} \text{diag}(\varphi^{i-T})_{1 \le i \le T} \Sigma_u^{1/2}$. Under the assumption that $u \in \mathbb{R}^T$ is concentrated, and under Assumptions 2 and 3, the out-of-sample risk decomposes as $\mathcal{R} = \mathcal{B}^2 + \mathcal{V} + \sigma^2$, where

$$ \mathbb{E}_{W, \mathbf{w}_{\text{in}}}[\mathcal{B}^2] \rightarrow \frac{1}{1-\alpha} \sum_{i=1}^{T} \frac{\kappa^2}{(\mu_i + \kappa)^2} \|\Theta_*^\top \Sigma_u^{1/2} v_i\|_2^2, \quad \mathbb{E}_{W, \mathbf{w}_{\text{in}}}[\mathcal{V}] \rightarrow \sigma^2 \frac{\alpha}{1-\alpha}, $$

with $\kappa := \lambda(1 + \delta)$ and $\alpha := \sum_{i=1}^{T} \frac{\mu_i}{N(\mu_i + \kappa)^2}$.

*Proof.* See Appendix B. We first show that $S^\top S$ is concentrated, and then compute the limit $M_\infty$ of its expectation using random matrix theory. Next, we prove that $\mathcal{R}$, viewed as a function of $S^\top S$, is Lipschitz with constant of order $O(1)$. This allows us to replace, in the limit, $\mathbb{E}[\mathcal{R}(S^\top S)]$ with $\mathcal{R}(M_\infty)$. $\square$

**Remark 2 (Difference with Ridge regression and Interpretation).** In the case of ridge regression ($\Sigma_u = \Sigma_{uz} = \Sigma_z$), we obtain the same expression but the same expression with $(\mu_i, v_i)$ the eigenvalue-eigenvector pairs of $\Sigma_u$ instead of $\Sigma_u^{1/2} \text{diag}(\varphi^{i-T})_{1 \le i \le T} \Sigma_u^{1/2}$. Which proves that Linear ESN is equivalent to ridge regression but with exponentially time-weighted version covariance, which down-weights older inputs. This induces a short-memory bias, distorting the input statistics and discarding part of the long-term information.
**Remark 3 (Optimal Regularization).** The asymptotic risk depends on the scalar $\kappa = \lambda(1+\delta)$, which can be optimized as a real-valued function independently of the fixed point $\delta$ (which itself depends on $\lambda$). Once the optimal $\kappa$ is found, the corresponding $\delta$ can be recovered from $\frac{\delta}{1+\delta} = \frac{1}{N} \sum_{i=1}^{T} \frac{\mu_i}{\mu_i+\kappa}$. This then yields the optimal $\lambda$.

In general, there is no closed-form expression for $\lambda^*$ for arbitrary $\Sigma_u$. However, when $\Sigma_u = I_T$, we obtain $\lambda^* = \frac{T}{N} \cdot \text{SNR}$, $\text{SNR} = \frac{\sum_{i=1}^{T} \|\Theta^{*\top} \mathbf{v}_i\|^2}{\sigma^2}$.

# 5 Experiments: Comparing Ridge Regression and Linear ESN

## 5.1 Double Descent

The double descent phenomenon describes how the test error decreases, peaks near the interpolation threshold, and then decreases again as model complexity grows. This behavior is well established for linear predictors, both empirically and theoretically.

**Key mechanism.** The factor $(1-\alpha)^{-1}$ in Theorem 1 diverges as $\alpha$ tends to 1. Recall that

$$ \alpha = \frac{1}{N} \sum_{i=1}^{T} \frac{\mu_i^2}{(\mu_i + \lambda(1 + \delta))^2}. $$

Double descent is typically observed for $\lambda$ close to 0 [18]. In this case,

$$ \alpha \xrightarrow{\lambda \to 0} \frac{1}{N} \sum_{i=1}^{r} 1 = \frac{r}{N}, $$

where $r$ is the rank of $\Sigma_u^{1/2} \text{diag}(\varphi^{i-T})_{1 \le i \le T} \Sigma_u^{1/2}$. In classical ridge regression with $\varphi = 1$, the feature covariance matrix is full rank, and double descent occurs when $N = T$. In contrast, for a Linear ESN, the matrix $\Sigma_u^{1/2} \text{diag}(\varphi^{i-T})_{1 \le i \le T} \Sigma_u^{1/2}$ is low rank, since $\varphi^{i-T}$ rapidly decays to 0 in numerical computations, and thus double descent is not observed as shown in Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: Comparison of test error curves for Linear ESNs and random projections as a function of $n/N$. For random projections, the test error exhibits a clear double descent peak near the interpolation threshold ($n/N \approx 1$). In contrast, Linear ESNs show no such peak due to their effective low-rank feature covariance, as explained above.</figcaption></figure>
## 5.2 Comparing ESNs and Ridge Regression

The Remark 2 indicates that Echo State Networks (ESNs) could surpass Ridge regression in scenarios characterized by *limited data and short temporal dependencies*. This advantage arises from the inherent inductive bias of ESNs: they assume that the relevant information is encoded in the recent history of the input, effectively emphasizing short-term temporal correlations.

Ridge regression, when applied directly to raw input vectors or generic features, lacks this temporal structure. In low-data regimes, this absence of an inductive bias often leads to overfitting or suboptimal generalization. ESNs, in contrast, implicitly filter past inputs through their recurrent reservoir, efficiently capturing patterns from the recent past and making better use of scarce data.

<figure><img src="image_4.png" /><figcaption>Figure 2: Comparison between ESNs and ridge regression. The test error is shown as a function of the sample size *N* (y-axis) and the signal decay rate *ρ* (x-axis), where $\theta^* = (\rho^t, t \in \{1, \dots, T\})$ for both models. Shaded regions indicate the standard deviation over multiple runs. The dashed line marks the separation threshold. ESNs exhibit superior performance in the limited-data, short-memory regime.</figcaption></figure>

In Figure 2, we confirm this intuition: ESNs consistently outperform Ridge regression when $T/N$ is small and when target function depends predominantly on short-term dependencies (small signal decay rate $\rho$). In this regime, the exponential time-weighting of ESNs effectively acts as a regularizer, improving generalization by attenuating the influence of distant, noisy inputs. However, as $N$ grows, Ridge regression eventually surpasses ESNs, since, as discussed in Remark 2, the ESNs inevitably discards part of the long-term information, which becomes increasingly useful when sufficient samples are available.

# 6 Conclusion

We introduced a general high-dimensional theory for ridge regression on arbitrary fixed representations, unifying and extending prior analyses of random features, pretrained models, and recurrent architectures. Applied to Echo State Networks, our framework yields simple closed-form predictions that match experiments, reveal their temporal inductive bias, and explain their absence of double descent. Our results show that Echo State Networks can outperform ridge regression in low-data, short-memory regimes due to their built-in temporal inductive bias, while ridge becomes advantageous as the sample size grows.
## References

[1] Radosław Adamczak. A note on the hanson-wright inequality for random vectors with dependencies. *Electronic Communications in Probability*, 2014.

[2] Johannes Alt, László Erdős, and Torben Krüger. Spectral radius of random matrices with independent entries. *Probability and Mathematical Physics*, 2021.

[3] Jimmy Ba, Murat A Erdogdu, Taiji Suzuki, Zhichao Wang, and Denny Wu. Learning in the presence of low-dimensional structure: A spiked random matrix perspective. In *Advances in Neural Information Processing Systems*, 2023.

[4] Francis Bach. High-dimensional analysis of double descent for linear regression with random projections. *SIAM Journal on Mathematics of Data Science*, pages 26-50, 2024.

[5] Romain Couillet and Zhenyu Liao. *Random matrix methods for machine learning*. Cambridge University Press, 2022.

[6] Romain Couillet, Gilles Wainrib, Hafiz Tiomoko Ali, and Harry Sevi. A random matrix approach to echo-state neural networks. In *Proceedings of The 33rd International Conference on Machine Learning*, 2016.

[7] Yatin Dandi, Ludovic Stephan, Florent Krzakala, Bruno Loureiro, and Lenka Zdeborová. Universality laws for gaussian mixtures in generalized linear models. In *Advances in Neural Information Processing Systems*, 2023.

[8] Lee H. Dicker. Ridge regression and asymptotic minimax estimation over spheres of growing dimension. *Bernoulli*, 2016.

[9] Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J. Tibshirani. Surprises in high-dimensional ridgeless least squares interpolation. *The Annals of Statistics*, 2022.

[10] Romain Ilbert, Malik Tiomoko, Cosme Louart, Ambroise Odonnat, Vasilii Feofanov, Themis Palpanas, and Ievgen Redko. Analysing multi-task regression via random matrix theory with application to time series forecasting. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2024.

[11] Herbert Jaeger. The "echo state" approach to analysing and training recurrent neural networks- with an erratum note'. *Bonn, Germany: German National Research Center for Information Technology GMD Technical Report*, 2001.

[12] Todd Kemp, Karl Mahlburg, Amarpreet Rattan, and Clifford Smyth. Enumeration of non-crossing pairings on bit strings, 2009.

[13] Michel Ledoux. *The Concentration of Measure Phenomenon*. American Mathematical Society, 2001.

[14] Cosme Louart. *Random matrix theory and concentration of the measure theory for the study of high dimension data processing*. PhD thesis, Université Grenoble Alpes [2020-....], 2023.

[15] Cosme Louart and Romain Couillet. Concentration of measure and large random matrices with an application to sample covariance matrices, 2021.
[16] Cosme Louart and Romain Couillet. A concentration of measure and random matrix approach to large-dimensional robust statistics. *The Annals of Applied Probability*, 2022.

[17] Cosme Louart, Zhenyu Liao, and Romain Couillet. A random matrix approach to neural networks. *The Annals of Applied Probability*, 28(2):1190–1248, 2018.

[18] Song Mei and Andrea Montanari. The generalization error of random features regression: Precise asymptotics and double descent curve, 2020.

[19] Jeffrey Pennington and Pratik Worah. Nonlinear random matrix theory for deep learning. In *Advances in Neural Information Processing Systems*, 2017.

[20] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In *Advances in Neural Information Processing Systems*, 2007.

[21] Dominic Richards, Jaouad Mourtada, and Lorenzo Rosasco. Asymptotics of ridge(less) regression under general source condition. In *Proceedings of The 24th International Conference on Artificial Intelligence and Statistics*, 2021.

[22] Mark Rudelson and Roman Vershynin. Hanson-wright inequality and sub-gaussian concentration, 2013.

[23] Dominik Schröder, Hugo Cui, Daniil Dmitriev, and Bruno Loureiro. Deterministic equivalent and error universality of deep random features learning. In *International Conference on Machine Learning*. PMLR, 2023.

[24] Mohamed El Amine Seddik, Cosme Louart, Mohamed Tamaazousti, and Romain Couillet. Random matrix theory proves that deep learning representations of gan-data behave as gaussian mixtures. In *International Conference on Machine Learning*, 2020.

[25] Terence Tao. *Topics in random matrix theory*. American Mathematical Society, 2011.

[26] Roman Vershynin. *High-Dimensional Probability: An Introduction with Applications in Data Science*. Cambridge University Press, 2018.

[27] Denny Wu and Ji Xu. On the optimal weighted $\ell_2$ regularization in overparameterized linear regression. In *Proceedings of the 34th International Conference on Neural Information Processing Systems*, 2020.
# A Proof of theorem 1

In this section, we provide the proof of Theorem 1.

## A.1 Decomposition

Recall from Equation (1) that We consider the standard linear model:

$$Y = \Theta_{*}^{\top} U + E,$$

where $\Theta_{*} \in \mathbb{R}^{T \times q}$, $U := [\mathbf{u}_1, \dots, \mathbf{u}_N] \in \mathbb{R}^{T \times N}$, $Y := [\mathbf{y}_1, \dots, \mathbf{y}_N] \in \mathbb{R}^{q \times N}$, and $E \in \mathbb{R}^{q \times N}$ has i.i.d. entries with mean zero and variance $\sigma^2$. The noise $E$ is assumed to be independent of $U$.

We are interested in the generalization performance of Echo State Network (ESN), which we study through its out-of-sample risk. For a new, independent test pair $(\mathbf{u}', \mathbf{y}')$, where $\hat{\mathbf{y}}'$ is the prediction, the risk is defined as:

$$\mathcal{R} := \frac{1}{q} \mathbb{E} \left[ \left\| \mathbf{y}' - \hat{\mathbf{y}}' \right\|_2^2 \right].$$

where the expectation is taken on everything that is random (from train and test). Substituting the model $\mathbf{y}' = \Theta_{*}^{\top} \mathbf{u}' + \epsilon'$, we obtain:

$$\mathcal{R} = \frac{1}{q} \mathbb{E} \left[ \left\| \Theta_{*}^{\top} \mathbf{u}' - \hat{\mathbf{y}}' \right\|_2^2 \right] + \sigma^2.$$

We define the excess risk as $\tilde{\mathcal{R}} := \mathcal{R} - \sigma^2$. Letting $\Sigma_u := \mathbb{E}[\mathbf{u}\mathbf{u}^{\top}] \in \mathbb{R}^{T \times T}$, we have:

$$q\tilde{\mathcal{R}} = \mathrm{Tr}(\Theta_{*}^{\top}\Sigma_{u}\Theta_{*}) - 2\mathbb{E}\left[\mathrm{Tr}(\hat{\mathbf{y}}'\mathbf{u}^{\top}\Theta^{*})\right] + \mathbb{E}\left[\left\|\hat{\mathbf{y}}'\right\|_{2}^{2}\right].$$

In the ESN framework, predictions are obtained via:

$$\hat{\mathbf{y}}' = \hat{W}_{\mathrm{out}}\mathbf{z}',$$

where $\mathbf{z}' \in \mathbb{R}^n$ is the final reservoir state computed from $\mathbf{u}'$, and the output weights $\hat{W}_{\mathrm{out}} \in \mathbb{R}^{q \times n}$ are estimated using ridge regression:

$$\hat{W}_{\mathrm{out}} = \frac{1}{N} Y Z^{\top} \left( \frac{1}{N} Z Z^{\top} + \lambda I_n \right)^{-1},$$

with $Z = [\mathbf{z}_1, \dots, \mathbf{z}_N] \in \mathbb{R}^{n \times N}$ and $Y = [\mathbf{y}_1, \dots, \mathbf{y}_N] \in \mathbb{R}^{q \times N}$.

We define the resolvent $Q := \left( \frac{1}{N} Z Z^{\top} + \lambda I_n \right)^{-1}$. The excess risk becomes:

$$q\tilde{\mathcal{R}} = \mathrm{Tr}(\Theta_{*}^{\top}\Sigma_{u}\Theta_{*}) - 2\mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}YZ^{\top}Qz'\mathbf{u}'^{\top}\Theta^{*}\right)\right] + \mathbb{E}\left[\left\|\frac{1}{N}YZ^{\top}Qz'\right\|_{2}^{2}\right].$$

Substituting $Y = \Theta_{*}^{\top} U + E$ yields:

$$q\tilde{\mathcal{R}} = \mathrm{Tr}(\Theta_{*}^{\top}\Sigma_{u}\Theta_{*}) - 2\mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}(\Theta_{*}^{\top}U + E)Z^{\top}Qz'\mathbf{u}'^{\top}\Theta^{*}\right)\right] \\ + \mathbb{E}\left[\left\|\frac{1}{N}(\Theta_{*}^{\top}U + E)Z^{\top}Qz'\right\|_{2}^{2}\right].$$
We expand the second and third terms:

$$ q\tilde{\mathcal{R}} = \mathrm{Tr}(\Theta_{*}^{\top}\Sigma_{u}\Theta_{*}) \\ - 2\mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}\Theta_{*}^{\top}UZ^{\top}Qz'u'^{\top}\Theta^{*}\right)\right] - 2\mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}EZ^{\top}Qz'u'^{\top}\Theta^{*}\right)\right] \\ + \mathbb{E}\left[\left\|\frac{1}{N}\Theta^{*\top}UZ^{\top}Qz'\right\|_{2}^{2}\right] + \mathbb{E}\left[\left\|\frac{1}{N}EZ^{\top}Qz'\right\|_{2}^{2}\right] \\ + 2\mathbb{E}\left[\left\langle\frac{1}{N}\Theta^{*\top}UZ^{\top}Qz', \frac{1}{N}EZ^{\top}Qz'\right\rangle\right] $$

By independence of $E$ and $u'$, and zero-mean noise, the cross terms vanish in expectation:

$$ \mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}EZ^{\top}Qz'u'^{\top}\Theta^{*}\right)\right] = 0, \quad \mathbb{E}\left[\left\langle\frac{1}{N}\Theta_{*}^{\top}UZ^{\top}Qz', \frac{1}{N}EZ^{\top}Qz'\right\rangle\right] = 0. $$

Hence, the excess risk simplifies to:

$$ q\tilde{\mathcal{R}} = \mathrm{Tr}(\Theta^{*\top}\Sigma_{u}\Theta_{*}) - 2\mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}\Theta_{*}^{\top}UZ^{\top}Qz'u'^{\top}\Theta^{*}\right)\right] \\ + \mathbb{E}\left[\left\|\frac{1}{N}\Theta_{*}^{\top}UZ^{\top}Qz'\right\|_{2}^{2}\right] + \mathbb{E}\left[\left\|\frac{1}{N}EZ^{\top}Qz'\right\|_{2}^{2}\right]. $$

Let $\Sigma_{uz} := \mathbb{E}[uz^\top]$ and $\Sigma_z := \mathbb{E}[zz^\top]$. We get the decomposition of the excess risk:

$$ q\tilde{\mathcal{R}} = \underbrace{\mathrm{Tr}(\Theta_{*}^{\top}\Sigma_{u}\Theta_{*}) - 2\mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}\Theta^{*\top}UZ^{\top}Q\Sigma_{uz}^{\top}\Theta^{*}\right)\right] + \mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}\Theta^{*\top}UQ\Sigma_{z}Q\Theta^{*}\right)\right]}_{qB^2} \\ + q \underbrace{\mathbb{E}\left[\sigma^2\left(Q\frac{Z^{\top}Z}{N^2}Q\Sigma_z\right)\right]}_{\nu}. $$

## A.2 Asymptotic regime

### A.2.1 Bias term

**Second Term** We are interested in finding the limit of $\mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}\Theta^{*\top}UZ^{\top}Q\Sigma_{uz}^{\top}\Theta^{*}\right)\right]$, as both the number of reservoirs $n$ and the number of samples $N$ tend to infinity proportionally (or equivalently, to determine the deterministic equivalent of $\frac{1}{N}UZ^{\top}Q$, in the language of random matrix theory).

This proof follows the structure of arguments in random matrix theory. We provide the complete proof here and refer back to it for the other terms. We begin by disentangling the resolvent $Q$ from the other random variables using Sherman-Morrison's identities (Lemma 1). We then apply the concentration of quadratic forms to approximate them by their expectations. Finally, we use known deterministic equivalents for the resolvents to conclude.

**Step 1: Disentangling dependence** Writing $UZ^{\top} = \sum_{i=1}^{N} u_i z_i^{\top}$ we have:

$$ \mathbb{E}\left[\mathrm{Tr}\left(\frac{1}{N}\Theta^{*\top}UZ^{\top}Q\Sigma_{uz}^{\top}\Theta^{*}\right)\right] = \frac{1}{N}\sum_{i=1}^{N}\mathbb{E}\left[\mathrm{Tr}\left(\Theta_{*}^{\top}u_i z_i^{\top}Q\Sigma_{uz}^{\top}\Theta^{*}\right)\right] $$

The random variable $(u_i, z_i)$ depends on $Q$, so we use the following lemma to decouple their dependence.
**Lemma 1.** [Sherman-Morrison Identities] Let $Q_{-i} := (\frac{1}{N} \sum_{j=1, j \neq i}^{N} \mathbf{z}_j \mathbf{z}_j^\top + \lambda I_n)^{-1}$ denote the resolvent with the i-th vector $\mathbf{z}_i$ removed. Then, the following identities hold:

$$Q = Q_{-i} - \frac{1}{N} \cdot \frac{Q_{-i} \mathbf{z}_i \mathbf{z}_i^\top Q_{-i}}{1 + \frac{1}{N} \mathbf{z}_i^\top Q_{-i} \mathbf{z}_i} \\ Q \mathbf{z}_i = \frac{Q_{-i} \mathbf{z}_i}{1 + \frac{1}{N} \mathbf{z}_i^\top Q_{-i} \mathbf{z}_i}$$

We use the above equation 1:

$$\frac{1}{N} \sum_{i=1}^{N} \mathbb{E} \left[ \mathrm{Tr} \left( \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} Q \Sigma_{uz}^{\top} \Theta^{*} \right) \right] = \frac{1}{N} \sum_{i=1}^{N} \mathbb{E} \left[ \mathrm{Tr} \left( \Theta_{*}^{\top} \frac{\mathbf{u}_{i} \mathbf{z}_{i}^{\top} Q_{-i}}{1 + \frac{1}{N} \mathbf{z}_{i}^{\top} Q_{-i} \mathbf{z}_{i}} \Sigma_{uz}^{\top} \Theta^{*} \right) \right]$$

We now have that $\mathbf{u}_i$ and $\mathbf{z}_i$ are independent of $Q_{-i}$.

**Step 2: Concentration** The term $\frac{1}{N} \mathbf{z}_i^\top Q_{-i} \mathbf{z}_i$ is still random and depends on $\mathbf{z}_i$, but we will approximate it by its expectation using a concentration argument.

Before proceeding, let us recall some useful lemmas on concentrated random variables and sub-Gaussian random variables which will be used in the proofs.

**Definition 5** (Sub-Gaussian and Sub-exponential random variables). A centered real random variable $x$ is said to be $K$-sub-Gaussian if

$$\mathbb{P}(|x| \ge t) \le 2 \exp \left( - \frac{t^2}{2K^2} \right) \quad \forall t \ge 0.$$

This definition extends to a random vector $\mathbf{z} \in \mathbb{R}^n$ if, for all unit vector $\mathbf{u}, \mathbf{u}^\top \mathbf{z}$ is $K$-sub-Gaussian. A centered real random variable $x$ is called $(\nu, b)$-sub-exponential if

$$\mathbb{P}(|x| \ge t) \le 2 \exp \left[ - \min \left( \frac{t^2}{2\nu^2}, \frac{t}{2b} \right) \right] \quad \forall t \ge 0.$$

**Lemma 2** (Basic facts on sub-Gaussian and sub-exponential variables). Throughout, $\lesssim$ hides an absolute constant. We have the following results.

(i) $\mathbf{z} \in \mathbb{R}^n$ is concentrated $\implies$ $\mathbf{z}$ is sub-Gaussian with a constant independent of $n$.

(ii) **Moments of a sub-Gaussian vector.** If $\mathbf{z} \in \mathbb{R}^n$ is $K$-sub-Gaussian, then for every integer $k \ge 1$

$$\mathbb{E} \|\mathbf{z}\|_2^k \lesssim (K \sqrt{n})^k.$$

(iii) **Moments of a sub-exponential variable.** Suppose $x \in \mathbb{R}$ is $(\nu, b)$-sub-exponential with $K = \max(\nu, b)$. Then for every integer $k \ge 1$

$$\mathbb{E} |x|^k \lesssim (Kk\sqrt{n})^k.$$

(iv) **Maximum of sub-exponential variables.** Let $\{x_i\}_{i=1}^m$ be independent $(\nu, b)$-sub-exponential variables with $K = \max(\nu, b)$. Then for every integer $k \ge 1$

$$\mathbb{E} \left[ \max_{1 \le i \le m} |x_i|^k \right] \lesssim (Kk \log m)^k.$$
Proof. See [26].

□

**Lemma 3.** If $\mathbf{z} \in \mathbb{R}^n$ is concentrated (Definition 3) and $\|\mathbb{E}[\mathbf{z}]\|_2 \le M$, then $\Sigma_z := \mathbb{E}[\mathbf{z}\mathbf{z}^\top]$ satisfies

$$ \|\Sigma_z\| \le K_0 + M^2, $$

for an absolute constant $K_0 > 0$ depending only on the concentration constants (hence independent of $n$).

Proof. Write $\mu := \mathbb{E}[\mathbf{z}]$. For any $\mathbf{a} \in \mathbb{R}^n$ with $\|\mathbf{a}\|_2 = 1$,

$$ \mathbf{a}^\top \Sigma_z \mathbf{a} = \mathbb{E}[(\mathbf{a}^\top \mathbf{z})^2] = \mathrm{Var}(\mathbf{a}^\top \mathbf{z}) + (\mathbf{a}^\top \mu)^2. $$

Since $f(\mathbf{x}) = \mathbf{a}^\top \mathbf{x}$ is 1-Lipschitz, concentration (Definition 3) gives

$$ \mathbb{P}(|\mathbf{a}^\top \mathbf{z} - \mathbf{a}^\top \mu| \ge t) \le Ce^{-ct^2} \quad \forall t \ge 0. $$

Hence, by tail integration,

$$ \mathrm{Var}(\mathbf{a}^\top \mathbf{z}) = \mathbb{E}[(\mathbf{a}^\top \mathbf{z} - \mathbf{a}^\top \mu)^2] = \int_0^\infty 2t \mathbb{P}(|\mathbf{a}^\top \mathbf{z} - \mathbf{a}^\top \mu| \ge t) dt \le 2C \int_0^\infty t e^{-ct^2} dt = \frac{C}{c}. $$

Moreover $(\mathbf{a}^\top \mu)^2 \le \|\mu\|_2^2 \le M^2$. Therefore

$$ \mathbf{a}^\top \Sigma_z \mathbf{a} \le \frac{C}{c} + M^2 \quad \text{for all } \|\mathbf{a}\|_2 = 1, $$

so $\|\Sigma_z\| \le C/c + M^2$. Setting $K_0 := C/c$ completes the proof. □

□

**Lemma 4 (Concentration of Quadratic forms).** Given a fixed matrix $A$, and $\mathbf{z}$ a concentrated random vector (definition 3) with a bounded expectation $\mathbb{E}[\mathbf{z}]$, we have

$$ \mathbb{E}\left[\left(\mathbf{z}^\top A \mathbf{y} - \mathbb{E}[\mathbf{z}^\top A \mathbf{y}]\right)^k\right] \le C \|A\|_F^k $$

for some constant $C > 0$.

Proof. Recall the Hanson-Wright inequality for a concentrated random vector $z$ (see Prop B.41 [14]): there exist constants $C', c' > 0$ such that for all $t > 0$,

$$ \mathbb{P}\left(\left|\mathbf{z}^\top Q \mathbf{z} - \mathbb{E}[\mathbf{z}^\top Q \mathbf{z}]\right| \ge t\right) \le C' \left(\exp\left(-\frac{c' t^2}{\|Q\|_F^2}\right) + \exp\left(-\frac{c' t}{\|Q\|_2}\right)\right). $$

Therefore, the random variable $z^\top Q z - \mathbb{E}[z^\top Q z]$ is sub-exponential and from Lemma 2 we get the result. □

**Lemma 5 (Operator-norm moment bound).** Let $A$ be a $p \times q$ random matrix whose columns $A_i$ are independent, sub-Gaussian vectors in $\mathbb{R}^q$ with common mean $\mu$ and covariance matrix $\Sigma$. If $\|\mu\|_2 = O(\sqrt{p})$ and $\|\Sigma\| = O(1)$ then for every integer $k \ge 1$ we have

$$ \mathbb{E}[\|A\|^k] = O((\sqrt{p} + \sqrt{q})^k). $$
Proof. Let $A' := (A - \mathbf{1}_p \boldsymbol{\mu}^{\top}) \Sigma^{-1/2}$. Then the rows of $A'$ are i.i.d. mean-zero, isotropic sub-Gaussian vectors. By [26, Theorem. 4.6.1], for absolute $c_0$, $C_0 > 0$ and all $t \ge 0$,

$$ \mathbb{P}(\|A'\| > C_0(\sqrt{p} + \sqrt{q}) + t) \le 2e^{-c_0t^2}. $$

so $\|A'\|$ is sub-exponential with parameters $(\nu, b) = (C_1(\sqrt{p} + \sqrt{q}), C_2)$. The triangle inequality gives

$$ \|A\| \le \sqrt{p} \|\boldsymbol{\mu}\|_2 + \|\Sigma\|^{1/2} \|A'\|. $$

Hence $\|A\|$ is sub-exponential with parameters

$$ (\tilde{\nu}, \tilde{b}) = (\sqrt{p} \|\boldsymbol{\mu}\|_2 + \|\Sigma\|^{1/2} C_1(\sqrt{p} + \sqrt{q}), \|\Sigma\|^{1/2} C_2) = O(\sqrt{p} + \sqrt{q}). $$

Applying Lemma 2 we get the desired result. $\square$

Here we proof that we can replace the quadratic form by its expectation in our expression. Let $\tilde{Q}_i := \mathbb{E}(Q_{-i})$, $\delta_i := \frac{1}{N}\mathbf{z}_i^{\top}Q_{-i}\mathbf{z}_i$ and $\tilde{\delta}_i := \mathbb{E}(\frac{1}{N}\mathbf{z}_i^{\top}Q_{-i}\mathbf{z}_i) = \frac{1}{N}\text{Tr}(\tilde{Q}_i\Sigma_z)$.

$$ \begin{aligned} & \left| \frac{1}{N} \sum_{i=1}^{N} \mathbb{E} \left[ \text{Tr} \left( \Theta_{*}^{\top} \frac{\mathbf{u}_i \mathbf{z}_i^{\top} Q_{-i}}{1 + \tilde{\delta}_i} \Sigma_{uz}^{\top} \Theta^{*} \right) \right] - \frac{1}{N} \sum_{i=1}^{N} \mathbb{E} \left[ \text{Tr} \left( \Theta_{*}^{\top} \frac{\mathbf{u}_i \mathbf{z}_i^{\top} Q_{-i}}{1 + \tilde{\delta}_i} \Sigma_{uz}^{\top} \Theta^{*} \right) \right] \right| \\ &= \left| \frac{1}{N} \sum_{i=1}^{N} \mathbb{E} \left[ \text{Tr} \left( \Theta_{*}^{\top} \mathbf{u}_i \mathbf{z}_i^{\top} Q_{-i} \frac{\delta_i - \tilde{\delta}_i}{(1 + \tilde{\delta}_i)(1 + \tilde{\delta}_i)} \Sigma_{uz}^{\top} \Theta^{*} \right) \right] \right| \\ &= \left| \frac{1}{N} \sum_{i=1}^{N} \mathbb{E} \left[ \text{Tr} \left( \Theta_{*}^{\top} \mathbf{u}_i \mathbf{z}_i^{\top} Q \frac{\delta_i - \tilde{\delta}_i}{1 + \tilde{\delta}_i} \Sigma_{uz}^{\top} \Theta^{*} \right) \right] \right| \end{aligned} $$

Let $D := \text{diag}(\frac{\delta_i - \tilde{\delta}_i}{1 + \tilde{\delta}_i}, i \in \{1, \dots, N\})$, we can write back the sum into a matrix form:

$$ \left| \frac{1}{N} \sum_{i=1}^{N} \mathbb{E} \left[ \text{Tr} \left( \Theta_{*}^{\top} \mathbf{u}_i \mathbf{z}_i^{\top} Q \frac{\delta_i - \tilde{\delta}_i}{\tilde{\delta}_i} \Sigma_{uz}^{\top} \Theta^{*} \right) \right] \right| = \left| \frac{1}{N} \mathbb{E} \left[ \text{Tr} \left( \Theta_{*}^{\top} U D Z^{\top} Q \Sigma_{uz}^{\top} \Theta^{*} \right) \right] \right| \\ \le \frac{1}{N} \mathbb{E} \left[ \|\Sigma_{uz}\| \|U\| \|Z\| \|Q\| \|D\| \|\Theta^{*} \Theta_{*}^{\top}\|_{F} \right] \\ = \frac{1}{N} \|\Sigma_{uz}\| \|\Theta^{*} \Theta_{*}^{\top}\|_{F} \mathbb{E} \left[ \|U\| \|Z\| \|Q\| \|D\| \right] $$

We begin by observing that $\|Q\| \le \frac{1}{\lambda}$, which follows from the fact that $(ZZ^{\top}/N + \lambda I_n) \succeq \lambda I_n$. Next, we note that $\|\Sigma_{uz}\| = O(1)$. To see this, consider the definition:

$$ \begin{aligned} \|\Sigma_{uz}\| &= \sup_{\|\mathbf{a}\|_2=1, \|\mathbf{b}\|_2=1} \mathbf{a}^{\top} \Sigma_{uz} \mathbf{b} \\ &= \sup_{\|\mathbf{a}\|_2=1, \|\mathbf{b}\|_2=1} \mathbf{a}^{\top} \mathbb{E}[\mathbf{u}\mathbf{z}^{\top}] \mathbf{b} \\ &\le \sup_{\|\mathbf{a}\|_2=1, \|\mathbf{b}\|_2=1} \sqrt{\mathbb{E}[(\mathbf{a}^{\top}\mathbf{u})^2]} \sqrt{\mathbb{E}[(\mathbf{z}^{\top}\mathbf{b})^2]}. \end{aligned} $$

The first term, $\sqrt{\mathbb{E}[(\mathbf{a}^{\top}\mathbf{u})^2]}$, does not depend on $n$. For the second term, we have $\mathbb{E}[(\mathbf{z}^{\top}\mathbf{b})^2] = \mathbf{b}^{\top}\Sigma_z\mathbf{b} \le \|\Sigma_z\|$, which is bounded by lemma 3. Hence, $\|\Sigma_{uz}\| = O(1)$ as claimed. Finally $\|\Theta_{*}\Theta^{*\top}\|_{F}$ is fixed and deterministic.
That means, we have:

$$ \frac{1}{N} \|\Sigma_{uz}\| \|\Theta^* \Theta_*^\top\|_F \mathbb{E} [\|U\| \|Z\| \|Q\| \|D\|] = O(\frac{1}{N} \mathbb{E} [\|U\| \|Z\| \|D\|]) $$

Using twice Cauchy-Schwarz, we have

$$ \begin{aligned} \frac{1}{N} \mathbb{E} [\|U\| \|Z\| \|D\|] &\le \frac{1}{N} \mathbb{E}[\|U\|^2 \|Z\|^2]^{\frac{1}{2}} \mathbb{E}[\|D\|^2]^{\frac{1}{2}} \\ &\le \frac{1}{N} \mathbb{E}[\|U\|^4]^{\frac{1}{4}} \mathbb{E}[\|Z\|^4]^{\frac{1}{4}} \mathbb{E}[\|D\|^2]^{\frac{1}{2}} \end{aligned} $$

Moving forward, we apply Lemma 5 to obtain bounds on the operator norms of U and Z: specifically,

$$ \mathbb{E}[\|U\|^k] = O(\sqrt{N}^k), \quad \text{and} \quad \mathbb{E}[\|Z\|^k] = O((\sqrt{N} + \sqrt{n})^k). $$

In addition, using Lemma 4 and using the fact that $\frac{1}{N}\mathbf{z}_i^\top Q_{-i} \mathbf{z}_i \ge 0$ so that $1 + \tilde{\delta}_i \ge 1$, we get

$$ \mathbb{E}[\|D\|^k] = O\left(\frac{\|Q\|_F^k}{N^k}\right) = O\left(\frac{(\sqrt{N}\|Q\|)^k}{N^k}\right) = O\left(\frac{1}{\sqrt{N}^k}\right), $$

Putting things together, we get that $\frac{1}{N}\mathbb{E}[\|U\|\|Z\|\|D\|] = O(\frac{1}{\sqrt{N}})$, thus finally we have

$$ \frac{1}{N} \sum_{i=1}^{N} \mathbb{E} \left[ \mathrm{Tr} \left( \Theta_{*}^{\top} \frac{\mathbf{u}_{i} \mathbf{z}_{i}^{\top} Q_{-i} \Sigma_{uz}^{\top} \Theta^{*}}{1 + \tilde{\delta}_{i}} \right) \right] \rightarrow \frac{1}{N} \sum_{i=1}^{N} \mathbb{E} \left[ \mathrm{Tr} \left( \Theta_{*}^{\top} \frac{\mathbf{u}_{i} \mathbf{z}_{i}^{\top} Q_{-i} \Sigma_{uz}^{\top} \Theta^{*}}{1 + \tilde{\delta}_{i}} \right) \right] $$

**Step 3: Limit** As a result of step 1 and 2, we can take the expectation of independent random variables:

$$ \mathbb{E}[\mathrm{Tr}(\Theta_{*}^{\top} \frac{\mathbf{u}_{i} \mathbf{z}_{i}^{\top} Q_{-i} \Sigma_{uz}^{\top} \Theta^{*}}{1 + \tilde{\delta}_{i}})] = \mathrm{Tr}(\Theta_{*}^{\top} \frac{\Sigma_{uz} \mathbb{E}[Q_{-i}]}{1 + \tilde{\delta}_{i}} \Sigma_{uz}^{\top} \Theta^{*}), $$

Using Sherman Morrison (lemma 1) we have $\|Q - Q_{-i}\|_F \rightarrow 0$. And since $\|\Sigma_{uz}^\top \Theta^* \Sigma_{uz}^\top \|_F$ is bounded, we have

$$ \mathrm{Tr}(\Theta_{*}^{\top} \frac{\Sigma_{uz} \mathbb{E}[Q_{-i}]}{1 + \tilde{\delta}_{i}} \Sigma_{uz}^{\top} \Theta^{*}) \rightarrow \mathrm{Tr}(\Theta_{*}^{\top} \frac{\Sigma_{uz} \mathbb{E}[Q]}{1 + \tilde{\delta}} \Sigma_{uz}^{\top} \Theta^{*}), $$

The limit of $\mathbb{E}[Q]$ is a classical result in random matrix theory.

**Lemma 6 ([16]).** Let $\bar{Q} := (\frac{\Sigma_z}{1+\delta} + \lambda I_n)^{-1}$, we have

$$ \mathbb{E}[Q - \bar{Q}] \rightarrow 0 $$

with $\delta$ such as $\tilde{\delta} - \delta \rightarrow 0$ and verifying the fixed-point equation:

$$ \delta - \frac{1}{N} \mathrm{Tr} \left( \left( \frac{\Sigma_z}{1 + \delta} + \lambda I_n \right)^{-1} \Sigma_z \right) \rightarrow 0 $$

Given that $\tilde{\Sigma} \Theta^* \Theta_*^\top \tilde{\Sigma}^\top$ is bounded in Frobenius norm and using the above limit, we have finally

$$ \mathbb{E} \left[ \mathrm{Tr} \left( \left( \frac{1}{N} Z^\top U \Theta_* \right)^\top Q \tilde{\Sigma} \Theta^* \right) \right] - \mathrm{Tr} \left( \tilde{\Sigma} \Theta_* \Theta_*^\top \frac{\tilde{\Sigma}^\top \bar{Q}}{1 + \delta} \right) \rightarrow 0 $$
**Third Term** Let's focus on determining the limit of

$$
\begin{aligned}
\mathbb{E} \left[ \left\| \left( \frac{1}{N} Z^\top U \Theta_* \right)^\top \left( \frac{1}{N} Z^\top Z + \lambda I_n \right)^{-1} \mathbf{z}_0 \right\|^2 \right] &= \frac{1}{N^2} \mathbb{E} \left[ \Theta_*^\top U^\top Z Q \mathbf{z}_0 \mathbf{z}_0^\top Q Z^\top U \Theta_* \right] \\
&= \frac{1}{N^2} \mathbb{E} \left[ \Theta_*^\top U^\top Z Q \Sigma_z Q Z^\top U \Theta_* \right]
\end{aligned}
$$

Let's decouple $Q$ from $U$ and $Z$ as we did above. We begin by writing:

$$ \frac{1}{N^2} \mathbb{E} \left[ \Theta_*^\top U^\top Z Q \Sigma_z Q Z^\top U \Theta_* \right] = \frac{1}{N^2} \sum_{i=1}^N \mathbb{E} \left[ \Theta_*^\top \mathbf{u}_i \mathbf{z}_i^\top Q \Sigma_z Q Z^\top U \Theta_* \right] $$

Using Sherman-Morrison identity (equation 1), we have

$$ \frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} Q \Sigma_{z} Q Z^{\top} U \Theta_{*} \right] = \frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \delta_{i}} \Sigma_{z} Q Z^{\top} U \Theta_{*} \right] $$

We use now the fact that $\delta_i$ is concentrated around its expectation $\tilde{\delta}_i$

$$ \frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} Q Z^{\top} U \Theta_{*} \right] - \frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} Q Z^{\top} U \Theta_{*} \right] \to 0 $$

To justify that, we rewrite the difference between the two above sums as $\frac{1}{N^2} \mathbb{E}[\Theta_*^\top U^\top \Delta_1 Z Q \Sigma_z Q Z^\top U \Theta_*]$ where $\Delta_1 = \text{diag}(\frac{\delta_i - \delta_i}{1+\delta_i}, i \in \{1, \dots, N\})$ and $\mathbb{E}[\|\Delta_1\|^k] = O((\log N)^k), k \in \mathbb{N}$. [Justification max of subgaussians is $\log N$ (Lemma 2)] We now do the same thing on the other side

$$
\begin{aligned}
\frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} Q Z^{\top} U \Theta_{*} \right] &= \frac{1}{N^2} \sum_{i,j=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} u_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} Q \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] \\
&= \frac{1}{N^2} \sum_{i,j=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j}}{1 + \tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] \\
&\to \frac{1}{N^2} \sum_{i,j=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} u_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j}}{1 + \tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right]
\end{aligned}
$$

We justify the last step by writing the difference between the last two terms as $\frac{1}{N^2} \mathbb{E}[\Theta_*^\top U^\top \Delta_2 Z Q \Sigma_z Q Z^\top \Delta_1 U \Theta_*]$, where $\Delta_2 = \text{diag}(\frac{1+\delta_i}{1+\delta_i}, i \in \{1, \dots, N\})$ and $\mathbb{E}[\|\Delta_2\|^k] = O(1), k \in \mathbb{N}$.

We split now the sum into two terms $i = j$ and $i \neq j$,

$$
\frac{1}{N^2} \sum_{i,j=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j}}{1 + \tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] = \frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \mathbf{z}_{i} \mathbf{u}_{i}^{\top} \Theta_{*} \right] \\
+ \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j}}{1 + \tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right]
$$

For the term $i = j$, we use the concentration of the quadratic form random variable $\xi_i := \frac{1}{N} \mathbf{z}_i^\top Q_{-i} \Sigma_z Q_{-i} \mathbf{z}_i$ around its expectation $\tilde{\xi}_i$,

$$ \frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-i}}{1 + \tilde{\delta}_{i}} \mathbf{z}_{i} \mathbf{u}_{i}^{\top} \Theta_{*} \right] - \frac{1}{N} \sum_{i=1}^{N} \frac{\tilde{\xi}_{i}}{(1 + \tilde{\delta}_{i})^2} \Theta_{*}^{\top} \Sigma_{u} \Theta_{*} \to 0 $$
The justification is the same as above, we write the difference as $\frac{1}{N} \sum_{i=1}^{N} \mathbb{E}[\Theta_{*}^{\top} U^{\top} \Delta_{3} U \Theta_{*}]$, where $\Delta_{3} = \text{diag}(\frac{\xi_{i}-\tilde{\xi}_{i}}{(1+\tilde{\delta}_{i})^{2}}, i \in \{1, \dots, N\})$ and $\mathbb{E}[\|\Delta_{3}\|^{k}] = O((\frac{\log N}{\sqrt{N}})^{k})$, $k \in \mathbb{N}$.

For the term $i \neq j$, we still have dependency between $Q_{-i}$ and $Q_{-j}$, so we reapply Sherman-Morrison (equation 1 this time):

$$ \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1+\tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j}}{1+\tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] = \kappa_{1} + \kappa_{2} + \kappa_{3} + \kappa_{4} $$

with

$$ \kappa_{1} = \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i,j}}{1+\tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-i,j}}{1+\tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] $$

$$ \kappa_{2} = -\frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{1}{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i} \mathbf{z}_{j} \mathbf{z}_{j}^{\top} Q_{-i,j}}{1+\tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j}}{1+\tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] $$

$$ \kappa_{3} = -\frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{1}{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i}}{1+\tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j} \mathbf{z}_{i} \mathbf{z}_{i}^{\top} Q_{-i,j}}{1+\tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] $$

$$ \kappa_{4} = \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{1}{N^2} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i} \mathbf{z}_{j} \mathbf{z}_{j}^{\top} Q_{-i,j}}{1+\tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j} \mathbf{z}_{i} \mathbf{z}_{i}^{\top} Q_{-i,j}}{1+\tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] $$

It is now straightforward that

$$ \kappa_{1} = \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{\Theta_{*}^{\top} \Sigma_{uz} \mathbb{E}[Q_{-i,j} \Sigma_{z} Q_{-i,j}] \Sigma_{uz}^{\top} \Theta_{*}}{(1+\tilde{\delta}_{i})(1+\tilde{\delta}_{j})} $$

For $\kappa_2$, we will redo Sherman-Morisson (equation 1)

$$ \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{1}{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i} \mathbf{z}_{j} \mathbf{z}_{j}^{\top} Q_{-i,j}}{1+\tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j}}{1+\tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] \\ = \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{1}{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i} \mathbf{z}_{j} \mathbf{z}_{j}^{\top} Q_{-i,j}}{1+\tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-i,j}}{1+\tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] \\ - \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{1}{N^2} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i,j} \mathbf{z}_{j} \mathbf{z}_{j}^{\top} Q_{-i,j}}{(1+\tilde{\delta}_{i})(1+\frac{1}{N} \mathbf{z}_{j}^{\top} Q_{-i,j} \mathbf{z}_{j})} \Sigma_{z} \frac{Q_{-i,j} \mathbf{z}_{i} \mathbf{z}_{i}^{\top} Q_{-i,j}}{(1+\tilde{\delta}_{j})(1+\frac{1}{N} \mathbf{z}_{i}^{\top} Q_{-i,j} \mathbf{z}_{i})} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] $$

The second term is negligible. In fact using Hanson Wright and the fact that $\|\boldsymbol{\mu}_z\|_2 = O(1)$ ($\mathbb{E}[\mathbf{z}_i A \mathbf{z}_j] = \boldsymbol{\mu}_z A \boldsymbol{\mu}_z$, $i \neq j$ with $\|A\| = O(1)$), we have $\mathbb{E}[\|\mathbf{z}_i A \mathbf{z}_j\|^k] = O(\sqrt{N}^k)$. We have 3 random variables of that type, so using Cauchy-Schwarz, we get that the second term is $O(\frac{\sqrt{N}^3}{N^2}) = o(1)$. We thus have

$$ \kappa_2 - \left( -\frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{1}{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} \frac{Q_{-i} \mathbf{z}_{j} \mathbf{z}_{j}^{\top} Q_{-i,j}}{1+\tilde{\delta}_{i}} \Sigma_{z} \frac{Q_{-j}}{1+\tilde{\delta}_{j}} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] \right) \to 0 $$
We then use the concentration of the random variable $\eta_{i,j} = \frac{1}{N}z_j Q_{-i,-j} \Sigma_z Q_{-i,j} z_j$ around its expectation $\tilde{\eta}_{i,j}$. To justify it, we take the difference

$$ 
\begin{aligned}
& \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{1}{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} Q_{-i,j} z_{j} \frac{\eta_{i,j} - \tilde{\eta}_{i,j}}{(1 + \tilde{\delta}_{i})(1 + \tilde{\delta}_{j})(1 + \tilde{\delta}_{j})} \mathbf{u}_{j}^{\top} \Theta_{*} \right] \\
&= \frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} Q_{-i} \sum_{\substack{j=1 \\ i \neq j}}^{N} \mathbf{z}_{j} \frac{\eta_{i,j} - \tilde{\eta}_{i,j}}{(1 + \tilde{\delta}_{i})(1 + \tilde{\delta}_{j})} u_{j}^{\top} \Theta_{*} \right] \\
&= \frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \Sigma_{uz} Q_{-i} \sum_{\substack{j=1 \\ i \neq j}}^{N} \mathbf{z}_{j} \frac{\eta_{i,j} - \tilde{\eta}_{i,j}}{(1 + \tilde{\delta}_{i})(1 + \tilde{\delta}_{j})} \mathbf{u}_{j}^{\top} \Theta_{*} \right] \\
&= \frac{1}{N^2} \sum_{i=1}^{N} \mathbb{E} \left[ \Theta_{*}^{\top} \Sigma_{uz} Q_{-i} Z^{\top} \Delta_{i} U \Theta_{*} \right]
\end{aligned}
$$

where $\Delta_i = \text{diag}(\frac{\eta_{i,j} - \tilde{\eta}_{i,j}}{(1+\tilde{\delta}_i)(1+\tilde{\delta}_j)}, j \in \{1, \dots, N\} \setminus \{i\}$ and $0$ for $i = j)$ and $\mathbb{E}[\|\Delta_i\|^k] = O((\frac{\log N}{\sqrt{N}})^k), k \in \mathbb{N}$.
This means that

$$ \kappa_2 - \left( -\frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{\tilde{\eta}_{i,j}}{(1 + \tilde{\delta}_i)(1 + \tilde{\delta}_j)} \mathbb{E} \left[ \Theta_{*}^{\top} \mathbf{u}_{i} \mathbf{z}_{i}^{\top} Q_{-i} \mathbf{z}_{j} \mathbf{u}_{j}^{\top} \Theta_{*} \right] \right) \to 0 $$

As we did above, this gives as that

$$ \kappa_2 - \left( -\frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \frac{\tilde{\eta}_{i,j}}{(1 + \tilde{\delta}_i)(1 + \tilde{\delta}_j)(1 + \tilde{\delta}_{i,j})} \Theta_{*}^{\top} \Sigma_{uz} \mathbb{E}[Q_{-i,j}] \Sigma_{uz}^{\top} \Theta_{*} \right) \to 0 $$

$\kappa_3$ is the same technique and $\kappa_4$ is negligible as we did above.

The limit of $\mathbb{E}[Q\Sigma_z Q]$ is also a classical result in random matrix theory :

**Lemma 7 ([5]).** Let $Q := \frac{(1+\delta)^2}{(1+\delta)^2 - \frac{1}{N} \text{Tr}(\Sigma_z Q \Sigma_z Q)} \bar{Q} \Sigma_z \bar{Q}$, we have

$$ \mathbb{E}[Q\Sigma_z Q - Q] \to 0. $$

Now as we did above in the step 3 limit A.2.1, we can replace everything with its limit, then we get the desired result after simplifications.

## A.3 Variance Term

In the same manner, we want the limit of

$$ \mathbb{E} \left[ \left( z_0^\top \left( \frac{Z^\top Z}{N} + \lambda I_n \right)^{-1} \frac{Z^\top Z}{N^2} \left( \frac{Z^\top Z}{N} + \lambda I_n \right)^{-1} z_0 \right) \right] = \mathbb{E} \left[ \text{Tr} \left( Q \frac{Z^\top Z}{N^2} Q \Sigma_z \right) \right], $$
We write:

$$
\begin{align*}
\mathbb{E} \left[ \mathrm{Tr} \left( Q \frac{Z^\top Z}{N^2} Q \Sigma_z \right) \right] &= \frac{1}{N^2} \sum_{i=1}^N \mathbb{E} \left[ \mathrm{Tr} \left( Q \mathbf{z}_i \mathbf{z}_i^\top Q \Sigma_z \right) \right] \\
&= \frac{1}{N^2} \sum_{i=1}^N \mathbb{E} \left[ \mathrm{Tr} \left( \frac{Q_{-i}}{1 + \frac{1}{N} \mathbf{z}_i^\top Q_{-i} \mathbf{z}_i} \mathbf{z}_i \mathbf{z}_i^\top \frac{Q_{-i}}{1 + \frac{1}{N} \mathbf{z}_i^\top Q_{-i} \mathbf{z}_i} \Sigma_z \right) \right]
\end{align*}
$$

As we did for the second term of the bias, we can here replace $\frac{1}{N} \mathbf{z}_i^\top Q_{-i} \mathbf{z}_i$ with its expectation in
the limit. Then using Lemma 7, we get the desired result after simplification.

B Proof of theorem 2

Let $\rho' := \varphi \cdot \rho(W_0)$, we have

$$
S = \left[ \left( \frac{W_0}{\rho'} \right)^{T-1} \mathbf{w}_{\text{in}}, \left( \frac{W_0}{\rho'} \right)^{T-2} \mathbf{w}_{\text{in}}, \dots, \left( \frac{W_0}{\rho'} \right)^0 \mathbf{w}_{\text{in}} \right] \in \mathbb{R}^{n \times T}
$$

denote the state matrix built from the input weight vector $\mathbf{w}_{\text{in}} \in \mathbb{R}^n$ and recurrent matrix $W \in \mathbb{R}^{n \times n}$.
We get $Z = SU$.

**Lemma 8 (Concentration of $S^\top S$).** Let $W_0 \in \mathbb{R}^{n \times n}$ have i.i.d. standard Gaussian entries and let $\mathbf{w}_{\text{in}} \sim \mathcal{N}(0, \frac{1}{n} I_n)$ and $\rho' = \varphi \cdot \rho(W_0)$ with fixed $\varphi < 1$, and let $S \in \mathbb{R}^{T \times n}$ be defined as in 7. Then, for some constant $C > 0$ depending only on $T$ and $\varphi$, we have

$$
\mathbb{P} \left( \max_{i,j \le T} |[S^\top S]_{ij} - \mathbb{E}[S^\top S]_{ij}| > \frac{C}{\sqrt{n}} \right) \le c_1 e^{-c_2 n}
$$

for some constants $c_1, c_2 > 0$. In particular,

$$
\|S^\top S - \mathbb{E}[S^\top S]\|_F = O_P\left(\frac{1}{\sqrt{n}}\right).
$$

*Proof.* For each $i, j \in \{0, \dots, T-1\}$, define

$$
A_{ij} := \left( \frac{W_0^\top}{\rho'} \right)^i \left( \frac{W_0}{\rho'} \right)^j \in \mathbb{R}^{n \times n},
$$

so that

$$
[S^\top S]_{ij} = \langle S_i, S_j \rangle = \mathbf{w}_{\text{in}}^\top A_{ij} \mathbf{w}_{\text{in}}.
$$

From [2] we have almost surely as $n$ tends to infinity $\rho(W_0)/\sqrt{n} \to 1$, and from Bai-Yin Theorem ([25] + Remark on i.i.d. standard random matrices), we have almost surely, $\|W_0\|/\sqrt{n} \to 2$. Fix $\varepsilon \in (0,1)$. Then, almost surely for all sufficiently large $n$,

$$
\rho(W_0) \ge (1 - \varepsilon)\sqrt{n}, \quad \|W_0\| \le (2 + \varepsilon)\sqrt{n}.
$$

Since $\rho' = \varphi\rho(W_0)$ with fixed $\varphi \in (0,1)$, it follows that

$$
\left\| \frac{W_0}{\rho'} \right\| = \frac{\|W_0\|}{\varphi \rho(W_0)} \le \frac{2 + \varepsilon}{\varphi(1 - \varepsilon)} =: C_{\varepsilon},
$$
for all large $n$. In particular, for each $i, j \le T - 1$,

$$ \|A_{ij}\| = \left\| \left(\frac{W_0^\top}{\rho'}\right)^i \left(\frac{W_0}{\rho'}\right)^j \right\| \le \|W_0/\rho'\|^{i+j} \le C_\epsilon^{2(T-1)} =: M, $$

and hence $\|A_{ij}\|_F \le \sqrt{n} M$.

Now write $X := \sqrt{n} \mathbf{w}_{\text{in}} \sim \mathcal{N}(0, I_n)$. Then

$$ [S^\top S]_{ij} - \mathbb{E}[S^\top S]_{ij} = \frac{1}{n} (X^\top A_{ij} X - \mathbb{E}[X^\top A_{ij} X]). $$

By the Hanson-Wright inequality ($W_0$ independent of $X$),

$$ \mathbb{P}(|X^\top A_{ij} X - \mathbb{E}[X^\top A_{ij} X| > u) \le 2 \exp\left(-c \min\left\{\frac{u^2}{\|A_{ij}\|_F^2}, \frac{u}{\|A_{ij}\|}\right\}\right). $$

Taking $u = nt$ and using $\|A_{ij}\| \le M$, $\|A_{ij}\|_F \le \sqrt{n} M$ gives

$$ \mathbb{P}(|[S^\top S]_{ij} - \mathbb{E}[S^\top S]_{ij}| > t) \le 2 \exp\left(-c \min\left\{\frac{nt^2}{M^2}, \frac{nt}{M}\right\}\right). $$

In particular, for $t = x/\sqrt{n}$,

$$ \mathbb{P}(|[S^\top S]_{ij} - \mathbb{E}[S^\top S]_{ij}| > \frac{x}{\sqrt{n}}) \le 2 \exp\left(-c \frac{x^2}{M^2}\right). $$

A union bound over all $i, j \le T - 1$ yields

$$ \mathbb{P}\left(\max_{i,j} \left| [S^\top S]_{ij} - \mathbb{E}[S^\top S]_{ij} \right| > \frac{x}{\sqrt{n}} \right) \le 2T^2 \exp\left(-c \frac{x^2}{M^2}\right), $$

□

**Lemma 9** (Limit of $\mathbb{E}[S^\top S]$). As $T$ is fixed, we have as $n \to \infty$,

$$ \mathbb{E}[S^\top S] \to \text{diag} \left( \varphi^{-(T-1)}, \varphi^{-(T-2)}, \dots, \varphi^{-1}, 1 \right) $$

*Proof.* We have

$$ [S^\top S]_{i,j} = \left\langle \left(\frac{W_0}{\rho'}\right)^{T-i} \mathbf{w}_{\text{in}}, \left(\frac{W_0}{\rho'}\right)^{T-j} \mathbf{w}_{\text{in}} \right\rangle, \quad 1 \le i, j \le T $$

Because $\mathbf{w}_{\text{in}}$ is independent of $W_0$, we get

$$ \mathbb{E}_{\mathbf{w}_{\text{in}}} [S^\top S]_{ij} = \frac{1}{n} (\rho')^{-2T+i+j} \text{Tr}(W_0^{\top (T-i)} W_0^{T-j}). $$

Since $W_0$ has i.i.d. $\mathcal{N}(0, 1)$, we have, using [12] [Prop + Remark 1.4, Prop 1.7], that

$$ \frac{1}{n} \text{Tr} \left( \left( \frac{W_0}{\sqrt{n}} \right)^{\top (T-i)} \left( \frac{W_0}{\sqrt{n}} \right)^{T-j} \right) \xrightarrow{n \to \infty} \delta_{ij} \quad \text{a.s.} $$

From [2], we have $\rho(W_0)/\sqrt{n} \to 1$ almost surely.

We insert this back to get:

$$ \mathbb{E}_{\mathbf{w}_{\text{in}}} [S^\top S]_{ij} \xrightarrow{n \to \infty} \delta_{ij} \varphi^{-(T-i)} \quad \text{a.s.} $$
Therefore:

$$ \mathbb{E}_{\mathbf{w}_{\text{in}}}[S^{\top}S] \xrightarrow{n \to \infty} \operatorname{diag}(\varphi^{-(T-1)}, \varphi^{-(T-2)}, \dots, \varphi^{-1}, 1) \quad \text{a.s.} $$

For every $n$, we have $\mathbb{E}_{\mathbf{w}_{\text{in}}}[S^\top S]_{ij}$ is almost surely bounded. In particular it is uniformly integral, that is for any $K > 1$,

$$ \sup_n \mathbb{E} \left[ \left| \mathbb{E}_{\mathbf{w}_{\text{in}}}[S^\top S]_{ij} \right| \cdot \mathbf{1}_{\{|\mathbb{E}_{\mathbf{w}_{\text{in}}}[S^\top S]_{ij}| > K\}} \right] = 0, $$

Thus we have the convergence in expectation, that is

$$ \mathbb{E}[S^\top S] \rightarrow \operatorname{diag}(\varphi^{-(T-1)}, \varphi^{-(T-2)}, \dots, \varphi^{-1}, 1) $$

□

Now that we have established the limit of $S^\top S$, we can proceed to compute the limit of the risk. As mentioned in Remark 1, under the hypothesis that $\mathbf{u}$ is concentrated, the vector $\mathbf{z}$, obtained as a Lipschitz function of $\mathbf{u}$ with a bounded Lipschitz constant, is also concentrated. We can thus apply Theorem 1, let $\tilde{S} := S\Sigma_u^{\frac{1}{2}}$ and $\tilde{\Theta}_* = \Sigma_u^{\frac{1}{2}}\Theta_*$ We have the asymptotic bias is given by:

$$ \mathcal{B}^2 \rightarrow \mathcal{B}_{\infty}^2 := \frac{1}{1-\alpha} \left( \operatorname{Tr}(\tilde{\Theta}^{*\top} \tilde{\Theta}^*) - \frac{2}{1+\delta} \operatorname{Tr}(\tilde{\Theta}^{*\top} [\tilde{S}^\top \bar{Q} \tilde{S}] \tilde{\Theta}^*) + \frac{1}{(1+\delta)^2} \operatorname{Tr}(\tilde{\Theta}^{*\top} [\tilde{S}^\top \bar{Q} \tilde{S}]^2 \tilde{\Theta}^*) \right) $$

where

$$ \bar{Q} := \left( \frac{\tilde{S}\tilde{S}^{\top}}{1+\delta} + \lambda I_n \right)^{-1}, \quad \delta = \frac{1}{N} \operatorname{Tr}[\tilde{S}^{\top}\bar{Q}\tilde{S}], \quad \alpha = \frac{\operatorname{Tr}([\tilde{S}^{\top}\bar{Q}\tilde{S}]^2)}{N(1+\delta)^2} $$

Since we know the expectation of $\tilde{S}^\top \tilde{S} \in \mathbb{R}^{T \times T}$, we will use the Woodbury identity to get it in the expression instead of $\tilde{S}\tilde{S}^\top$:

$$ \tilde{S}^{\top} \bar{Q} \tilde{S} = \tilde{S}^{\top} \tilde{S} \left( \frac{\tilde{S}^{\top} \tilde{S}}{1+\delta} + \lambda I_T \right)^{-1} $$

That makes $\mathcal{B}_{\infty}^2$ a function of $\tilde{S}^\top \tilde{S}$. Our goal is to use the fact that the map $\tilde{S}^\top \tilde{S} \mapsto \mathcal{B}_{\infty}^2(\tilde{S}^\top \tilde{S})$ is Lipschitz in Frobenius norm with Lipschitz constant $L = O(1)$, so that we may formally replace $\tilde{S}^\top \tilde{S}$ by the limit of its expectation inside the global expectation. That is, we want (by defining $M_{\infty} := \Sigma_u^{\frac{1}{2}} \operatorname{diag}(\varphi^{(i-T)}) \Sigma_u^{\frac{1}{2}}$):

$$ \begin{aligned} |\mathcal{B}_{\infty}^2(\tilde{S}^{\top} \tilde{S}) - \mathcal{B}_{\infty}^2(M_{\infty})| & \le |\mathcal{B}_{\infty}^2(\tilde{S}^{\top} \tilde{S}) - \mathcal{B}_{\infty}^2(\mathbb{E}[\tilde{S}^{\top} \tilde{S}])| + |\mathcal{B}_{\infty}^2(\mathbb{E}[\tilde{S}^{\top} \tilde{S}]) - \mathcal{B}_{\infty}^2(M_{\infty})| \\ & \le L \|\tilde{S}^{\top} \tilde{S} - \mathbb{E}[\tilde{S}^{\top} \tilde{S}]\|_F + L \|\mathbb{E}[\tilde{S}^{\top} \tilde{S}] - M_{\infty}\|_F \end{aligned} $$

So that

$$ \begin{aligned} |\mathbb{E}[\mathcal{B}_{\infty}^2(\tilde{S}^{\top} \tilde{S})] - \mathcal{B}_{\infty}^2(M)| & \le L \mathbb{E} \left[ \left\| \tilde{S}^{\top} \tilde{S} - \mathbb{E}[\tilde{S}^{\top} \tilde{S}] \right\|_F \right] + L \left\| \mathbb{E}[\tilde{S}^{\top} \tilde{S}] - M \right\|_F \\ & \rightarrow 0. \end{aligned} $$

Let us define $M := \tilde{S}^\top \tilde{S} \in \mathbb{R}^{T \times T}$ and let's prove the Lipschitzness of the map $M \mapsto \mathcal{B}_{\infty}^2(M)$. We have

$$ \mathcal{B}_{\infty}^2(M) = \frac{1}{1-\alpha(M)} \left( t_1 - \frac{2}{1+\delta(M)} t_2(M) + \frac{1}{(1+\delta(M))^2} t_3(M) \right) $$
where

$$
\begin{align*}
t_1 &:= \operatorname{Tr}(\tilde{\Theta}^* \tilde{\Theta}^*), & t_2(M) &:= \operatorname{Tr}(\tilde{\Theta}^* A(M) \tilde{\Theta}^*), \\
t_3(M) &:= \operatorname{Tr}(\tilde{\Theta}^* A(M)^2 \tilde{\Theta}^*), & A(M) &:= M \left( \frac{M}{1+\delta} + \lambda I_T \right)^{-1}, \\
\delta(M) &:= \frac{1}{N} \operatorname{Tr}(M \left( \frac{M}{1+\delta(M)} + \lambda I_T \right)^{-1}), & \alpha(M) &:= \frac{\operatorname{Tr}(A(M)^2)}{N(1+\delta(M))^2}.
\end{align*}
$$

**Definition 6** (Lipschitz constant and Uniform bound). *The Lipschitz constant of a function f is defined as*

$$
L_f := \inf \{ L > 0 : |f(x) - f(y)| \le L \|x - y\| \quad \forall x, y \}.
$$

The uniform bound of f is defined as

$$
B_f := \sup_x |f(x)|.
$$

Let us first start by showing the Lipschitzness of the map $M \mapsto \delta(M)$.

**Lemma 10** (Uniform bound and Lipschitzness of $\delta$). Let $\lambda > 0$ and integers $N > T$. For $M \succeq 0 \in \mathbb{R}^{T \times T}$ and $\delta \ge 0$ define

$$
\varphi(\delta, M) := \frac{1}{N} \operatorname{Tr} \left( \left( \frac{M}{1+\delta} + \lambda I_T \right)^{-1} M \right) - \delta.
$$

For each $M \ge 0$, let $\delta(M) \ge 0$ be any solution of $\varphi(\delta(M), M) = 0$ (e.g. [15]). Then

$$
0 \le \delta(M) \le \frac{T}{N-T} \quad \text{and} \quad \left|\delta(M_1) - \delta(M_2)\right| \le \frac{\sqrt{T}}{\lambda(N-T)} \|M_1 - M_2\|_F,
$$

i.e. $B_{\delta} \le T/(N - T)$ and the map $M \mapsto \delta(M)$ is Lipschitz (w.r.t. $\| \cdot \|_F$) with constant $L_{\delta} \le \sqrt{T}/(\lambda(N - T))$.

*Proof.* Fix $M \ge 0$ and set

$$
Q := \left( \frac{M}{1+\delta} + \lambda I_T \right)^{-1}, \quad P := Q - \frac{1}{1+\delta} Q M Q.
$$

Note that $Q \ge 0$ and $Q^{-1} \ge \lambda I_T$, hence $\|Q\| \le \lambda^{-1}$.

**Uniform bound.** At a fixed point $\varphi(\delta(M), M) = 0$, letting $\{\sigma_i\}_{i=1}^r$ be the nonzero eigenvalues of $M$ with $r = \operatorname{rank}(M) \le T$,

$$
\delta = \frac{1}{N} \sum_{i=1}^{r} \frac{\sigma_i}{\sigma_i/(1+\delta) + \lambda} = \frac{1}{N} \sum_{i=1}^{r} \frac{(1+\delta)\sigma_i}{\sigma_i + \lambda(1+\delta)} \le \frac{r}{N}(1+\delta) \le \frac{T}{N}(1+\delta).
$$

Rearranging gives $\delta \le T/(N - T)$.
**Lipschitzness.** Differentiate $\varphi$:

$$ \partial_{\delta}\varphi = \frac{1}{N(1+\delta)^2} \mathrm{Tr}(QMQM) - 1, \quad d_M\varphi[H] = \frac{1}{N} \mathrm{Tr}(PH) \quad (\forall H = H^\top). $$

From the spectral decomposition of $M$,

$$ \mathrm{Tr}(QMQM) = \sum_{i:\sigma_i>0} \frac{\sigma_i^2}{(\sigma_i/(1+\delta) + \lambda)^2} \le r(1+\delta)^2 \le T(1+\delta)^2, $$

so at a fixed point,

$$ -\partial_{\delta}\varphi(\delta(M), M) \ge 1 - \frac{T}{N} = \frac{N-T}{N} > 0. $$

Using $Q^{-1} = \frac{M}{1+\delta} + \lambda I_T$,

$$ \frac{1}{1+\delta} QMQ = Q(Q^{-1} - \lambda I_T)Q = Q - \lambda Q^2, $$

hence $P = \lambda Q^2$. Therefore

$$ \|P\| \le \lambda \|Q\|^2 \le \lambda \cdot \lambda^{-2} = \lambda^{-1}, \quad \|P\|_F \le \sqrt{T} \|P\| \le \frac{\sqrt{T}}{\lambda}, $$

and thus

$$ |d_M\varphi[H]| = \frac{1}{N} |\mathrm{Tr}(PH)| \le \frac{1}{N} \|P\|_F \|H\|_F \le \frac{\sqrt{T}}{\lambda N} \|H\|_F. $$

By the implicit function theorem,

$$ D\delta(M)[H] = -\frac{d_M\varphi(\delta(M), M)[H]}{\partial_{\delta}\varphi(\delta(M), M)}, $$

whence

$$ |D\delta(M)[H]| \le \frac{\frac{\sqrt{T}}{\lambda N}}{\frac{N-T}{N}} \|H\|_F = \frac{\sqrt{T}}{\lambda(N-T)} \|H\|_F. $$

Taking a supremum over $\|H\|_F = 1$ yields $L_\delta \le \sqrt{T}/(\lambda(N-T))$. $\square$

**Lemma 11** (Uniform bound and Lipschitzness of $A$). Let

$$ A(M) := M \left( \lambda I_T + \frac{M}{1+\delta(M)} \right)^{-1}, \quad M \ge 0, $$

where $\delta(M)$ is defined as in Lemma 10. Then

$$ B_A \le 1 + B_\delta, \quad \text{and} \quad L_A \le \frac{1}{\lambda} + L_\delta. $$

*Proof.* Throughout, $\leq$ the Loewner order. For brevity write, for $s \in (0, 1]$,

$$ Q_s(M) := (\lambda I_T + sM)^{-1}, \quad Q(M) := Q_{s(M)}(M). $$

Note that $M \ge 0$ implies $0 < \lambda I_T \le \lambda I_T + sM$, hence $\|Q_s(M)\| \le \lambda^{-1}$ and $0 \le Q_s(M) \le \lambda^{-1} I_T$. Using $X(\lambda I_T + sX)^{-1} = \frac{1}{s}(I_T - \lambda(\lambda I_T + sX)^{-1})$ for any $X \ge 0$, we have

$$ A(M) = MQ(M) = \frac{1}{s(M)}(I_T - \lambda Q(M)). $$
**Uniform bound.** From $0 \le Q(M) \le \lambda^{-1}I_T$ we get $0 \le I_T - \lambda Q(M) \le I_T$. We have

$$0 \le A(M) \le \frac{1}{s(M)} I_T = (1 + \delta(M)) I_T,$$

hence $\|A(M)\| \le 1 + \delta(M)$ and therefore

$$B_A \le 1 + B_\delta.$$

**Lipschitzness.** Let $M_1, M_2 \ge 0$, and set $s_i := s(M_i)$ and $\delta_i := \delta(M_i)$. Decompose

$$A(M_1) - A(M_2) = \underbrace{M_1(\lambda I_T + s_1 M_1)^{-1} - M_2(\lambda I_T + s_1 M_2)^{-1}}_{(I)} + \underbrace{M_2(\lambda I_T + s_1 M_2)^{-1} - M_2(\lambda I_T + s_2 M_2)^{-1}}_{(II)}.$$

(I) Fixed $s$ part is $1/\lambda$-Lipschitz. For any $X, Y \ge 0$ and fixed $s \in (0, 1]$,

$$X(\lambda I_T + sX)^{-1} - Y(\lambda I_T + sY)^{-1} = \frac{\lambda}{s} [(\lambda I_T + sY)^{-1} - (\lambda I_T + sX)^{-1}].$$

By the resolvent identity,

$$(\lambda I_T + sY)^{-1} - (\lambda I_T + sX)^{-1} = (\lambda I_T + sY)^{-1} s(Y - X) (\lambda I_T + sX)^{-1}.$$

Taking norms and using $\|(\lambda I_T + sZ)^{-1}\| \le \lambda^{-1}$ for $Z \ge 0$ gives

$$\|X(\lambda I_T + sX)^{-1} - Y(\lambda I_T + sY)^{-1}\| \le \frac{\lambda}{s} \cdot \frac{1}{\lambda} \cdot s \cdot \frac{1}{\lambda} \|X - Y\| = \frac{1}{\lambda} \|X - Y\|.$$

Thus $\|(I)\| \le \frac{1}{\lambda}\|M_1 - M_2\|$.

(II) Varying $s$. Using again the resolvent identity in the $s$-parameter,

$$Q_{s_1}(M) - Q_{s_2}(M) = (s_2 - s_1) Q_{s_1}(M) M Q_{s_2}(M),$$

and multiplying by $M$ on the left,

$$MQ_{s_1}(M) - MQ_{s_2}(M) = (s_2 - s_1) MQ_{s_1}(M) MQ_{s_2}(M).$$

We have $\|MQ_s(M)\| \le \frac{1}{s}$, therefore,

$$\|MQ_{s_1}(M) - MQ_{s_2}(M)\| \le |s_2 - s_1| \|MQ_{s_1}(M)\| \|MQ_{s_2}(M)\| \le \frac{|s_2 - s_1|}{s_1 s_2}.$$

With $s_i = (1 + \delta_i)^{-1}$ we have the exact identity

$$\frac{|s_2 - s_1|}{s_1 s_2} = \left| \frac{1}{s_1} - \frac{1}{s_2} \right| = |(1 + \delta_1) - (1 + \delta_2)| = |\delta(M_1) - \delta(M_2)|.$$

Hence

$$\|(II)\| \le |\delta(M_1) - \delta(M_2)| \le L_\delta \|M_1 - M_2\|.$$

Combining (I) and (II) yields

$$\|A(M_1) - A(M_2)\| \le \left(\frac{1}{\lambda} + L_\delta\right) \|M_1 - M_2\|,$$

so $L_A \le \frac{1}{\lambda} + L_\delta$, as claimed.

□
**Lemma 12** (Bounds and Lipschitzness of $t_2$ and $t_3$). We have

$$B_{t_2} \leq B_A t_1, \quad B_{t_3} \leq B_A^2 t_1.$$

Moreover, $t_2$ and $t_3$ are Lipschitz in Frobenius norm with

$$L_{t_2} \leq t_1 L_A \quad L_{t_3} \leq 2t_1 B_A L_A$$

*Proof.* Set $G := \tilde{\Theta}^{*\top}\tilde{\Theta}^* \geq 0$, so that $\mathrm{Tr}(G) = t_1$.

**Uniform bounds.** Since $A(M) \geq 0$,

$$t_2(M) = \mathrm{Tr}(GA(M)) = \mathrm{Tr}(G^{1/2}A(M)G^{1/2}) \leq \|A(M)\| \mathrm{Tr}(G) \leq B_A t_1.$$

Similarly,

$$t_3(M) = \mathrm{Tr}(GA(M)^2) \leq \|A(M)^2\| \mathrm{Tr}(G) = \|A(M)\|^2 t_1 \leq B_A^2 t_1.$$

**Lipschitzness.** For $M_1, M_2$,

$$|t_2(M_1) - t_2(M_2)| = |\mathrm{Tr}(G(A(M_1) - A(M_2)))| \leq t_1 \|A(M_1) - A(M_2)\|.$$

Using the Lipschitz property of A gives

$$|t_2(M_1) - t_2(M_2)| \leq t_1 L_A \|M_1 - M_2\|_F.$$

Note

$$A(M_1)^2 - A(M_2)^2 = (A(M_1) - A(M_2))A(M_1) + A(M_2)(A(M_1) - A(M_2)).$$

Hence

$$\|A(M_1)^2 - A(M_2)^2\| \leq (\|A(M_1)\| + \|A(M_2)\|) \|A(M_1) - A(M_2)\| \leq 2B_A \|A(M_1) - A(M_2)\|.$$

Therefore

$$|t_3(M_1) - t_3(M_2)| = |\mathrm{Tr}(G(A(M_1)^2 - A(M_2)^2))| \leq t_1 \|A(M_1)^2 - A(M_2)^2\| \leq 2t_1 B_A L_A \|M_1 - M_2\|_F.$$

□

**Lemma 13** (Product of bounded Lipschitz functions). Let $f, g: D \to \mathbb{R}$ satisfy

$$|f(M)| \leq B_f, \quad |g(M)| \leq B_g \quad (\forall M \in D),$$

Then the product $h(M) := f(M)g(M)$ is Lipschitz with constant

$$L_h \leq B_f L_g + B_g L_f.$$

**Lemma 14** (Uniform bound and Lipschitzness of $M \mapsto (1 - \alpha(M))^{-1}$). Let $N > T$, $\lambda > 0$, and $M \geq 0 \in \mathbb{R}^{T \times T}$. Define $\alpha(M) := \frac{\mathrm{Tr}(A(M)^2)}{N(1+\delta(M))^2}$. Then the map $M \mapsto \frac{1}{1-\alpha(M)}$ is uniformly bounded and Lipschitz (with respect to $\| \cdot \|_F$), with

$$B_{(1-\alpha)^{-1}} \leq \frac{N}{N-T}, \quad L_{(1-\alpha)^{-1}} \leq \left(\frac{N}{N-T}\right)^2 \frac{1}{N} \left[2\sqrt{T}(1+B_\delta)L_A + 2T(1+B_\delta)^2 L_\delta\right].$$
Proof. Uniform bound. From Lemma 11, $\|A(M)\| \le 1 + \delta(M)$. Hence

$$ \alpha(M) = \frac{\mathrm{Tr}(A(M)^2)}{N(1 + \delta(M))^2} \le \frac{T \|A(M)\|^2}{N(1 + \delta(M))^2} \le \frac{T}{N}. $$

Since $N > T$, this gives $0 \le \alpha(M) \le T/N < 1$, and therefore

$$ \frac{1}{1 - \alpha(M)} \le \frac{1}{1 - T/N} = \frac{N}{N - T}. $$

Lipschitzness. Write $\alpha(M) = \frac{1}{N}g(M)h(M)$, where

$$ g(M) := \mathrm{Tr}(A(M)^2), \quad h(M) := (1 + \delta(M))^{-2}. $$

We have

$$ L_g \le 2\sqrt{T}(1 + B_\delta)L_A, \quad |g(M)| \le T(1 + B_\delta)^2, \quad L_h \le 2L_\delta, \quad |h(M)| \le 1. $$

By Lemma 13,

$$ L_{\alpha} \le \frac{1}{N}(B_g L_h + B_h L_g) = \frac{1}{N}[2\sqrt{T}(1 + B_{\delta})L_A + 2T(1 + B_{\delta})^2 L_{\delta}]. $$

Now, the function $x \mapsto (1-x)^{-1}$ has derivative $(1-x)^{-2}$, so over $[0, T/N]$ the Lipschitz factor is at most $(N/(N-T))^2$. Hence

$$ L_{(1-\alpha)^{-1}} \le \left( \frac{N}{N-T} \right)^2 L_{\alpha}. $$

This proves the claim.

<figure><img src="image_5.png" /></figure>

Now, using Lemma 13 and the fact that each function is Lipschitz and bounded, we obtain that $M \mapsto \mathcal{B}_{\infty}^2(M)$ is Lipschitz, with constant $L$ of order $O(1)$. That is, as we said above (by defining $M_{\infty} := \Sigma_u^{\frac{1}{2}} \mathrm{diag}(\varphi^{(i-T)}) \Sigma_u^{\frac{1}{2}}$), we have

$$ \begin{aligned} |\mathcal{B}_{\infty}^2(\tilde{S}^{\top}\tilde{S}) - \mathcal{B}_{\infty}^2(M_{\infty})| & \le |\mathcal{B}_{\infty}^2(\tilde{S}^{\top}\tilde{S}) - \mathcal{B}_{\infty}^2(\mathbb{E}[\tilde{S}^{\top}\tilde{S}])| + |\mathcal{B}_{\infty}^2(\mathbb{E}[\tilde{S}^{\top}\tilde{S}]) - \mathcal{B}_{\infty}^2(M_{\infty})| \\ & \le L \|\tilde{S}^{\top}\tilde{S} - \mathbb{E}[\tilde{S}^{\top}\tilde{S}]\|_F + L \|\mathbb{E}[\tilde{S}^{\top}\tilde{S}] - M_{\infty}\|_F \end{aligned} $$

Thus

$$ \begin{aligned} |\mathbb{E}[\mathcal{B}_{\infty}^2(\tilde{S}^{\top}\tilde{S})] - \mathcal{B}_{\infty}^2(M)| & \le L \mathbb{E} \left[ \|\tilde{S}^{\top}\tilde{S} - \mathbb{E}[\tilde{S}^{\top}\tilde{S}]\|_F \right] + L \|\mathbb{E}[\tilde{S}^{\top}\tilde{S}] - M\|_F \\ & \to 0. \end{aligned} $$

By plugging in $M_{\infty}$ and simplifying, we obtain the desired result.