Bulk Spectra of Truncated Sample Covariance Matrices
Subhroshekhar Ghosh Soumendu Sundar Mukherjee Himasish Talukdar


§ ABSTRACT 
Determinantal Point Processes (DPPs), which originate from quantum and statistical physics, are known for modelling diversity. Recent research <cit.> has demonstrated that certain matrix-valued U-statistics (that are truncated versions of the usual sample covariance matrix) can effectively estimate parameters in the context of Gaussian DPPs and enhance dimension reduction techniques, outperforming standard methods like PCA in clustering applications. This paper explores the spectral properties of these matrix-valued U-statistics in the null setting of an isotropic design. These matrices may be represented as X L X^⊤, where X is a data matrix and L is the Laplacian matrix of a random geometric graph associated to X. The main mathematically interesting twist here is that the matrix L is dependent on X. We give complete descriptions of the bulk spectra of these matrix-valued U-statistics in terms of the Stieltjes transforms of their empirical spectral measures. The results and the techniques are in fact able to address a broader class of kernelised random matrices, connecting their limiting spectra to generalised Marčenko-Pastur laws and free probability.


Design of a Standard-Compliant Real-Time Neural Receiver for 5G NR
    
Reinhard Wiesmayr^,∗, Sebastian Cammerer^†, Fayçal Aït Aoudia^†, Jakob Hoydis^†
Jakub Zakrzewski^†, and Alexander Keller^†

^†NVIDIA, ^ETH Zurich, contact: scammerer@nvidia.com

^∗Work done during an internship at NVIDIA.

This work has received financial support from the European Union under Grant Agreement 101096379 (CENTRIC). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission (granting authority). Neither the European Union nor the granting authority can be held responsible for them.



    July 2023
===========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





§ INTRODUCTION

The explosion of large-scale data, often referred to as “big data”, has transformed industries, research fields, and everyday life in the recent years. The phenomenon of massive scale data has called for new approaches to modelling and analysis. In particular, the question of diverse samples to enable a more parsimonious representation of data has led to connections with statistical physics, wherein models of strongly repulsive particle systems have been leveraged to augment the diverseness of features in machine learning procedures.

A key model in that respect is that of determinantal point processes or DPPs. A DPP is a probability distribution over subsets of a given ground set, such that the probability of a subset is proportional to the determinant of a kernel matrix corresponding to the subset. DPPs are known for their ability to model diversity, making them useful for selecting a set of items that are spread out over the feature space. Originating in quantum and statistical physics, DPPs have quickly grown to have an increasing impact as a significant component of a machine learning toolbox based on negative dependence.

A major parametric model of DPPs that has attracted attention in recent years is that of the Gaussian Determinantal Processes, abbrv. GDP <cit.>. In particular, it was shown in <cit.> that a certain matrix-valued statistic

    Σ̂ = 1/2n^2∑_1 ≤ i, j ≤ n(X_i - X_j≤ r)   (X_i - X_j)(X_i - X_j)^⊤,

where (X_i)_i=1^n ⊂ℝ^p are data points and r is a suitably chosen threshold, effectively performs parameter estimation in the GDP model.
Σ̂ may be viewed as a certain truncation of the sample covariance matrix (based on pairwise distances between data points), Further, it was demonstrated empirically in <cit.> that this matrix-valued test statistic can be leveraged as an ansatz to build dimension reduction tools that arguably  outperform standard PCA based methods, especially in the context of clustering applications.



An understanding of this latter phenomenon would require an understanding of the spectrum of the matrix in (<ref>). In this paper, we take a first step towards this by studying its bulk spectrum, modelling the data as i.i.d. centered random variables. This is the most fundamental and basic setting in which one first needs to understand the behaviour and properties of the matrix in (<ref>).

In fact, we are able to analyse a broader class of matrix-valued statistics which widely generalises Σ̂, by incorporating a general class of kernel functions K(X_i,X_j) in lieu of the distance based cutoff function (X_i - X_j≤ r). We perform a detailed analysis of the bulk spectrum of this broad class of kernelised random matrices and obtain a concrete description of their limiting spectral distributions as the size of the dataset n and the dimension p go to ∞ in a way such that p/n → c ∈ (0, ∞) (the so-called proportional asymptotics regime). In particular, in the smooth case, where the kernelised interaction is a suitably regular function of their mutual interaction, we can explicitly characterise the limiting spectral distribution as a parameterised family of Marčenko-Pastur laws (see Theorem <ref>). In the non-smooth case, we obtain a certain generalised Marčenko-Pastur law as the limiting spectral distribution (see Theorem <ref>). See Figure <ref> for histograms of the bulk spectra of Σ̂ for different choices of the threshold r. See Figure <ref> for an example of a kernelised version of Σ̂ with a smooth kernel. We also obtain the limiting spectral distribution in the semi-high-dimensional regime where p / n → 0 and p ≫√(n) (see Theorem <ref>).



The main insight that goes into analysing the spectrum of the matrix in (<ref>) is to represent it as a matrix-valued Rayleigh quotient:

    1/n^2 X L X^T,

where L is the Laplacian matrix of the random geometric graph on n vertices whose edges are given by (X_i - X_j≤ r). It is to be noted that matrices of the form X A X^⊤, where A is a positive semi-definite matrix independent of X have been studied in the literature in great detail. For example, in <cit.>, the authors considered the matrix 1/nXAX^⊤, where X is a n × p (note that the roles of n and p are reversed in their notation but this is only a cosmetic difference) matrix of i.i.d. entries with zero mean and unit variance and A is a diagonal matrix having some deterministic limiting spectral distribution μ_A. If X and A are independent and p/n→ y, then their result says that the above-mentioned matrix has a deterministic limiting spectral distribution, whose Stieltjes transform s is described as the unique solution, in the upper half plane ℂ^+ = {z ∈ℂ :  z > 0}, of the equation

    s(z) = 1/-z + y∫t d μ_A(t)/1+ts(z),

for z ∈ℂ^+. The resulting limiting spectral distribution is called a generalised Marčenko-Pastur Law which also admits the following free probabilistic interpretation: it is the free multiplicative convolution of μ_A and the Marčenko-Pastur law. In a more recent work, under certain additional assumptions, <cit.> obtained local laws for the matrices 1/nX A X^⊤ and 1/n A^1/2XX^⊤ A^1/2, where A is a deterministic matrix. The crucial difference of our model from these existing works is that the Laplacian matrix L is dependent on X. As such we need use careful decoupling arguments to analyse its spectrum. Obtaining local laws under our setting is an interesting direction for future research.

It may also be observed that Σ̂, and its kernelised generalizations, belong to the wider family of matrix-valued U-statistics. As such, our results also contribute to the burgeoning theory of matrix-valued U-statistics and their applications. For instance, the spectrum of a matrix-valued Kendall's τ statistic was studied recently by <cit.>. They showed that if X_1, X_2, …, X_n are i.i.d. p-dimensional random vectors with independent entries from a continuous distribution, then the empirical spectral distribution of the matrix-valued Kendall's τ statistic, defined as

    τ = 1/n2∑_1≤ i <j ≤ nsign(X_i - X_j)  sign(X_i - X_j)^⊤,

converges weakly to 1/3 + 2/3Y, in probability, where Y follows the standard Marčenko-Pastur distribution (here the sign function is applied componentwise). The proof heavily relies on a matrix version of the Hoeffding decomposition for U-statistics. Although, the matrix (<ref>) is also a matrix-valued U-statistic, the presence of the cutoff factor makes a direct use of Hoeffding decomposition difficult. Instead, we directly analyse the Stieltjes transform of the empirical spectral distribution.

