# Asymptotic Gaussian Fluctuations of Eigenvectors in Spectral Clustering

Hugo Lebeau, Florent Chatelain, Romain Couillet

**Abstract**—The performance of spectral clustering relies on the fluctuations of the entries of the eigenvectors of a similarity matrix, which has been left uncharacterized until now. In this letter, it is shown that the signal + noise structure of a general spike random matrix model is transferred to the eigenvectors of the corresponding Gram kernel matrix and the fluctuations of their entries are Gaussian in the large-dimensional regime. This CLT-like result was the last missing piece to precisely predict the classification performance of spectral clustering. The proposed proof is very general and relies solely on the rotational invariance of the noise. Numerical experiments on synthetic and real data illustrate the universality of this phenomenon.

**Index Terms**—Spectral clustering, central limit theorem, kernel matrix, spike eigvector, Gaussian fluctuations.

## I. INTRODUCTION

SPECTRAL clustering is a popular unsupervised classification technique which finds applications in many domains, such as image segmentation [1], text mining [2], and as a general purpose method for data analysis [3], [4], [5]. It relies on the spectrum of a suitably chosen similarity matrix to perform dimensionality reduction before applying a standard clustering algorithm such as K-means. Consider, e.g., the following toy example where $n$ vectors $\mathbf{x}_1, \dots, \mathbf{x}_n \in \mathbb{R}^p$ are separated in two clusters $C^+$, $C^-$ centered around $+\mu, -\mu$ respectively, i.e., $\mathbf{x}_i = \pm\mu + \mathbf{w}_i$ where $\mathbf{w}_i \sim N(0, I_n)$. Then, the dominant eigenvector $\hat{v}$ of the Gram kernel matrix $K = \frac{1}{p}[\mathbf{x}_i^\top \mathbf{x}_j]_{1 \le i,j \le n}$ is an information-theoretically optimal estimator [6] of the vector $\frac{1}{\sqrt{n}}\mathbf{j}$ such that $j_i = \pm 1$ if $\mathbf{x}_i \in C^\pm$. In this case, clustering is achieved with the trivial decision rule $\mathbf{x}_i \to C^\pm$ if $\hat{v}_i \ge 0$.

The achievable performances of spectral clustering can be theoretically predicted thanks to the study of random matrix models corresponding to similarity matrices. For this purpose, random matrix theory offers powerful tools [7], [8], [9]. In particular, it allows to derive the limiting spectral distribution of the kernel matrix and to predict the position of isolated eigenvalues in spiked random matrix models [10], [11], [12]. The latter are of particular importance as, in a wide range of problems, the information of interest can be modeled as a low-rank signal corrupted with noise. In our previous toy example, the data matrix $\mathbf{X} = [\mathbf{x}_1 \ \dots \ \mathbf{x}_n]$ is a rank-one perturbation $\mu \mathbf{j}^\top$ of a noise matrix $\mathbf{W}$ with i.i.d. $N(0, 1)$ entries. In order to theoretically predict the error rate of

spectral clustering for a given signal-to-noise ratio, one must therefore study the behavior of the dominant eigenvectors of the similarity matrix. Tools such as the ones used in [12], [13] allow to express the quality of their alignment with the true underlying signal, i.e., $\frac{1}{\sqrt{n}}|\mathbf{j}^\top \hat{v}|$. Although this tells us when an estimation of the signal is possible (depending on the signal-to-noise ratio) and its efficiency, a precise characterization of the fluctuations of the entries of spiked eigenvectors still lacks to rigorously predict the error rate of spectral clustering. Indeed, in our toy example, the expected error rate $\mathbb{P}(\hat{v}_i j_i < 0)$ cannot be expressed unless the law of $\hat{v}$ is known.

Yet, it is often stated that the entries of $\hat{v}$ have Gaussian fluctuations in the large-dimensional regime, so that $\mathbb{P}(\hat{v}_i j_i < 0)$ is a Gaussian integral. In [14], this result is formally stated but no proof is given. Hence, we fill this missing gap with a rigorous proof of this phenomenon for a general spiked random matrix model. Although we stick to a simple signal + noise model here, the proposed proof is not restricted to Gaussian noise (in fact, the noise only needs to be rotationally invariant) and can easily be adapted to most standard spiked models (such as, notably, the general model considered in [11]). Our result and its proof thus support a wide range of previous works studying the performance of spectral algorithms. The demonstration can be summarized in two simple facts 1) an eigenvector of the kernel matrix can be decomposed into the sum of a deterministic signal part and a random noise part 2) the random part is uniformly distributed on a certain sphere, hence any finite subset of its entries tends to a centered Gaussian vector in the large-dimensional limit.