The rest of the paper is organised as follows. In Section <ref>, we describe the model under consideration and recall preliminaries of random matrices. In Section <ref> we state our main results and work out some examples. We also provide brief proof sketches of our main results in this section. Section <ref> gives detailed proofs of all the results. Finally, in Appendix <ref>, we collect some useful results from matrix analysis and concentration of measure which are used throughout the paper.



§ THE MODEL

Suppose w_ij, i ∈ [p], j ∈ [n] are i.i.d. random variables on some probability space (Ω, ℱ, ℙ). Assume w_11 = 0,  (w_11) = σ^2 and w^4_11 < ∞. Define the p-dimensional vectors X_j= (w_1j, w_2j,… , w_pj)^⊤, j = 1, 2, …, n. X is the p × n matrix with X_j's as columns. Also define X̅ = 1/n∑_i = 1^n X_i.

We consider two asymptotic regimes:

    
  * The proportional asymptotics regime: p/n→ c ∈ (0, ∞).
    
  * The semi-high-dimensional regime: p/n→ 0.


Suppose that K_p:ℝ^p ×ℝ^p → [0,1] is a function symmetric in its coordinates, that is K_p (U, V) = K_p (V, U). Let A denote the n × n symmetric random matrix with entries A_ij = K_p (X_i, X_j).

The Empirical Spectral Distribution (ESD) of a real symmetric matrix Y_n × n is defined as

    μ_Y = 1/n∑_i = 1^nδ_λ_i,

where λ_1, λ_2, …, λ_n are the eigenvalues of Y. The weak limit of the ESD (defined almost surely or in probability depending on the context) is called the Limiting Spectral Distribution (LSD).

In this paper we are interested in the truncated covariance matrix

    M = 1/2n^2∑_1 ≤ i, j ≤ n A_ij (X_i - X_j)(X_i - X_j)^⊤,

which is a generalisation of the estimator in (<ref>).

Notice that if K ≡ 1, then

    M = 1/2n^2∑_1 ≤ i, j ≤ n (X_i - X_j)(X_i - X_j)^⊤ = 1/n∑_i = 1^n (X_i - X̅)(X_i - X̅)^⊤,

which is the sample-covariance matrix of the observations X_1, …, X_n. As it is a rank-1 perturbation of the matrix 1/n X X^⊤ (which will also be called the sample-covariance matrix), they share the same LSD.

In the proportional asymptotic regime, it is well known that the sample-covariance matrix 1/n X X^⊤ has as its LSD the Marčhenko-Pastur distribution _c, σ^2 with parameters (c, σ^2). Recall that when c ∈ (0, 1], _c, σ^2 has density

    d_c, σ^2(x) = 1/2 πσ^2√((b - x)(x - a))/cx𝕀_(a, b)dx,

where a = σ^2 (1 - √(c))^2 and b = σ^2(1 + √(c))^2. When c > 1, _c, σ^2 has a mass of (1 - 1/c) at 0, the remaining part has the same density as above, i.e.

    _c, σ^2 = (1 - 1/c) δ_0 + 1/cν,

where

    dν(x) = 1/2 πσ^2√((b - x)(x - a))/cx𝕀_(a, b)dx.

We will see that when the kernel K_p is Lipschitz and the entries w_ij satisfy some regularity conditions, then the LSD of M is a scaled Marčenko-Pastur law (see Theorem <ref>). However, if K_p is non-smooth, then a different LSD emerges, which is a generalised Marčenko-Pastur law (see Theorem <ref>).

On the other hand, in the semi-high-dimensional regime, one requires different scaling and centering. In is well known that the ESD of √(n/p)(1/nXX^⊤ - I) converges to the standard semi-circle law (see <cit.>, p. 864). The semi-circle law _ϖ^2 with variance ϖ^2 > 0 is defined as

    d_ϖ^2 (x) = 1/2 πϖ^2√(4ϖ^2 - x^2)  (|x| ≤ 2ϖ)   dx.

For ϖ = 1, we have the standard semi-circle law.

In our setup, we prove that under certain conditions on the moments of A_12, in the regime

    p ≫√(n),

the ESD of

    E = √(n/p)(M - α_pσ^2 I)

converges weakly to a semi-circle law with parameters depending on A_12^2 and σ^2, almost surely. We also describe the Stieltjes transform of the limiting distribution (see Theorem <ref>).



§ MAIN RESULTS



 §.§ The non-smooth case

Let d : ℝ^2 →ℝ be a symmetric function such that

    𝔼|d(w_11,w_12)|^3< ∞.

Typical examples of d(x,y) are (x-y)^2 or |x-y|. We define d_p:ℝ^p×ℝ^p→ℝ, d_p(x,y)=∑_i=1^p d(x_i,y_i). Let ϕ_p:ℝ→ [0,1] be monotonic and potentially dependent on p. For the purpose of the first theorem we shall assume that K_p has the following form

    K_p(x,y)=ϕ_p(d_p(x,y)).

Notice that this class of kernels includes the indicator kernel (x-y≤ r_p) and the Gaussian kernel 1- exp(-x-y^2/2τ^2_p), where r_p, τ_p are suitable constants.

We shall also require some limiting properties of the sequence of functions (ϕ_p). We state them now. First fix the following notations:

    m_1    = [d(X_11,X_12)]; 
    
        m_2    =Var(d(X_11,X_12)); 
    
        m_2^(1)   = Var(𝔼[d(X_11,X_12)|X_11]); 
    
        m_2^(2)   = 𝔼Var(d(X_11,X_12)|X_11).

Define the functions ψ_p,ϕ̃_p: ℝ→ℝ as ψ_p(x)=pm_1+√(pm_2)x and ϕ̃_p=ϕ_p∘ψ_p.

Suppose there exists ϕ̃ such that for any ϵ>0,

    (|ϕ̃_p-ϕ̃|>ϵ)→ 0

as p →∞ where  is the Lebesgue measure. In other words, ϕ̃_p converges in Lebesgue measure to ϕ̃.


Our first theorem describes the LSD of M in terms of its Stieltjes transform. Throughout the paper z will denote a complex number with u= z and v= z, i.e. z = u + ι v. Recall that the Stieltjes transform S_μ of a probability measure μ on ℝ is a complex function defined for z∈ℂ^+ as follows:

    S_μ(z) := ∫dμ(x)/x-z.

Suppose {μ_n}_n ≥ 1, μ are probability measures on ℝ with Stieltjes transforms {S_μ_n}_n ≥ 1 and S_μ, respectively. It is well known that S_μ_n→ S_μ pointwise on ℂ^+ if and only if μ_n →μ weakly (see, e.g., <cit.>). Moreover, if {μ_n}_n ≥ 1 are random probability measures and μ is a deterministic probability measure, then S_μ_n (z) → S_μ (z) almost surely for each fixed z ∈ℂ^+ if and only if μ_n →μ almost surely.

Suppose that Assumption <ref> holds. Then the ESD of M converges weakly to a deterministic distribution, almost surely. Moreover, if s(z) is the Stieltjes transform of the limiting distribution, then s(z) is the unique solution in ℂ^+ of the following equation:

    1 + zs(z) = 𝔼_ζ[σ^2 s(z) ζ/1+cσ^2 s(z) ζ],

where

    ζ = _Z_2[ϕ̃(√(m_2^(1)/m_2)Z_1 + √(m_2^(2)/m_2)Z_2) ],

with Z_1, Z_2 being i.i.d. N(0,1) random variables.


Suppose ξ_1, ξ_2, …, ξ_n are i.i.d. bounded random variables such that ξ_i and X_j are independent for all i ≠ j. Further assume that ξ_1 converges in distribution to some variable ζ. Then, the proof of Theorem <ref> will show that the LSD of 1/n∑_i=1^n ξ_i X_i X_i^⊤ is given by (<ref>). This result is known if we further assume that ξ_i is independent of X_i, but here we allow them to depend.