In this letter, we consider a general signal + noise random matrix model and briefly recall known results regarding its limiting spectral distribution and the behavior of its dominant eigenvalues and eigenvectors. Then, we show that the entries of the kernel eigenvectors indeed have Gaussian fluctuations in the large-dimensional regime. We present a short, self-contained and general proof which is our main contribution. Finally, we illustrate this result with numerical experiments on synthetic and real data.

**Simulations.** Python codes to reproduce simulations are available in the following GitHub repository https://github.com/HugoLebeau/asymptotic_fluctuationsspectral_clustering.

## II. MODEL AND MAIN RESULT

### A. Notations

The symbols $a, \boldsymbol{a}$ and $\mathbf{A}$ respectively denote a scalar, a vector and a matrix. Their entries are $a_i$ and $A_{i,j}$. The set of positive integers below $n$ is $[n] = \{1, \dots, n\}$. The

<footnote>H. Lebeau and R. Couillet are with Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LIG, 38 000 Grenoble, France, (e-mail: {hugo.lebeau, romain.couillet}@univ-grenoble-alpes.fr).</footnote>

<footnote>F. Chatelain is with Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, 38 000 Grenoble, France (e-mail: florent.chatelain@grenoble-inp.fr).</footnote>
cardinality of a set $\mathcal{E}$ is $|\mathcal{E}|$. Given an ordered set of indices $\mathcal{I} = (i_1, \dots, i_{|\mathcal{I}|})$, $[a]_{\mathcal{I}}$ is the vector $[a_{i_1} \dots a_{i_{|\mathcal{I}|}}]^{\top}$. The norm of a vector $\mathbf{x} \in \mathbb{R}^n$ is $\|\mathbf{x}\| = \sqrt{\mathbf{x}^{\top}\mathbf{x}}$. The unit sphere in $\mathbb{R}^n$ is $\mathbb{S}^{n-1} = \{\mathbf{x} \in \mathbb{R}^n \mid \|\mathbf{x}\| = 1\}$. The $n \times n$ identity matrix is $\mathbf{I}_n$. For a real number $x \in \mathbb{R}$, $[x]^+ = \max(0, x)$ and $\delta_x$ is the Dirac measure at $x$. The imaginary unit is denoted i. If the random variable $X$ follows the law $\mathcal{L}$, we write $X \sim \mathcal{L}$. The convergence in distribution of a sequence of random variables $(X_n)_{n \ge 0}$ to $\mathcal{L}$ is denoted $X_n \xrightarrow{n \to +\infty} \mathcal{L}$. Its almost sure convergence to $L$ is denoted $X_n \xrightarrow{n \to +\infty, \text{a.s.}} L$. The multivariate normal distribution with mean $\mu$ and covariance $\Sigma$ is denoted $\mathcal{N}(\mu, \Sigma)$. The set of eigenvalues of a square matrix $\mathbf{A}$ is its spectrum, $\text{Sp} \mathbf{A}$. Given two real-valued sequences $(u_n)_{n \ge 0}$ and $(v_n)_{n \ge 0}$, we write $u_n = \mathcal{O}(v_n)$ if $|u_n/v_n|$ is bounded as $n \to +\infty$ and $u_n \asymp v_n$ if $u_n/v_n \to 1$.

## B. Spiked Matrix Model

Consider the following statistical model

$$ \mathbf{X} = \mathbf{P} + \mathbf{W} \in \mathbb{R}^{p \times n}, \quad \mathbf{P} = \mathbf{L}\mathbf{V}^{\top} \quad (1) $$