If M is represented as a matrix-valued Rayleigh quotient 1/n^2 X L X^⊤, our proof will show that the ESD of 1/nL will converge weakly to ζ as defined in (<ref>). Now, a moment's thought will reveal that the equations (<ref>) and (<ref>) are equivalent once we make the necessary adjustments for the scaling. In other words, even though X and L are dependent and L is not diagonal, a generalised Marčenko-Pastur law emerges as the LSD.

[Indicator kernel]
We first consider kernel

    K_p(x, y) = (x - y≤ r_p),

where r_p is an appropriate threshold. This gives us the estimator (<ref>) that motivated the present study. In this case, we will assume that w_11 is Gaussian. At the least, we shall need that α_p K_p(X_1, X_2) converges to a nonzero quantity as p →∞. The suitable choice for r_p turns out to be

    r_p^2= ((2p+2√(2p)z_α) σ^2+o(√(p))).

For α∈(0,1), z_α=Φ^-1(α), where Φ is the distribution function of the standard normal variable. Observe that X_1-X_2^2/2σ^2∼χ^2_p. Using the central limit theorem, U=X_1-X_2^2/2σ^2-p/√(2p) converges, in distribution, to a standard Gaussian. Notice that

    α_p   = ℙ(X_1-X_2≤ r_p) 
       =ℙ(X_1-X_2^2/2σ^2-p/√(2p)≤r_p^2/2σ^2-p/√(2p)) 
       = ℙ(U≤ z_α+ o(1)) 
       = α + o(1),

as p→∞. Now, in this case, ϕ_p (t) = I(t ≤ r_p^2) and d(x, y) = (x-y)^2. An easy computation shows that

    m_1 = 2 σ^2,       m_2 = 8 σ^4, 
    
        m_2^(1)= 2 σ^4,       m_2^(2) = 6 σ^4.

Then ϕ̃_p and ϕ̃ turn out to be as follows:

    ϕ̃_p(t)    = I( t ≤ z_α +o(1)), 
    ϕ̃ (t)    = I(t ≤ z_α).

Hence the distribution of ζ may be described as

    ζ   = 1/√(2 π)∫_ℝϕ̃(1/2Z+ √(3)/2t)e^-t^2/2dt 
       =  Φ(-1/√(3)Z +2/√(3) z_α)
       d=Φ(1/√(3)Z +2/√(3) z_α),

where Z ∼ N(0,1).

In the next example, however, Theorem <ref> can not be applied.
[Gaussian kernel]
Let us now consider the Gaussian kernel

    K_p(x,y) = 1- e^-x-y^2/2pτ^2,

where we have taken τ_p^2 =2 p τ^2. Here, ϕ_p (t) = 1 - e^-t/2pτ^2 and d(x,y) = (x-y)^2. Let us also assume that w_11 is Gaussian. One can calculate that

    ϕ̃_p(t) = 1-e^-σ^2/τ e^-√(2)σ^2 t/√(p)τ^2,

which does not satisfy Assumption <ref>.


In order to examine the Gaussian kernel, one needs to use the smoothness properties of the Gaussian kernel which was missing in our analysis. This leads us to the next theorem where the smoothness of the kernel is crucially used.


 §.§ The smooth case

We say that K_p is Lipschitz with Lipschitz constant κ_p (or κ_p-Lipschitz in short) if

    |K_p(x_1, y_1) - K_p(x_2, y_2)| ≤κ_p(x_1-x_2+y_1-y_2).

In particular, if for p ∈ℕ, ϕ_p : ℝ→ [0,1] is a κ_p-Lipschitz function and one takes K_p(x,y) = ϕ_p(x-y), then one can check that K_p is also an κ_p-Lipschitz kernel. Following <cit.>, we will make the following assumption on w_11.
[LC class property]
Let ω > 0. We require w_11 to satisfy either of the following three conditions:

    
  (a) w_11d=φ(Z) for some Lipschitz function φ with φ_Lip≤ω, where Z is a standard normal variable.
    
  (b) w_11 has density uniformly bounded below by 1/ω.
    
  (c) X_1 is strongly log-concave with curvature ≥ 1/ω^2.


Define α_p =  K(X_1, X_2).

Suppose that w_11 satisfies Assumption <ref> with parameter ω_p. Let K_p be an κ_p-Lipschitz kernel, such that κ_pω_p =o(1/√(log(n))) and α_p →α. Then, the ESD of M converges weakly to _c, α^2σ^2, almost surely.



Let us revisit Example <ref> in light of Theorem <ref>. Take w_11∼ N(0, σ^2). In this case, one may take κ_p = √(2)/e√(p)τ, α_ p = 1- (1 + 2σ^2/pτ^2)^-p/2 and α = 1 - e^-σ^2 / τ^2. Thus, by Theorem <ref>, μ_M _c, σ^2 (1 - e^-σ^2/τ^2)^2, almost surely.




 §.§ The semi-high-dimensional regime

In Theorem <ref> if we take c to be 0, then the Stieltjes transform of the LSD turns out to be 1/ασ^2 - z which corresponds to the measure δ_ασ^2. This shows that E is the right matrix to look at.
Define β_p ^2  A_12^2. Indeed, we show that

    If α_p →α, β_p^2 →β^2 and w_11^8 < ∞, in the regime p ≫√(n), the ESD of E converges weakly to _βσ^2, almost surely.



    Note that without additional assumptions on |α_p - α|, we can not determine the convergence of
    
    E' = √(n/p)(M - ασ^2 I),

    hence the appearance of α_p in (<ref>) instead of α. For instance, if  |α_p - α| = O(1/√(p)), then we can replace E by E' in Theorem <ref>.


In the set up of Example <ref>, β_p^2=α_p →α. Thus, the LSD of E, in this case is _√(α)σ^2. As for Example <ref>, β_p turns out to be 1 - 2(1-2 σ^2/p τ^2)^-p/2 + (1-4σ^2/pτ^2)^-p/2 and the parameter for  changes accordingly.

In the remainder of this section, we briefly sketch the proofs of Theorems <ref>,  <ref> and <ref>.



 §.§ Proof sketches

In this section, we give brief sketches of the proofs.


  §.§.§ The non-smooth case

We first observe that the matrix M can be written in the form 1/n^2XLX^⊤ where L is Laplacian matrix of a random geometric graph. Since the entries of L are bounded, it turns out to be sufficient to consider only 1/n^2 X D X^⊤, where D = (L). Note that D_ii = ∑_j≠ i K_p(X_i, X_j). This has rather high dependence on X_i and low dependence on the X_j's, for j ≠ i. Observe that conditional on X_i, { K_p(X_i, X_j) -  [K_p(X_i, X_j) | X_i]}_j≠ i are bounded, centered i.i.d. random variables. Using Hoeffding's inequality, we show that the sum ∑_j≠ i (K_p(X_i, X_j)-  [K_p(X_i, X_j) | X_i]) is negligible. It thus boils down to finding the LSD of the matrix

    M̅ = ∑_i=1^n ξ_i X_i X_i^⊤,

where the ξ_i's are i.i.d and independent of X_j, j≠ i. One can still not use the existing results in the literature since ξ_i is not independent of X_i. However, ξ_i involves a sum of certain i.i.d. variables. This enables us to apply a Berry-Esseen bound, eventually providing us with the distributional convergence of ξ_i. We now look at the Stieltjes transform S_M̅ of  M̅. Using standard martingale techniques, we deduce that S_M̅(z) - [S_M̅(z)] → 0 almost surely for each fixed z ∈^+. With this our goal becomes to get a recursive formula for [S_M̅]. This is done by using the Sherman-Morrison formula and various perturbation inequalities for matrices coupled with the fact that ξ_i converges in distribution as p →∞.



  §.§.§ The smooth case