with $\mathbf{L} \in \mathbb{R}^{p \times K}$ and $\mathbf{V} = [\mathbf{v}_1 \dots \mathbf{v}_K] \in \mathbb{R}^{n \times K}$ such that $\mathbf{V}^{\top}\mathbf{V} = \mathbf{I}_K$. It models a low-rank signal $\mathbf{P}$ corrupted by additive Gaussian noise $W_{i,j} \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, 1)$. In a spectral clustering perspective, $K$ represents the number of classes and $\mathbf{P} = \mathbf{M}\mathbf{J}^{\top}$ where $\mathbf{M} = [\mu_1 \dots \mu_K]$ is a matrix gathering the $K$ cluster means and $J_{i,k} = 1$ if $\mathbf{x}_i$ is in the $k$-th cluster (i.e., $\mathbf{x}_i = \mu_k + \mathbf{w}_i$) and 0 otherwise. This is congruent with model (1): define the $K \times K$ diagonal matrix $\mathbf{D}$ such that $D_{k,k} = n_k$ where $n_k$ is the number of samples belonging to the $k$-th cluster, then $\mathbf{L} = \mathbf{M}\mathbf{D}^{1/2}$ and $\mathbf{V} = \mathbf{J}\mathbf{D}^{-1/2}$.

Given model (1), we are interested in the reconstruction of $\mathbf{V}$ from the dominant eigenvectors of the Gram kernel matrix $\mathbf{K} = \frac{1}{p}\mathbf{X}^{\top}\mathbf{X}$. We study this problem in the regime where $K$ is fixed and $p, n \to +\infty$ at the same rate, i.e., $0 < c^{\text{def}} \le \lim p/n < +\infty$. This models the fact that, in practice, the number of samples $n$ is comparable to the number of features $p$ and they are both large. Moreover, we make the following assumptions.

**Assumption 1.** All classes are of comparable size, i.e., $\liminf n_k/n > 0$ as $p, n \to +\infty$ for all $k \in [K]$.

**Assumption 2.** $\lim_{p, n \to +\infty} \max_{\substack{1 \le i \le n \\ 1 \le k \le K}} \sqrt{n}V_{i,k}^2 = 0$.

**Assumption 3.** As $n \to +\infty$, the eigenvalues $\ell_1 \ge \dots \ge \ell_K > 0$ of $\frac{1}{n}\mathbf{L}^{\top}\mathbf{L}$ are not degenerate (i.e., have multiplicity one) and the columns of $\mathbf{V}$ are ordered accordingly.

Assumption 3 is only to simplify the presentation of the results so it is not necessary, but often verified in practice. However, Assumption 2 states that $\mathbf{V}$ must be delocalized, i.e., not sparse. It is naturally verified for spectral clustering as a result of Assumption 1 since $\mathbf{V} = \mathbf{J}\mathbf{D}^{-1/2}$, but the results presented below concern the statistical model (1), which is more general and also encompasses PCA for example [13].

<figure><img src="image_1.png" /><figcaption>Fig. 1. Empirical Spectral Distribution (ESD) of $\mathbf{K} = \frac{1}{p}\mathbf{X}^{\top}\mathbf{X}$ and Marčenko-Pastur Distribution (MP). The green dashed lines are the positions $\xi_k$ of isolated eigenvalues predicted by Theorem 1. Experimental setting: $n = 1000$, $p = 2000$, $K = 3$, $(n_1, n_2, n_3) = (333, 334, 333)$, $(\|\mu_1\|, \|\mu_2\|, \|\mu_3\|) = (3, 4, 5)$.</figcaption></figure>

## C. Eigenvalue Distribution and Spiked Eigenvalues

We briefly recall known results on model (1) in order to set the ground for our main result in Theorem 2.

Firstly, the empirical spectral distribution of $\mathbf{K}$, that is $\frac{1}{n} \sum_{\lambda \in \text{Sp} \mathbf{K}} \delta_{\lambda}$, converges weakly almost surely to the Marčenko-Pastur distribution $\mu_{\text{MP}} = [1-c]^{+} \delta_0 + \nu$ where $\nu$ has density supported on $[E_-, E_+]$ with $E_{\pm} = (1 \pm \sqrt{c^{-1}})^2$ [15], [7], [8], [9]. In other words, as $p, n \to +\infty$, the histogram of eigenvalues of $\mathbf{K}$ approaches $\mu_{\text{MP}}$, as depicted in Figure 1.

Due to the low-rank perturbation $\mathbf{P}$, the $K$ dominant eigenvalues of $\mathbf{K}$ may isolate themselves from the bulk characterized by the Marčenko-Pastur distribution if their corresponding signal-to-noise ratios (the eigenvalues of $\frac{1}{n}\mathbf{L}^{\top}\mathbf{L}$) are large enough — they are then called spikes. Their behavior is specified in the following theorem, which is a particular case of Theorem 2 in [13].

**Theorem 1 (Spikes).** Let $(\lambda_k, \hat{v}_k)_{k \in [K]}$ denote the dominant eigenvalue-eigenvector pairs of $\mathbf{K}$ such that $\lambda_1 \ge \dots \ge \lambda_K$. Then, for all $k \in [K]$,

$$ \lambda_k \xrightarrow{\text{a.s.}, p, n \to +\infty} \xi_k^{\text{def}} = \begin{cases} \frac{(\ell_k+c)(\ell_k+1)}{\ell_k c} & \text{if } \ell_k > \sqrt{c} \\ E_+ & \text{otherwise} \end{cases}, \\ |\mathbf{v}_k^{\top} \hat{v}_k|^2 \xrightarrow{\text{a.s.}, p, n \to +\infty} \zeta_k^{\text{def}} = \begin{cases} 1 - \frac{\ell_k+c}{\ell_k(\ell_k+1)} & \text{if } \ell_k > \sqrt{c} \\ 0 & \text{otherwise} \end{cases}. $$

This result states that $\lambda_k$ leaves the bulk if, and only if, $\ell_k > \sqrt{c}$ (this is a well-known phase transition phenomenon [16]) and further gives its almost sure asymptotic position $\xi_k$. This is illustrated in Figure 1. Moreover, Theorem 1 indicates the almost sure asymptotic alignment $\zeta_k$ of the corresponding eigenvector $\hat{v}_k$ with the underlying signal $v_k$. This is depicted in Figure 2 (top row).

It should be noted that the previous results are not restricted to Gaussian noise: up to a control on the moments of the distribution, they can be generalized thanks to an “interpolation trick” [17, Corollary 3.1]. In addition, a similar spectral behavior is observed with non-i.i.d. noise following the realistic assumption that it is concentrated [18], [19].

## D. Fluctuations of Spiked Eigenvectors Entries

The convergence of $|\mathbf{v}_k^{\top} \hat{v}_k|^2$ to $\zeta_k$ stated in Theorem 1 is an important result which justifies the use of the dominant
<figure><img src="image_2.png" /><figcaption>Fig. 2. Dominant eigenvectors of $K = \frac{1}{p}X^\top X$. Top: Coordinates of $\hat{v}_k$ (blue) and the underlying signal $\sqrt{\zeta_k}v_k$ (orange) with $\zeta_k$ given in Theorem 1. The dotted orange lines are the $\pm 1\sigma$-error curves deduced from Theorem 2. Bottom: Histogram of the entries of $\hat{v}_k - \sqrt{\zeta_k}v_k$ (blue) and probability density function of $\mathcal{N}(0, \frac{1-\zeta_k}{n})$ (orange). **Experimental setting:** like in Figure 1.</figcaption></figure>

eigenvectors of $K$ as estimators of the underlying signal $V$. Yet, it is not enough to characterize its reconstruction performance. Indeed, the fluctuations of the entries of $\hat{v}_k$ must be known to fully characterize how it is aligned with $v_k$.