As in the non-smooth case, we start by writing M in the Rayleigh quotient form 1/n^2 X L X^⊤. Our goal is to show that M has the same LSD as 1/n^2 X ℒ X^⊤ where ℒ =  L = α_n (nI-J), where J is the n× n matrix with all entries equal to 1. This is done by using matrix perturbation inequalities and Theorem 1 from <cit.>. Since J is of rank 1, we have a further simplification: we may just consider the matrix α_n/nX X^⊤, which clearly has _c, α^2 σ^2 as its LSD.



  §.§.§ The p/n → 0 case.

As in the proof of Theorem <ref>, we can show that the LSD of E is the same as that of

    E̅_1 = ∑_i=1^n ξ_i (X_i X_i^⊤ - σ^2 I).

We analyse this matrix via its Stieltjes transform and various matrix perturbation inequalities.



§ PROOFS



 §.§ Proof of Theorem  <ref>

Since zeroing out the diagonal entries of A does not change the Laplacian L = D - A, we may redefine A (with a slight abuse of notation) as follows:

    A    = (((1 - δ_ij)A_ij)), 
    
        D    = (A ), 
    
        L    = D - A,

where  is the vector with each entry equal to 1. Notice that L is the Laplacian matrix corresponding to the weighted adjacency matrix A. The basic observation that will help us find the LSD of M is the following matrix Rayleigh quotient representation:

    M = 1/n^2 X L X^⊤.

The idea is to divide M into three parts such that one part determines the LSD and the rest are negligible. In order to do so, we further decompose D into two diagonal matrices. Define ξ_i =  [K_p(X_i,V) | X_i], where V∼ X_1 is independent of X_1, X_2,… X_n. Notice that

    D_ii   = ∑_j≠ i K_p(X_i,X_j) 
       = (n-1)ξ_i+∑_j≠ i(K_p(X_i,X_j)- [K_p(X_i,X_j)|X_i]).

Define the following quantities:

    ξ    = (ξ_1,ξ_2, …, ξ_n)^⊤,
    ξ'_ij    = K_p(X_i,X_j)- [K_p(X_i,X_j)|X_i],
    ξ'_i        = ∑_j≠ iξ'_ij,
    ξ'          = (ξ'_1, ξ'_2,…,ξ'_n)^⊤,
    
        D_1           = (n-1)diag(ξ),
    
        D_2           = diag(ξ').

With these notations set, one has D = D_1 + D_2. Now decompose M as

    M=1/n^2 X D_1 X^⊤ + 1/n^2 X D_2 X^⊤ - 1/n^2 X A X^⊤.

Define M̃ = 1/n^2 X D_1 X^⊤. We will show that M has the same LSD as M̃. Towards that end, let d_W_2 denotes the 2-Wasserstein distance between probability measures μ_1 and μ_2 possessing finite second moments:

    d_W_2(μ_1,μ_2) := inf√( (Z_1 - Z_2)^2),

where the infimum is taken over all possible couplings of (Z_1, Z_2) with marginals Z_1 ∼μ_1 and Z_2 ∼μ_2.

    We have d_W_2(μ_M, μ_M̃)  0.


Using the Hoffman-Wielandt inequality (see Lemma <ref>) and the facts that

    G H_≤min{G_H_, G_H_},

and

    X^⊤_ = X_ = √(X^⊤ X_),

we have the following estimate

    d_W_2(μ_M, μ_M̃)    ≤1/√(p)M - M̃_
       = 1/n^2 √(p)X(D_2-A) X^⊤_
       ≤1/n^2 √(p)X D_2 X^⊤_ +1/n^2 √(p)X  A X^⊤_
       ≤1/n^2X D_2 X^⊤_+ 1/n^2 √(p)X^⊤ X_A_
       ≤1/n X^⊤ X_ ( 1/nD_2_+  1/n √(p)A_ ).

Since the fourth moment of the entries is finite, using Theorem 3.1 of <cit.>,

    1/nX^⊤ X_ (1+√(c))^2σ^2.

Let us now consider A_ and D_2_. Since each entry of A is bounded by 1, A_≤ n. Notice that D_2_ = max_i |ξ'_i|. Fix i∈[n]. Conditional on X_i, {ξ'_ij}_j≠ i are i.i.d. random variables. Moreover, |ξ'_i|≤ 1. Therefore by Hoeffding's inequality, for any t > 0,

    ℙ(|ξ'_i|>t|X_i)≤ e^-t^2/2(n-1).

Since the right hand side does not depend on X_i, the above bound also holds unconditionally. Now, the exchangeability of the ξ'_i's yields

    ℙ(1/nD_2_>√(6log n/n))    ≤ nℙ(|ξ_1'|>√(6nlog n))
       ≤ ne^-6nlog n /2(n-1)
       ≤1/n^2.

It follows by the Borel-Cantelli Lemma that 1/nD_2_ 0. We conclude that

    d_W_2(μ_M, μ_M̃)  0.

This completes the proof.

In fact, a further simplification is possible.

    Suppose that M̅ := 1/n∑_i=1^n ξ_i X_i X_i^⊤. Then d_W_2(μ_M̃, μ_M̅)  0.


Since M̃= n-1/n^2∑_i=1^n ξ_i X_i X_i^⊤, we have

    d_W_2(μ_M̃ , μ_M̅) ≤M̃ - M̅_≤1/n^2X X^⊤_ 0,

where for the second inequality we have used the fact that |ξ_i| ≤ 1 so that the matrix

    XX^⊤ - ∑_i = 1^n ξ_i X_i X_i^⊤ = ∑_i = 1^n (1 - ξ_i) X_i X_i^⊤

is positive semi-definite.


By virtue of Lemmas <ref> and <ref>, it is enough to find the LSD of M̅. Let

    S_n(z)1/p (M̅-zI)^-1

be the Stieltjes transform of M̅. We shall now show that S_n(z) converges almost surely as n→∞, for each fixed z ∈ℂ^+. The following lemma shows that it suffices to find the limit of S_n(z).


For each fixed z ∈ℂ^+,

    S_n(z) -  S_n(z)  0.



The proof uses the well-known martingale technique in random matrix theory (see, e.g., the proof of Theorem 3.10 of <cit.>). Define ℱ_0 = {ϕ, Ω} and ℱ_k = σ(X_1, X_2, …, X_k) for k ∈ [n]. By _k, k= 0, 1, …, n, we denote the conditional expectation operator given ℱ_k. Then

    S_n(z) -  S_n(z) = ∑_k=0^n (_k [1/p(M̅ - zI)^-1] - _k-1[1/p(M̅ - zI)^-1]).

Call the k-th summand above γ_k. Then {(γ_k, ℱ_k)}_k=1^n is a martingale difference sequence. Suppose

    M̅_k = M̅ - 1/nξ_k X_k X_k^⊤.

Notice that

    _k (M̅_k - zI)^-1  = _k-1(M̅_k - zI)^-1.

By Lemma <ref>(b),

    |(M̅ - zI)^-1 - (M̅_k - zI)^-1| ≤1/v.

Define S_nk=1/p(M̅_k - zI)^-1. Therefore

    |γ_k| = |_k [S_n(z) - S_nk(z)] -_k-1 [S_n(z) - S_nk(z)]| ≤ 2/v.

Thus γ_k is a bounded martingale difference sequence. By Lemma <ref>,

    |S_n(z) -  S_n(z)|^4 ≤K_4/p^4( ∑_k=1^n |γ_k|^2 )^2 ≤4 K_4 n^2/v^4 p^4 = O(n^-2).

Now an application of the Borel-Cantelli lemma gives us (<ref>).



    Suppose V, W are i.i.d. copies of X_1. Define
    
    ζ_p= 𝔼[ϕ_p(d(V,W))|V].

    Then ζ_p ζ, where
    
    ζ = _Z_2[ϕ̃(√(m_2^(1)/m_2)Z_1 + √(m_2^(2)/m_2)Z_2) ],

    where Z_1, Z_2 are i.i.d. N(0,1) random variables.



Write V=(V_1,V_2,…,V_p)^⊤ and W=(W_1,W_2,…,W_p)^⊤. Let us fix the following notations:

    m_1'(v)    = 𝔼[d(V_1,W_1)|V_1=v], 
    
        m_2'(v)   = Var(d(V_1,W_1)|V_1=v),
    
        m_3'(v)    = 𝔼[|d(V_1,W_1)-m_1'(V_1)|^3|V_1=v].

Define

    T=∑_i=1^p d(V_i,W_i)-∑_i=1^p m_1'(V_i)/√(∑_i=1^p m_2'(V_i)).

Applying the Berry-Esseen theorem (see, e.g., <cit.>) on {d(V_i,W_i)}_i=1^p, conditional on V,

    sup_x|ℙ(T≤ x|V)-Φ(x)|≤ C_1 Q/√(p),

where C_1 is an absolute constant and

    Q = (1/p∑ _i=1^p m_2'(V_i))^-3/2(1/p∑_i=1^p m_3'(V_i))^1/2.

By the SLLN,

    Q m_2^(2)^-3/2 (m_3)^1/2.

For notational convenience, let us define

    a_p(V)   =∑_i=1^p m_1'(V_i),
    
    b_p(V)   = √(∑_i=1^p m_2'(V_i)).

Suppose Z∼ N(0,1) is independent of V and W. Without loss of generality assume that ϕ_p is increasing. Then,

    |𝔼[ϕ_p    (d(V,W))|V] - 𝔼[ϕ_p(a_p(V)+b_p(V)Z)|V]|
       = |∫_0^1ℙ(ϕ_p(a_p(V)+b_p(V)T)>t|V) dt -∫_0^1ℙ(ϕ_p(a_p(V)+b_p(V)Z)>t|V) dt|
       ≤∫_0^1 |ℙ(T>ϕ_p^-1(t)-a_p(V)/b_p(V)|V)-ℙ(Z>ϕ_p^-1(t)-a_p(V)/b_p(V)|V)| dt
       ≤ C_1 Q/√(p)→ 0

almost surely. Notice that

    𝔼[ϕ_p(a_p(V)+b_p(V)Z)|V]=1/√(2 π)∫_ℝϕ_p(a_p(V)+b_p(V)t)e^-t^2/2dt.

From the relation between ϕ and ϕ̃, we have

    ϕ_p(a(V)+b(V)t)    =ϕ̃_p∘ψ_p^-1(a(V)+b(V)t)
       = ϕ̃_p (a(V)+b(V)t-pm_1/√(pm_2))  
       = ϕ̃_p( a_p'(V)+b_p'(V)t),

where

    a_p'(V)   = a(V)-pm_1/√(pm_2),
    
        b_p'(V)    = b(V)/√(pm_2).

Using the central limit theorem on {m_1'(V_i)}_i=1^p, as p→∞, we have

    a_p(V)-pm_1/√(pm_2^(1)) N(0,1),

and consequently,

    a_p'(V)  N(0,m_2^(1)/m_2).

On the other hand, an application of the SLLN shows that

    b_p'(V) √(m_2^(2)/m_2).

We are interested in the weak limit of

    1/√(2 π)∫_ℝϕ̃_p(a_p'(V)+b_p'(V)t)e^-t^2/2dt

as p →∞. Suppose (by an abuse of notation) that a_p', b_p' are sequences of real numbers. Fix ϵ>0 and define the event

    G = {|ϕ̃_p(a_p'+b_p't)dt - ϕ̃(a_p'+b_p't)|≤ϵ}.

Then

    ∫_ℝ|ϕ̃_p    (a_p'+b_p't)dt - ϕ̃(a_p'+b_p't)|e^-t^2/2dt
       ≤∫_G|ϕ̃_p(a_p'+b_p't)dt - ϕ̃(a_p'+b_p't)|e^-t^2/2dt+∫_G^c|ϕ̃_p(a_p'+b_p't)dt - ϕ̃(a_p'+b_p't)|e^-t^2/2dt.

The first integral above is bounded by √(2π)ϵ and the second one by (|ϕ̃_p-ϕ̃|>ϵ). Taking a supremum over a_p', b_p', we get

    sup_a_p', b_p'∫_ℝ|ϕ̃_p(a_p'+b_p't)dt - ϕ̃(a_p'+b_p't)|e^-t^2/2dt ≤√(2π)ϵ + (|ϕ̃_p-ϕ̃|>ϵ).

Now, sending p→∞ and noticing that ϵ is arbitrary, we obtain

    sup_a_p', b_p'∫_ℝ|ϕ̃_p(a_p'+b_p't)dt - ϕ̃(a_p'+b_p't)|e^-t^2/2dt → 0.

as p→∞. It is enough to consider the weak limit of

    1/√(2 π)∫_ℝϕ̃_p(a_p'(V)+b_p'(V)t)e^-t^2/2dt

as p→∞.

Next we prove that if a_p'→ a' and b_p' → b' (both deterministic sequences) as p →∞, then

    ∫_ℝϕ̃(a_p'+b_p't)e^-t^2/2dt →∫_ℝϕ̃(a'+b't)e^-t^2/2dt,

in other words, the map (a',b')↦∫_ℝϕ̃(a'+b't)e^-t^2/2 is continuous. Since ϕ̃ is monotonic, it is continuous almost everywhere. Thus, ϕ̃(a_p+b_pt)→ϕ̃(a+bt) for almost every t. An application of the bounded convergence theorem then proves our claim. Combining this with the convergence of a'(V) and b'(V), we get

    1/√(2π)∫_ℝϕ̃(a_p'(V)+b_p(V)'t)e^-t^2/2dt d⟶1/√(2 π)∫_ℝϕ̃(√(m_2^(1)/m_2)Z+ √(m_2^(2)/m_2)t)e^-t^2/2dt,

where Z∼ N(0,1). This completes the proof.



𝔼 S_n(z) satisfies the following approximate recursion:

    1 + z𝔼S_n(z) = 𝔼_ζ[σ^2ζ𝔼S_n(z)/1+c_nσ^2ζ𝔼S_n(z)] + o(1).



    Suppose ξ_0∼ξ_1, X_0∼ X_1 are mutually independent and they are independent of the {X_i}_i=1^n. Define

    M̅'   =M̅+1/nξ_0X_0X_0^⊤,
    
        S_n'(z)   = 1/p (M̅'-zI)^-1.

Then

    p    =  (M̅'-zI)  (M̅'-zI)^-1
       =1/n∑_i=0^nξ_i X_i^⊤ (M̅'-zI)^-1 X_i- z  (M̅'-zI)^-1.

Taking expectation on both sides and using the exchangeability of the summands, we have

    p   =(1+1/n) 𝔼(ξ_0 X_0^⊤ (M̅'-zI)^-1X_0) - pz𝔼S_n'(z).

Dividing by n and noticing that |ξ_0 X_0^⊤ (M̅'-zI)^-1X_0|≤ 1+|z|/v, we get

    c_n+ c_nz𝔼S'_n(z)=1/n𝔼(ξ_0 X_0^⊤ (M̅'-zI)^-1X_0)+o(1).

Observe that

    |S'_n(z)-S_n(z)|≤v/p.

Therefore

    c_n+c_nz𝔼S_n(z)=1/n𝔼(ξ_0X_0^⊤ (M̅'-zI)^-1X_0) +o(1).

An application of the Sherman-Morrison formula on the first term yields

    1/nξ_0X_0^⊤ (M̅'-zI)^-1X_0= 1/nξ_0X_0^⊤ (M̅-zI)^-1X_0/1+1/nξ_0X_0^⊤ (M̅-zI)^-1X_0.

Then

    |1/nξ_0    X_0^⊤ (M̅'-zI)^-1X_0-c_nσ^2ξ_0𝔼S_n(z)/1+c_nσ^2ξ_0𝔼S_n(z)|
       =|1/nξ_0X_0^⊤ (M̅-zI)^-1X_0/1+1/nξ_0X_0^⊤ (M̅-zI)^-1X_0-1/nξ_0𝔼(X_0^⊤ (M̅-zI)^-1X_0)/1+1/nξ_0𝔼(X_0^⊤ (M̅-zI)^-1X_0)|
       ≤1/nξ_0|X_0^⊤ (M̅-zI)^-1X_0-𝔼(X_0^⊤ (M̅-zI)^-1X_0)|/|1+1/nξ_0X_0^⊤ (M̅-zI)^-1X_0||1+1/nξ_0𝔼(X_0^⊤ (M̅-zI)^-1X_0)|
       ≤|z|^2/n|X_0^⊤ (M̅-zI)^-1X_0-𝔼(X_0^⊤ (M̅-zI)^-1X_0)|/|z+z1/nξ_0X_0^⊤ (M̅-zI)^-1X_0||z+z1/nξ_0𝔼(X_0^⊤ (M̅-zI)^-1X_0)|
       ≤|z|^2/v^21/n|X_0^⊤ (M̅-zI)^-1X_0-𝔼(X_0^⊤ (M̅-zI)^-1X_0)|,

where we have used Lemma <ref>(h) to lower bound the denominator. Since X_0 is independent of (M̅ -zI)^-1, by Lemma <ref>,

    |X_0^⊤ (M̅-zI)^-1X_0-𝔼(X_0^⊤ (M̅-zI)^-1X_0)| ≤ C√((M̅ -zI)^-1)≤C √(p)/√(v),

for some constant C > 0, which only depends on the second and fourth moments of w_11.

Thus

    𝔼1/nξ_0X_0^⊤ (M̅'-zI)^-1X_0=𝔼c_nσ^2ξ_0𝔼S_n(z)/1+c_nσ^2ξ_0𝔼S_n(z)+o(1).

Also,

    |𝔼c_nσ^2ξ_0𝔼S_n(z)/1+c_nσ^2ξ_0𝔼S_n(z)   -𝔼c_nσ^2ζ𝔼S_n(z)/1+c_nσ^2ζ𝔼S_n(z)|
       ≤σ^2 c_n|z|^2|𝔼S_n(z)| 𝔼|ξ_0-ζ|/|z+zc_nσ^2ξ_0𝔼S_n(z)||z+zc_nσ^2ζ_0𝔼S_n(z)|
       ≤σ^2 c_n |z|^2/v^21/v𝔼|ξ_0-ζ|.

Lemma <ref> coupled with the Skorohod representation theorem gives us that 𝔼|ξ_0-ζ|→ 0 as p→∞. Combining everything, we get the desired approximate functional for S_n(z):

    c_n+c_nz𝔼S_n(z)=𝔼_ζ[c_nσ^2ζ𝔼S_n(z)/1+c_nσ^2ζ𝔼S_n(z)] + o(1).

This completes the proof.



The equation

    1+zs(z)=𝔼_ζ[σ^2ζ s(z)/1+cσ^2ζ s(z)],

z∈ℂ^+, has a unique solution for s(z) in ℂ^+, where ζ is defined in (<ref>).


Suppose (<ref>) has two distinct solutions s_1, s_2 in ℂ^+. Fix some z∈ℂ^2 such that s_1(z) ≠ s_2(z). Note for j=1, 2,

    s_j = 1/-z + σ^2 ζ/1 + c σ^2 s_j ζ.

Notice

    Im(s_j) = v - Im(σ^2 ζ/1 +c σ^2 s_j ζ)/|-z + σ^2 ζ/1 + c σ^2 s_j ζ|^2 < [ c σ^4 Im(s_j) ζ^2/|1 + c σ^2 s_j ζ|^2]/|-z + σ^2 ζ/1 + c σ^2 s_j ζ|^2.

So,

    [ c σ^4 ζ^2/|1 + c σ^2 s_j ζ|^2]/|-z + σ^2 ζ/1 + c σ^2 s_j ζ|^2 > 1.

Also,

    s_1 - s_2 = c σ^4 (s_1 - s_2) ζ^2/(1 + c σ^2 s_1 ζ)(1 + c σ^2 s_2 ζ)/(-z + σ^2 ζ/1 + c σ^2 s_1 ζ)(-z + σ^2 ζ/1 + c σ^2 s_2 ζ).

Cancelling s_1 - s_2 from both sides and applying the Cauchy-Schwarz inequality on the right hand side we get

    1 ≤[ c σ^4 ζ^2/|1 + c σ^2 s_1 ζ|^2]/|-z + σ^2 ζ/1 + c σ^2 s_1 ζ|^2[ c σ^4 ζ^2/|1 + c σ^2 s_2 ζ|^2]/|-z + σ^2 ζ/1 + c σ^2 s_2 ζ|^2,

which contradicts (<ref>).


Now we are ready to prove Theorem <ref>.

Applying Lemmas <ref> and <ref>, it is enough to consider M̅ instead of M.  Since |𝔼S_n(z)|≤1/v, using the Bolzano-Weierstrass theorem, 𝔼S_n(z) has a convergent subsequence. Consider any such subsequence |𝔼S_n_k(z)|→ s(z). Then, by Lemma <ref>, s(z) will satisfy the limiting equation

    1+zs(z)=𝔼σ^2ζ s(z)/1+cσ^2ζ s(z).

Lemma <ref> shows that this equation has a unique solution for s(z) in ℂ^+.

This shows that all the subsequential limits of 𝔼S_n(z) are same proving that 𝔼S_n(s)→ s(z) as n→∞ where s(z) is the unique solution in ℂ^+ of the aforementioned equation. We can thus conclude that, almost surely, μ_M̅ converges weakly to some sub-probability measure μ with Stieltjes transform s. In order to show that μ is indeed a probability measure, by Prohorov's theorem it is sufficient to show that ∫ x^2 𝔼μ_M̅(dx) is uniformly bounded. By the trace-moment formula, the last quantity is the same as 1/p𝔼(M̅^2) which we compute next.

    1/p𝔼(M̅^2)    = 1/pn^2( ∑_i=1^n ξ_i^2 X_iX_i^⊤ X_i X_i^⊤ +∑_i ≠ jξ_i ξ_j X_i X_i^⊤ X_j X_j ^⊤)
       = 1/p n^2 (nX_1^4 + n(n-1)  (X_1^⊤ X_2)^2)
       = O(1),

where the last step follows from the fact that X_1^4 = O(p^2) and (X_1^⊤ X^2)^2 = (X_1^⊤ X_2) = O(p). This completes the proof of Theorem <ref>.




 §.§ Proof of Theorem <ref>

We again begin with the matrix Rayleigh quotient representation:

    M = 1/n^2 X^⊤ L X^⊤.


Observe that 𝔼[A] = α_n, p (J - I) and 𝔼[D]= (n - 1) α_n, p I, where J = 11^⊤. Define

    = [L] = [D] - [A] = α_n, p ((n - 1) I - (J - I)) = α_n, p (nI - J).


Further, let

    F    =  1/n (A - α_n, p (J - I)), 
    
        G    = 1/n (D - (n - 1)α_n, p I) = (F ).

With the above notation, 1/nL = 1/n - F + G. Thus

    M = 1/n^2 X L X^⊤   = 1/n X (1/n + F - G) X^⊤
       = 1/n^2 X X^⊤ + 1/n X (F - G) X^⊤
       = M̃ + 1/n X (F - G) X^⊤,

where M̃= 1/n^2 X  X^⊤. We will show that the M and M̃ have the same LSD. Indeed,

    d_W_2(μ_M, μ_M̃)    ≤1/√(n)M - M̃_
       = 1/n √(n)X (F - G) X^⊤_
       ≤1/n √(n)X_F - G_F X^⊤_
       = 1/nX X^⊤_1/√(n)F - G_
       ≤1/nX X^⊤_1/√(n)(F_ + G_) 
       ≤1/nX X^⊤_ (F_ + 1/√(n)F1) 
       ≤2/nX X^⊤_F_.


Under the finite fourth moment condition, recall that 1/nX X ^⊤_→σ^2(1 + √(c))^2 almost surely. We now show that F_ 0. Because of Assumption <ref>, we can apply Theorem 1 of <cit.> on A to get that

    (F_≥ 2 L ωσ(C + t/√(n))) ≤exp(-t^2 / C^2)

for some C > 0. We can simplify this as follows:

    (F_≥ 2 L ωσ t) ≤exp(- C'^2 t^2 ).

for some C'>0. Now choosing t_n= 1/C'√(log n), we note that

    (F_≥2 L ωσ√(log n)/C') ≤1/n^2.

Combining this with the fact that L ω = o(1/√(log n)), we conclude that F_→ 0 almost surely. This implies that d_W_2(μ_M, μ_M̃)→ 0 almost surely, that is M and M̃ have the weak limit almost surely if it exists. Notice that ℒ = α_pn I - α_p J. Since, J is of rank 1, by Lemma <ref>,

    d_KS(μ_M̃,μ_M̂) ≤1/n

where d_KS denotes the Kolmogorov-Smirnov distance and M̂ = α_p/n XX^⊤. Since

    d_W_2(μ_M̂, μ_α/nX X^⊤) ≤|α_p-α|/nX X^⊤_ 0,

we can only consider α_p/nX X^⊤. It is well known that the last matrix has _c, α^2σ^2 as its LSD. Combining everything, we conclude that

    μ_M _c, α^2σ^2,

almost surely. This completes the proof.



 §.§ Proof of Theorem <ref>


We define A, D, L as in the proof of Theorem <ref> and

    Ẽ√(n/p)(n-1/n^2∑_i=1^n ξ_i X_i X_i^⊤ - α_p σ^2 I),    E̅   √(n/p)(1/n∑_i=1^n ξ_i X_i X_i^⊤ - α_p σ^2 I).

One can prove analogues of Lemmas <ref> and <ref> in this set up so that d_W_2(μ_E, μ_Ẽ)  0 and d_W_2(μ_Ẽ , μ_E̅)  0, respectively (the analogue of the Bai-Yin result in this regime can be found in <cit.>). For the purpose of finding the LSD, it is thus enough to consider E̅. Now, E̅ can be further decomposed as E̅_1 + E̅_2, where

    E̅_1 = 1/√(np)∑_i=1^n ξ_i (X_i X_i^⊤ - σ^2 I),     and    E̅_2 = σ^2/√(np)∑_i=1^n (ξ_i - α_p) I.

Since the (ξ_i - α_p)'s are centered, independent and bounded, invoking Hoeffding's inequality, E̅_2_≤σ^2/√(p) w.p. ≥ 1 - 2e^-√(p)/2, implying that d_W_2 (μ_E̅, μ_E̅_1)  0. In order to find the LSD of E̅_1, we proceed with the Stieltjes transform once again. We first show (<ref>). Our technique remains same for this part, except that we need analogues of equations (<ref>) and (<ref>). First we define

    S̃_n(z) = 1/p( E̅_1 - zI)^-1, 
       E̅_1k = E̅_1 -1/√(np)ξ_k (X_k X_k^⊤ - σ^2 I).

By Lemma <ref> (b) and (c),

    |(E̅_1 - zI)^-1 -(E̅_1k - zI)^-1| ≤p + X_k^2/v^2 √(np).

Now, by virtue of Lemma <ref> for l=4,

    |S̃_n(z) - S̃_n(z)|^4    ≤K_4/p^4( ∑_k=1^n (2(p + X_k^2)/√(np))^2 )^2 
       ≤16 K_4 n^2/v^4 n^2 p^6( ∑_k=1^n (2p^2 + 2X_k^4))^2 
       ≤64 K_4 n^2/v^4 n^2 p^6(np^2 + ∑_k=1^n X_k^4)^2 
       ≤128 K_4 n^2/v^4 n^2 p^6(n^2p^4 + (∑_k=1^n X_k^4)^2) 
       ≤128 K_4 n^2/v^4 n^2 p^6(n^2p^4 + nX_1^8 + n(n-1) (X_1^4)^2 ) 
       =O(1/p^2).

Here the last line follows from the facts that X_1^8 = O(p^4) and X_1^4 = O(p^2). Now, by Borel-Cantelli lemma, (<ref>) follows. So, it is enough to consider S̃_n(z). Define s̃_n(z) = S̃_n(z). By Lemma <ref>, we get

    σ^4 β_n^2 s̃_n^2(z) + z s̃_n(z) + 1 = o(1),

which has the solutions

    s̃_n(z) = -z ±√(z^2 - 4σ^4 β_n^2 + o(1))/2σ^4 β_n^2.

Here √(z'), for any z' ∈ℂ\ℝ, denotes the square root in the upper-half plane. Clearly, one must take the + sign in  (<ref>), otherwise the right hand side has negative imaginary part which is not allowed. Now taking limit as n →∞, we get

    lim_n →∞s̃_n(z) = -z ±√(z^2 - 4σ^4 β^2)/2σ^4 β^2,

which the Stieltjes transform of the semi-circle law with variance β^2 σ^4. This completes the proof of Theorem <ref>.


   Let S̃_n(z) and s̃_n(z) be defined as in Theorem <ref>. Then, s̃_n(z) satisfies the following approximate functional equation
   
    σ^4 β_n^2 s̃_n^2(z) + z s̃_n(z) + 1 = o(1).



To this end, take an independent copy (ξ_0, X_0) of (ξ_1, X_1). Define 
    E̅'_1 = 1/√(np)∑_i=0^n ξ_i (X_i X_i^⊤ - σ^2 I),
            S'_n = 1/p( E̅'_1 - zI)^-1.

Proceeding as in the proof of Theorem <ref>, we can get

    1 + z s̃'_n(z) = n+1/p√(np)ξ_0 (X_0^⊤ (E̅'_1 - zI)^-1 X_0 - σ^2 (E̅'_1 - zI)^-1),

where s̃'_n(z) = S̃'_n(z). Since E̅'_1_≤1/v,

    |ξ_0 (X_0^⊤ (E̅'_1 - zI)^-1 X_0 - σ^2 (E̅'_1 - zI)^-1)| ≤1/vX_0^2 + p σ^2 /v≤2p σ^2 /v.

Thus, for the asymptotic purposes one can change the n+1 by n in the right hand side of (<ref>). By Lemma <ref>(b), (c),

    |(E̅'_1 - zI)^-1 - (E̅_1 - zI)^ -1| ≤1/v^2(1/√(np)X_0^2 + σ^2 √(p/n)),

where the expectation of the right hand side is 2 σ^2/v^2√(p/n). This also implies that |s̃'_n(z) - s̃_n(z)| ≤2 σ^2/v^21/√(np). Now, applying Lemma <ref>(b),

    |X_0^⊤ (E̅'_1 - zI)^-1 X_0 - X_0^⊤(E̅_1 + 1/√(np)ξ_0 X_0 X_0^⊤ - zI )^-1 X_0 | ≤σ^2/v √(np)X_0^2,

where the right hand side has expectation σ^2/v√(p/n). Using these perturbation results, we have the following upgrade of (<ref>):

    1 + zs̃_n(z) = √(n)/p√(p)ξ_0 (X_0^⊤( E̅_1 + 1/√(np)ξ_0 X_0 X_0^⊤ - zI )^-1 X_0 - σ^2 (E̅_1 - zI)^-1) + o(1).

We apply Sherman-Morrison formula (see Lemma <ref>(i)) on the first term in the right hand side to get

    1/√(np)ξ_0 X_0^⊤( E̅_1 + 1/√(np)ξ_0 X_0 X_0^⊤ - zI )^-1 X_0 = 1/√(np)ξ_0 X_0^⊤( E̅_1 - zI )^-1 X_0/1 + 1/√(np)ξ_0 X_0^⊤( E̅_1 - zI )^-1 X_0.

Notice we can write (<ref>) in the form 1 + zs̃_n(z) = -  B_1 +  B'_1 +  B_2 -  B'_2 +o(1), where

    B_1    = σ^2/p^2ξ_0^2 (E̅_1 - zI)^-1 X_0^⊤ ( E̅_1 - zI )^-1 X_0,
    
         B'_1    = B_1 1/√(np)ξ_0 X_0^⊤ ( E̅_1 - zI )^-1 X_0/1 + 1/√(np)ξ_0 X_0^⊤ ( E̅_1 - zI )^-1 X_0, 
    
         B_2    = √(n)/p √(p)ξ_0 R,
    
         B'_2    = 1/p^2ξ_0^3 R X_0^⊤ (E̅_1 -zI)^-1 X_0/1 + 1/√(np)ξ_0 X_0^⊤ (E̅_1 -zI)^-1 X_0, 
    
         R    = X_0^⊤ (E̅_1 - zI)^-1 X_0 - σ^2 (E̅_1 - zI)^-1.

We show that among the four, only B_1 contributes to the equation. Note

    B_1 = σ^4/p^2ξ_0^2 ((E̅_1 - zI)^-1)^2 + σ^2/p^2ξ_0^2 R(E̅_1 - zI)^-1 = σ^4 ξ_0^2 S̃_n^2(z)+σ^4 S̃_n(z) R/p.

By using Lemma <ref>,

    |R|^2 ≤ C_2 ((E̅_1 - zI)^-2) ≤C_2 p/v^2,

where C_2 only depends on the fourth moment of w_11. This additionally shows that R= O_P(1). Also, (S̃_n(z)) → 0 because, S̃_n(z) is a bounded sequence that converges almost surely. Thus,

    B_1 = σ^4 (ξ_0^2 S̃_n^2(z)) + o(1) = σ^4 ξ_0^2 S̃_n^2(z) +o(1) =σ^4 β_n^2 s̃_n^2(z) + o(1).

By (<ref>) and (<ref>), B_2 = o(1). Notice by Lemma <ref>(h), Im(z + z/√(np)ξ_0 X_0^⊤ (E̅_1 -zI)^-1 X_0) ≥ v, |X_0^⊤ (E̅_1 -zI)^-1 X_0| ≤X_0^2/v and (E̅_1 -zI)^-1≤ p/v. So,

    |B'_1| ≤σ^2 |z| /v^3 p √(np)X_0^4 = σ^2 |z|/v^3 p √(np) (p  w_11^4 +p(p-1)σ^4) = o(1),
        |B'_2| ≤|z|/v^2 p^2 [|R| X_0^2] ≤|z|/v^2 p^2 ( R^2)^1/2 (X_0^4)^1/2 = o(1).

Now, (<ref>) can be written as

    σ^4 β_n^2 s̃_n^2(z) + z s̃_n(z) + 1 = o(1)

completing the proof of Lemma <ref>




§ ACKNOWLEDGEMENTS

SG was supported in part by the MOE grants R-146-000-250-133, R-146-000-312-114, A-8002014-00-00 and MOE-T2EP20121-0013. SSM was partially supported by the INSPIRE research grant DST/INSPIRE/04/2018/002193 from the Dept. of Science and Technology, Govt. of India, and a Start-Up Grant from Indian Statistical Institute.

abbrvnat





§ AUXILIARY RESULTS

Here we collect lemmas and results borrowed from the literature. First we define some notations.

    Mat_n(ℂ) := The set of all  n × n  matrices with complex entries.

For A ∈Mat_n(ℂ), define the Hilbert-Schmidt norm of A by

    A_ := √(∑_1 ≤ i ≤ j ≤ n|A_ij|^2).

For x ∈ℝ^n, let x = √(∑_i = 1^n x^2_i). The Operator norm of A is defined as

    A_ := sup_x = 1Ax.

For a matrix A with eigenvalues λ_1, …, λ_n, let F_A(x) := 1/n∑_i = 1^n 𝕀_[λ_i, ∞)(x) be the empirical distribution function associated with the eigenvalues.  Let 𝕊_n denotes the set of all permutations of the set {1, 2, …, n}.

Let A,B ∈Mat_n(ℂ) are two normal matrices, with eigenvalues λ_1(A),λ_2(A), …, λ_n(A) and λ_1(B),λ_2(B), …, λ_n(B) respectively. Then we have

    min_σ∈𝕊_n∑_i = 1^n | λ_i(A) - λ_σ(i)(B)|^2 ≤A - B ^2_.

An immediate consequence of this is that

    d_W_2(μ_A, μ_B)^2 ≤A - B^2/n.



Let A, B ∈Mat_n(ℂ) are two Hermitian matrices. Then

    sup_x ∈ℝ |F_A(x) - F_B(x)| ≤rank(A - B)/n.



Let C ∈Mat_p(ℝ) be a symmetric, positive semi-definite matrix, y ∈ℝ^p, z = u+ι v ∈ℂ^+, ε >0, then

    
  (a) (C - zI)^-1_≤ 1/v;
    
  (b) |(C + y y^⊤ -zI)^-1 - (C - zI)^-1| ≤min{1/v, y^2/v^2};
    
  (c) |(C + ε I -zI)^-1 - (C -zI)^-1| ≤p ε/v^2;
    
  (d) (C + ε I -zI)^-1 - (C -zI)^-1_≤ε/v^2;
    
  (e) |y^⊤ (C+y y^⊤ -zI)^-1 y| ≤ 1 + |z|/v;
    
  (f) (z(C- zI)^-1) ≥ 0;
    
  (g) ((C - zI)^-1) >0;
    
  (h) (z y^⊤ (C- zI)^-1y) ≥ 0;
    
  (i) y^⊤ (c+ y y^⊤ - zI)^-1 y = y^⊤ (C-zI)^-1 y/1 + y^⊤ (C-zI)^-1 y.


In Lemma <ref>, parts (a) and (e)-(i) can be found in <cit.>. Parts (b)-(d) follow from the matrix identity A^-1 -B^-1 = A^-1(B-A)B^-1 provided A and B are invertible. We also state the following result for quadratic form of random vectors <cit.>.

Let A ∈Mat_p(ℂ) be non-random and Y = (Y_1, Y_2, …, Y_p) ∈ℂ^p be a random vectors of independent entries. Suppose Y_i = 0, |Y_i|^2 = 1 and |Y_i|^l≤ν_l for l= 3, 4, …, L, for some L ≥ 4. Then, for all l ≤ L/2, there exists C_l>0 (depends only on l), such that

    |Y^* A Y - A|^l ≤ C_l((ν_4 (AA^*))^l/2 + ν_2l ( (AA^*))^l/2).



Let X_k be a complex martingale difference sequence with respect to some filtration {ℱ_k}_k ≥ 1. Then, for l ≥ 1,

    |∑_k=1^n X_k|^l ≤ K_p (∑_k=1^n |X_k|^2 )^l/2,

whenever sup_k  |X_k|^l < ∞.