Consider, e.g., the multi-class spectral clustering problem with $P = MJ^\top$. Here, $[v_k]_i = J_{i,k}/\sqrt{n_k}$. Hence, $x_i$ is classified in the $k$-th class if $[\hat{v}_k]_i > [\hat{v}_{k'}]_i$ for all $k' \neq k$. The reconstruction performance thus depends on the probability of correct classification $\mathbb{P}([\hat{v}_k]_i > [\hat{v}_{k'}]_i | J_{i,k} = 1)$. In the theorem below, we show that the entries of $\hat{v}_k$ asymptotically have Gaussian fluctuations around those of $v_k$ with variance $(1 - \zeta_k)/n$, as illustrated in the bottom row of Figure 2.

**Theorem 2.** For all finite ordered set of indices $\mathcal{I} = (i_1, \dots, i_{|\mathcal{I}|}) \subset [n]$ and $k \in [K]$,

$$ \frac{\sqrt{n} [\hat{v}_k - \sqrt{\zeta_k} v_k]_{\mathcal{I}}}{\sqrt{1 - \zeta_k}} \xrightarrow{D} \mathcal{N}(0, I_{|\mathcal{I}|}) $$

with $\hat{v}_k$ such that $v_k^\top \hat{v}_k \ge 0$ (otherwise, consider $-\hat{v}_k$).

This result invokes the quantity $\zeta_k$ introduced in Theorem 1, which quantifies the alignment of $\hat{v}_k$ with $v_k$. Theorem 2 specifies that $[\hat{v}_k]_\mathcal{I}$ behaves like $\mathcal{N}(\sqrt{\zeta_k}[v_k]_\mathcal{I}, \frac{1-\zeta_k}{n}I_{|\mathcal{I}|})$ in the large-dimensional regime. That is, the more $\hat{v}_k$ is aligned with $v_k$ (i.e., the closer $\zeta_k$ is to 1), the more it concentrates around $\sqrt{\zeta_k}v_k$, since the variance is $(1 - \zeta_k)/n$. Furthermore, the entries of $[\hat{v}_k]_\mathcal{I}$ are asymptotically independent for any finite ordered set of indices $\mathcal{I}$. In the multi-class spectral clustering problem considered above, since $\hat{v}_k$ and $\hat{v}_{k'}$ are asymptotically independent if $k' \neq k$ [20, Theorem 4], Theorem 2 yields

$$ \mathbb{P}([\hat{v}_k]_i > [\hat{v}_{k'}]_i | J_{i,k} = 1) \asymp \Phi \left( \sqrt{\frac{n}{n_k} \frac{\zeta_k}{2 - (\zeta_k + \zeta_{k'})}} \right) $$

where $\Phi: x \mapsto \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} e^{-\frac{t^2}{2}} dt$ is the Gaussian cumulative distribution function.

We prove Theorem 2 in Section III below. The proof hinges on the rotational invariance of the noise (Lemma 1). In fact, it does not need the entries of $W$ to be distributed according to the Gaussian law, but only that its distribution be invariant under isometries. This makes it a very general proof, which can easily be adapted to most standard spiked models as those discussed, e.g., in [9, §2.5.4].

### III. PROOF OF MAIN RESULT

Consider the tangent-normal decomposition

$$ \hat{v}_k = \sum_{\kappa=1}^{K} \tau_{\kappa} v_{\kappa} + \sqrt{1 - \|\tau\|^2} \hat{v}_k^{\sharp} \quad (2) $$

where $\hat{v}_k^\sharp = (I_n - VV^\top) \frac{\hat{v}_k}{\sqrt{1 - \|\tau\|^2}}$ is a unit-norm vector orthogonal to the span of $V$ and $\tau = [\tau_1 \dots \tau_K]^\top$ with $\tau_\kappa = v_\kappa^\top \hat{v}_k$ measuring the cosine between $v_\kappa$ and $\hat{v}_k$. Let $O$ be an $n \times n$ orthogonal matrix such that $OV = V$ — i.e., a rotational symmetry about the span of $V$ — and $\tilde{K} \stackrel{\text{def}}{=} OKO^\top$. Then, since $K = \frac{1}{p}X^\top X$ and $X = LV^\top + W$,

$$ \tilde{K} = \frac{1}{p} \left( [OV] L^\top L [OV]^\top + [OV] L^\top [WO^\top] + [WO^\top]^\top L [OV]^\top + [WO^\top]^\top [WO^\top] \right). $$

**Lemma 1.** $W$ and $WO^\top$ are identically distributed.

*Proof.* The distribution of $[WO^\top]_{i,j} = \sum_{k=1}^n W_{i,k}O_{j,k}$ is $\mathcal{N}(0, 1)$ and Cov($[WO^\top]_{i,j}$, $[WO^\top]_{i',j'}$) is 1 if $(i, j) = (i', j')$ and 0 otherwise. Hence $[WO^\top]_{i,j}$ i.i.d. $\mathcal{N}(0, 1)$. $\square$

According to the previous lemma, $\tilde{K}$ follows the same model as $K$ since $OV = V$. Therefore, its $k$-th dominant eigenvector can likewise be decomposed as

$$ \tilde{v}_k = \sum_{\kappa=1}^{K} \tilde{\tau}_{\kappa} v_{\kappa} + \sqrt{1 - \|\tilde{\tau}\|^2} \tilde{v}_k^{\sharp} $$

with $\tilde{\tau}_\kappa = v_\kappa^\top \tilde{v}_k$ and $\tilde{v}_k^\sharp = (I_n - VV^\top) \frac{\tilde{v}_k}{\sqrt{1 - \|\tilde{\tau}\|^2}}$ identically distributed to $\hat{v}_k^\sharp$. Yet, $\tilde{v}_k = O\hat{v}_k$. Thus, $\hat{v}_k^\sharp$ and $O\hat{v}_k^\sharp$ are identically distributed for all $n \times n$ orthogonal matrix $O$ such that $OV = V$. Consequently, denoting $\eta$ the probability distribution of $\hat{v}_k^\sharp$ and $V^\perp = \{w \in \mathbb{R}^n | V^\top w = 0\}$, then, for all $x, y \in \mathbb{S}^{n-1} \cap V^\perp$, we have $d\eta(x) = d\eta(Ox) = d\eta(y)$ with $O = I_n - \frac{(x-y)(x-y)^\top}{1-x^\top y}$ satisfying $OV = V$. This shows that $\hat{v}_k^\sharp$ is uniformly distributed on $\mathbb{S}^{n-1} \cap V^\perp$.

Then, $\hat{v}_k^\sharp$ can be written as $Uu$ where $u$ is uniformly distributed on $\mathbb{S}^{n-1-K} \subset \mathbb{R}^{n-K}$ and $U \in \mathbb{R}^{n \times (n-K)}$ is such that $U^\top U = I_{n-K}$ and $U^\top V = 0$ (the columns of $U$ form
an orthonormal basis of $V^\perp$ in $\mathbb{R}^n$). We use the following theorem to identify the asymptotic distribution of $\sqrt{n}[\hat{v}_k^\sharp]_\mathcal{I}$.

**Theorem 3** ([21], [22]). The characteristic function of a vector $w$ uniformly distributed on $\mathbb{S}^{n-1}$ is given by $\varphi_w(t) \stackrel{\text{def}}{=} \mathbb{E}[e^{i t^\top w}] = \Omega_n(\|t\|)$ where $\Omega_n$ is such that $r \mapsto \Omega_n(r\sqrt{n})$ converges uniformly in $r \ge 0$ to $r \mapsto e^{-r^2/2}$ as $n \to +\infty$.

Let $t \in \mathbb{R}^n$ be such that $t_i = 0$ if $i \notin \mathcal{I}$. The characteristic function of $\sqrt{n}[\hat{v}_k^\sharp]_\mathcal{I}$ is

$$ \varphi_{\sqrt{n}[\hat{v}_k^\sharp]_\mathcal{I}}(t_{i_1}, \dots, t_{i_{|\mathcal{I}|}}) = \mathbb{E} \left[ e^{i \sqrt{n} t^\top U u} \right] \\ = \Omega_{n-K} (\sqrt{n} \|U^\top t\|) $$

and $\|U^\top t\| = \sqrt{\|t\|^2 - \|V^\top t\|^2} = \|t\| + \mathcal{O}(\|V^\top t\|^2)$. According to Assumption 2, $\sqrt{n}\|V^\top t\|^2 \to 0$ as $p, n \to +\infty$, thus $\Omega_{n-K}(\sqrt{n}\|U^\top t\|) = \Omega_{n-K}(\sqrt{n-K}\|t\| + \epsilon_n)$ with $\epsilon_n \to 0$ as $p, n \to +\infty$ and

$$ \left| \Omega_{n-K} (\sqrt{n-K} \|t\| + \epsilon_n) - e^{-\frac{1}{2} \|t\|^2} \right| \le \left| \Omega_{n-K} (\sqrt{n-K} \|t\| + \epsilon_n) - e^{-\frac{1}{2} \left[ \|t\| + \epsilon_n / \sqrt{n-K} \right]^2} \right| + \left| e^{-\frac{1}{2} \left[ \|t\| + \epsilon_n / \sqrt{n-K} \right]^2} - e^{-\frac{1}{2} \|t\|^2} \right| $$

As $p, n \to +\infty$, the first term vanishes from the uniform convergence given in Theorem 3 and the second term vanishes by continuity. Therefore, $\varphi_{\sqrt{n}[\hat{v}_k^\sharp]_\mathcal{I}}(t_{i_1}, \dots, t_{i_{|\mathcal{I}|}}) \to e^{-\frac{\|t\|^2}{2}}$ and, by Lévy's continuity theorem [23], we can conclude that $\sqrt{n}[\hat{v}_k^\sharp]_\mathcal{I} \xrightarrow{D} \mathcal{N}(0, I_{|\mathcal{I}|})$ as $p, n \to +\infty$. And, finally, given decomposition (2), Theorem 1 and the independence of $\hat{v}_k$ and $\hat{v}_{k'}$ if $k' \ne k$ [20, Theorem 4],

$$ \hat{v}_k = \sqrt{\zeta_k} v_k + \sqrt{1 - \zeta_k} \hat{v}_k^\sharp + \varepsilon \quad \text{with} \quad \|\varepsilon\| \xrightarrow{\substack{\text{a.s.} \\ p, n \to +\infty}} 0. $$

This concludes the proof of Theorem 2.

## IV. NUMERICAL EXPERIMENTS

To illustrate this result, we conduct a first experiment on synthetic data following model (1) with $K = 3$ classes of equal size and $(\|\mu_1\|, \|\mu_2\|, \|\mu_3\|) = (3, 4, 5)$. The $x_i$'s are ordered by class. The spectral distribution of $K$ is plotted in Figure 1 and Figure 2 shows the dominant eigenvectors with the histograms of residuals $\hat{v}_k - \sqrt{\zeta_k} v_k$. We observe a very good fit of the latter to the probability density function of $\mathcal{N}(0, \frac{1-\zeta_k}{n})$ — the $\hat{v}_k$'s exactly correspond to a deterministic signal $\sqrt{\zeta_k} v_k$ corrupted by additive centered Gaussian noise. The signal + noise structure of model (1) has been transferred to the spectral estimator of $V$.

Then, we conduct a second experiment on the Fashion-MNIST dataset [24] consisting of 28 × 28 images of clothes separated in 10 classes of size 7000 each. We select two classes $k_1, k_2$ and perform binary spectral clustering using the dominant eigenvector of $K = \frac{1}{p}X^\top X$ where the columns of $X$ are the 784 pixels of the images from classes $k_1$ and $k_2$. The dimension of $X$ is thus 784 × 14000. Here, we assume a similar model as our toy example in the introduction: $X = \mu j^\top + W$ where $j_i = \pm 1$ depending on the class of the $i$-th image. Thus, according to Theorem 2,

<figure><img src="image_3.png" /><figcaption>Fig. 3. Observed (upper right, blue) and predicted (lower left, orange) classification accuracies of binary spectral clustering on the Fashion-MNIST dataset [24].</figcaption></figure>

the $i$-th entry of the dominant eigenvector $\hat{v}$ asymptotically follows $\mathcal{N}(\sqrt{\zeta} \frac{j_i}{\sqrt{n}}, \frac{1-\zeta}{n})$ and, given $j$, $\zeta$ can be estimated as $(\sum_{i=1}^n \frac{j_i}{\sqrt{n}} \hat{v}_i)^2$. We can then compare the observed accuracy $\frac{1}{n} \sum_{i=1}^n \hat{v}_{i,j_i>0}$ to the one expected from Theorem 2, $\mathbb{P}(\hat{v}_{i,j_i>0}) \approx \Phi(\sqrt{\frac{\zeta}{1-\zeta}})$. The results are presented in Figure 3 for each pair of classes $k_1, k_2$.

We find a very good agreement between the observed and predicted accuracies, regardless of whether the problem is easy (e.g., Trouser vs Sandal) or hard (e.g., Bag vs Ankle Boot). This observation confirms the general scope of Theorem 2: starting from real data $X$ which is clearly not Gaussian, the normal distribution naturally emerges in the fluctuations of the entries of the large-dimensional eigenvector $\hat{v}$.

## V. CONCLUSION

After recalling known results on spectral clustering under a general signal + noise random matrix model, we have shown that the entries of spiked eigenvectors have Gaussian fluctuations in the large-dimensional regime. This formalizes and clearly states a result which is often implicitly assumed in many problems, without ever being actually proven. The proposed proof relies solely on the rotational invariance of the noise. It is thus very general and can easily be extended to most standard spike models. Numerical experiments have demonstrated the universality of this phenomenon: the Gaussian behavior of the entries of spike eigenvectors can even be observed on real unprocessed data. This allows to accurately predict the classification performance of spectral clustering. An interesting problem for future work is to understand how these results extend to more exotic spike models such as [25].