Asymptotic properties of spiked eigenvalues and eigenvectors of signal-plus-noise matrices with their applications

Consider a signal-plus-noise model with the form of

$\displaystyle\mathbf{X}_{n}=\mathbf{A}_{n}+\bm{\Sigma}^{1/2}\mathbf{W}_{n}\in\mathbb{R}^{p\times n},$

(1)

where $\mathbf{A}_{n}$ is the signal matrix with finite rank, $\mathbf{W}_{n}$ consists of i.i.d. random variables, and $\bm{\Sigma}$ accounts for the covariance structure in the noise. Such a model is popular in many fields including machine learning (Yang et al.,, 2016), matrix denoising (Nadakuditi,, 2014) or signal processing (Vallet et al.,, 2012). When $\bm{\Sigma}$ is an identity matrix, there has been a huge amount of work on eigenvalues and eigenvectors for such signal-plus-noise type matrices. To name a few, Loubaton and Vallet, (2011) derived the almost sure limits of eigenvalues, Ding, (2020) obtained the limits and convergent rates of the leading eigenvalues and eigenvectors, Bao et al., (2021) showed the distributions of the principal singular vectors and singular subspaces. When $\bm{\Sigma}_{n}$ is set to be a diagonal matrix, Hachem et al., (2013) investigated the limiting behavior of the random bilinear form of the sample covariance matrix under a separable model, which includes the case of $\bm{\Sigma}$ being diagonal in (1). When the signal-to-noise ratio tends to infinity, i.e., the ratio of the spectral norm of the signal part to the noise part tends to infinity, Cape et al., (2019) also considered the asymptotic properties of spiked eigenvectors under Model (1). By imposing Gaussianity on $\mathbf{W}_{n}$, Han et al., (2021) provided a eigen-selected spectral clustering method with theoretical justifications.

However, the assumptions that $\bm{\Sigma}$ is an identity or diagonal matrix, and the signal-to-noise ratio tends to infinity, seem to be restricted and hard to verify in practice. In this paper, we aim to investigate the asymptotic properties of the eigenvalues of $\mathbf{X}_{n}\mathbf{X}_{n}^{\top}$, as well as both the left and right spiked singular vectors of $\mathbf{X}_{n}$ under the regime where $p/n\rightarrow c>0$, with mild regularity conditions towards $\bm{\Sigma}$ and $\mathbf{A}_{n}$, and mild moment assumptions on $\mathbf{W}_{n}$. To the best of our knowledge, we first systematically study the properties of eigenvalues and eigenvectors of Model (1) under such mild conditions. Specifically, we consider

$$\mathbf{S}_{n}:=\mathbf{X}_{n}\mathbf{X}_{n}^{\top}=(\mathbf{A}_{n}+\bm{\Sigma}^{1/2}\mathbf{W}_{n})(\mathbf{A}_{n}+\bm{\Sigma}^{1/2}\mathbf{W}_{n})^{\top}$$

(2)

and

$\displaystyle\tilde{\mathbf{S}}_{n}:=\mathbf{X}_{n}^{\top}\mathbf{X}_{n}=(\mathbf{A}_{n}+\bm{\Sigma}^{1/2}\mathbf{W}_{n})^{\top}(\mathbf{A}_{n}+\bm{\Sigma}^{1/2}\mathbf{W}_{n}).$

(3)

In order to obtain the asymptotic properties of spiked eigenvectors of $\mathbf{S}_{n}$ and $\tilde{\mathbf{S}}_{n}$, we analyze the quadratic forms involving the resolvents $Q_{n}(z)$ and $\tilde{Q}_{n}(z)$ of matrices $\mathbf{S}_{n}$ and $\tilde{\mathbf{S}}_{n}$ defined as

$\displaystyle Q_{n}(z)=(\mathbf{S}_{n}-z\mathbf{I})^{-1}$

(4)

and

$\displaystyle\tilde{Q}_{n}(z)=(\tilde{\mathbf{S}}_{n}-z\mathbf{I})^{-1},$

(5)

respectively, where $z\in\mathcal{C}^{+}$ and $\mathbf{I}$ refer to an identity matrix with comparable sizes. The study on the spiked eigenvalues leverages the main results developed in Liu et al., (2022).

To demonstrate the use of the theoretical results, we consider applications in spectral clustering. When each column of $\mathbf{A}_{n}$ can be only chosen from a finite number of distinct unknown deterministic vectors, (1) can be regarded as a collection of samples generated from a mixture model. Thus, in a vector form, the $i$-th column of Model (1) can be written as

$$\mathbf{x}_{i}=\mathbf{a}_{i}+\bm{\Sigma}^{1/2}\mathbf{w}_{i},\in\mathbb{R}^{p}$$

(6)

where $\mathbf{a}_{i}=\bm{\mu}_{s}/\sqrt{n}$ for some $s\in\{1,\ldots,K\}$ if $i\in\mathcal{V}_{s}\subseteq\{1\ldots,n\}$. The normalized constant $\sqrt{n}$ in $\mathbf{a}_{i}$ is to unify the Assumption 1 below. Here $\cup_{s=1}^{k}\mathcal{V}_{s}=\{1\ldots,n\}$ and $\mathcal{V}_{s}\cap\mathcal{V}_{t}=\emptyset$ for any $s\neq t$, and $K$ actually refers to the number of the different distributions (i.e., clusters) in a mixture model. One should also note that the labels are unknown in clustering problems. Numerous literatures investigate mixture models. In statistics, Redner and Walker, (1984) considered the clustering problem for the Gaussian mixture model in low dimensional cases, while Cai et al., (2019) considered the high dimensional cases. Some classical techniques about clustering were also proposed in past decades; see e.g., MacQueen et al., (1967), Bradley et al., (1999), Kaufman and Rousseeuw, (1987), Maimon and Rokach, (2005) and Duda and Hart, (1973). In empirical economics, mixture models are used to introduce unobserved heterogeneity. An important example of this setup from the econometrics literature is Keane and Wolpin, (1997), which investigated the clustering problem in labor markets. Such models also arise in analyzing some classes of games with multiple Nash equilibria. See for example, Berry and Tamer, (2006), Chen et al., (2014) and others.

Our main theoretical contribution is to precisely characterize the first-order limits of the eigenvalues and eigenvectors of $\mathbf{S}_{n}$ and $\widetilde{\bf S}_{n}$. There are two observations that can be obtained based on our main theoretical results that are somewhat surprising, as they exhibit some overlaps with findings in the literature, albeit in different scopes of problems. The first is that the limits of the spiked eigenvalues of $\mathbf{S}_{n}$ coincide with these of the sample covariance matrices without the signal part, by letting the population covariance matrix be $\mathbf{A}_{n}\mathbf{A}_{n}^{\top}+\bm{\Sigma}$, see the discussion below Theorem 2. The second is that, the spiked right singular vectors of $\mathbf{X}_{n}$ have an intrinsic block structure if $\mathbf{A}_{n}$ contains a finite number of distinct deterministic factors, even for a moderate signal-to-noise ratio. Our Corollary 1 precisely quantifies the deviation of the right singular vector from a vector with entries having a group structure. This finding is highly relevant to the field of spectral clustering, which has been extensively discussed in the literature. It is worth noting that many existing studies assume strong moment conditions on the noise and consider scenarios where the signal-to-noise ratio tends to infinity.

As applications, we propose a method to estimate the number of clusters by leveraging the asymptotic limits of sample eigenvalues. We also discuss how the theoretical results intuitively explain why the spiked eigenvectors have clustering power in the context of spectral clustering.

The remaining sections are organized as follows. In Section 2 we state the main results of the quantities of (4) and (5) under Model 1. Section 3 includes the applications in terms of clustering and classification. In Section 4, we demonstrate some numerical results regarding the applications mentioned in Section 3. All the proofs are relegated to the Appendix part and the supplementary materials.

Conventions: We use $C$ to denote generic large constant independent of $n$, and its value may change from line to line. $a\wedge b=\min\{a,b\}$, $\mathbf{1}$ and $\mathbf{I}$ refer to a vector with all entries being one and an identity matrix with a comparable size, respectively. We let $\|\cdot\|$ denote the Euclidean norm of a vector or the spectral norm of a matrix.

In this section, we mainly investigate the limits of the eigenvalues and eigenvectors of $\mathbf{S}_{n}$ and $\widetilde{\bf S}_{n}$ defined in (2) and (3), respectively. We first impose some mild conditions on $\mathbf{W}_{n}$ for establishing the asymptotic limits of the eigenvalues and eigenvectors:

We consider the high dimensional setting specified by the following assumption.

Note that when $\mathbf{A}_{n}$ has a bounded rank, the limiting spectral distribution of $\mathbf{X}_{n}\mathbf{X}_{n}^{\top}$ is the same as that of the model by setting $\mathbf{A}_{n}$ as a zero matrix. This can be directly concluded by the rank inequality, see Theorem A.43 of Bai and Silverstein, (2010). However, to investigate the limiting behaviors of the spiked eigenvalues and eigenvectors under Model (1), more assumptions on $\mathbf{A}$ and $\bm{\Sigma}$ are required.

The key technical tool is the deterministic equivalents of $Q_{n}$ and $\tilde{Q}_{n}$ in (4) and (5). We introduce it first as it requires weaker assumptions than the main results on spiked eigenvectors and eigenvalues, and may be of independent interest. For any $z\in\mathcal{C}^{+}$, let $\tilde{r}_{n}(z)\in\mathcal{C}^{+}$ be the unique solution to the equation

$\displaystyle z=-\frac{1}{\tilde{r}_{n}}+c_{n}\int\frac{tdF^{\mathbf{R}_{n}}(t)}{1+t\tilde{r}_{n}},$

(8)

where $F^{\mathbf{R}_{n}}(t)$ is the empirical spectral distribution of $\mathbf{R}_{n}$. Proposition 1 below provides the deterministic equivalence of $Q_{n}$ and $\tilde{Q}_{n}$.

Based on Proposition 1, we now first focus on the eigenvectors corresponding to the spiked eigenvalues of $\mathbf{R}_{n}$. The following assumption is needed.

Let $\hat{\mathbf{v}}_{k}\in\mathbb{R}^{p}$ and $\hat{\mathbf{u}}_{k}\in\mathbb{R}^{n}$ be the eigenvector associated with the $k$-th largest (spiked) eigenvalue of $\mathbf{S}_{n}$ and $\tilde{\mathbf{S}}_{n}$, respectively. The following theorem characterizes the asymptotic behaviours of $\hat{\mathbf{v}}_{k}$ and $\hat{\mathbf{u}}_{k}$.

The asymptotic behaviour of the spiked eigenvalues is also considered, and thus some more notations are also required. Similar to Bai and Yao, (2012), for the spiked eigenvalue $\gamma$ outside the support of $H$ and $\gamma\neq 0$, we define

$$\varphi(\gamma)=z(-\frac{1}{\gamma})=\gamma\left(1+c\int\frac{t\mathrm{d}H(t)}{\gamma-t}\right),$$

(18)

where $z$ is regarded as the function defined in (8) with its domain extended to the real line. As defined in Bai and Yao, (2012), a spiked eigenvalue $\gamma$ is called a distant spike if $\varphi^{\prime}(\gamma)>0$ which is coincident to Assumption 4, and a close spike if $\varphi^{\prime}(\gamma)\leq 0$. Note that $\mathbf{S}_{n}$ and $\tilde{\mathbf{S}}_{n}$ share the same nonzero eigenvalues, and we denote by $\lambda_{1}\geq\ldots\geq\lambda_{p\wedge n}>0.$

We verify that our main results match with the corresponding parts in Ding, (2020). As the first $K$ eigenvectors of $\mathbf{R}_{n}=\mathbf{A}_{n}\mathbf{A}_{n}^{\top}+\mathbf{I}_{n}$ are the same as the left singular vectors of $\mathbf{A}_{n}$, the singular value decomposition of $\mathbf{A}_{n}$ can be written as $\mathbf{A}_{n}=\sum_{i=1}^{K}d_{i}\xi_{i}\zeta_{i}^{\top}$ where $\xi_{i}$ is the eigenvector associated with $\mathbf{R}_{n}$. Then $\gamma_{k}=d_{k}^{2}+1$ for $k=1,\cdots,K$, and $\gamma_{k}=1$ for $(K+1)\leq k\leq p$. Theorem 2 implies that

$$\lambda_{k}\overset{P}{\rightarrow}\varphi(\gamma_{k})=(d_{i}^{2}+1)(1+d_{i}^{-1}c).$$

By taking $\mathbf{u}_{n}=\zeta_{k}$ in (16), we find $\mathbf{u}_{n}^{\top}\mathbf{A}_{n}^{\top}\xi_{k}=d_{k}$, and $\eta_{k}=1-c_{n}d_{k}^{-4}+O(n^{-1})$, thus

$$\zeta_{k}^{\top}\hat{\mathbf{u}}_{k}\hat{\mathbf{u}}_{k}^{\top}\zeta_{k}-(d_{k}^{4}-c_{n})/[d_{k}^{2}(1+d_{k}^{2})]=O_{P}(n^{-1/2}).$$

These limits coincide with $p(d_{k})$ and $a_{2}(d_{k})$ defined in (2.6) and (2.9) of Ding, (2020), respectively.

One may wonder whether the asymptotic distributions of the spiked eigenvalues and eigenvectors of $\mathbf{X}_{n}\mathbf{X}_{n}^{\top}$ as those of $\mathbf{R}_{n}^{1/2}\mathbf{W}_{n}\mathbf{W}_{n}\mathbf{R}_{n}^{1/2}$, given that their first order limits coincide. Several recent studies investigated the latter model, including (Jiang and Bai,, 2021; Zhang et al.,, 2022; Bao et al.,, 2022), etc. Through simulations, we observe different asymptotic variances between the two models, as indicated by Table 5.

The aforementioned theoretical results are all built on $\mathbf{S}_{n}$ or $\tilde{\mathbf{S}}_{n}$ that refer to the noncentral covariance matrices. In some situations, the centered versions are also of interest. Specifically, we consider the corresponding covariance matrices

$$\bar{\mathbf{S}}_{n}=(\mathbf{X}_{n}-\bar{\mathbf{X}}_{n})(\mathbf{X}_{n}-\bar{\mathbf{X}}_{n})^{\top},$$

and

$$\tilde{\bar{\mathbf{S}}}_{n}=(\mathbf{X}_{n}-\bar{\mathbf{X}}_{n})^{\top}(\mathbf{X}_{n}-\bar{\mathbf{X}}_{n}),$$

where $\bar{\mathbf{X}}_{n}=\bar{\mathbf{x}}_{n}\mathbf{1}^{\top}$ and $\bar{\mathbf{x}}_{n}=\sum_{k=1}^{n}\mathbf{x}_{k}/n$. Let $\Phi=\mathbf{I}-\mathbf{1}\mathbf{1}^{\top}/n$ and denote the spectral decomposition of $\bar{\mathbf{R}}_{n}=\mathbf{A}_{n}\Phi\mathbf{A}_{n}^{\top}+\bm{\Sigma}$ by $\bar{\mathbf{R}}_{n}=\sum_{k=1}^{p}\bar{\gamma}_{k}\bar{\xi}_{k}\bar{\xi}_{k}^{\top}$, where $\bar{\gamma}_{1}>\ldots>\bar{\gamma}_{\bar{K}}>\ldots\geq\bar{\gamma}_{p}$. Here $\bar{K}$ may be equal to $K$ or $K-1$ in some different cases. Moreover, define the corresponding resolvent $\bar{Q}_{n}(z)$ and $\bar{\tilde{Q}}_{n}(z)$ of matrix $\bar{\mathbf{S}}_{n}$ and $\bar{\tilde{\mathbf{S}}}_{n}$, respectively:

$\displaystyle\bar{Q}_{n}(z)=(\bar{\mathbf{S}}_{n}-z\mathbf{I})^{-1},\quad\bar{\tilde{Q}}_{n}(z)=(\bar{\tilde{\mathbf{S}}}_{n}-z\mathbf{I})^{-1}.$

Based on the given notations, we established the corresponding results for the centralized sample covariance matrices:

Relying on Proposition 2, we also have the following conclusion for the spiked eigenvalues and the corresponding eigenvectors of $\bar{\mathbf{S}}_{n}$ and $\bar{\tilde{\mathbf{S}}}_{n}$.

In this section, based on the results in Section 2, we aim to develop some potential applications. Spectral clustering has been used in practice frequently in data science and the theoretical underpinning of such a method has received extensive interest in recent years; see e.g., (Couillet et al.,, 2016; Zhou and Amini,, 2019; Löffler et al.,, 2021), etc. This section is to have a deep insight into the spectral clustering based on the Model (6). Moreover, we also propose a new criterion to estimate the number of clusters. Recalling (6), for any $i\in\mathcal{V}_{s}$, there is $\mbox{E}\mathbf{x}_{i}=\bm{\mu}_{s}/\sqrt{n}$, where $s=1,\ldots,K$. Let $\mathbf{N}=[\bm{\mu}_{1},\ldots,\bm{\mu}_{K}]/\sqrt{n}\in\mathbb{R}^{p\times K}$, $\mathbf{H}=[\mathbf{h}_{1},\ldots,\mathbf{h}_{K}]\in\mathbb{R}^{n\times K}$, $\mathbf{h}_{s}=(\mathbf{h}_{s}(1),\ldots,\mathbf{h}_{s}(n))^{\top}\in\mathbb{R}^{n}$, where $\mathbf{h}_{s}(i)=1$ if $i\in\mathcal{V}_{s}$ and $\mathbf{h}_{s}(i)=0$ otherwise. In a matrix form, write

$$\mathbf{X}_{n}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]=\mathbf{N}\mathbf{H}^{\top}+\bm{\Sigma}^{1/2}\mathbf{W}_{n}$$

Notice that

$$\mbox{E}(\tilde{\mathbf{S}}_{n})=\mathbf{H}\mathbf{N}^{\top}\mathbf{N}\mathbf{H}^{\top}+\frac{\mbox{tr}\bm{\Sigma}}{n}\mathbf{I}_{n}.$$

The block structure of $\mbox{E}(\tilde{\mathbf{S}}_{n})$ (except the diagonal positions) is similar to that of stochastic block models (SBM). This motivates one to use spectral clustering for high dimensional data with different means across groups.

To do the clustering, it is of interest to estimate the number of clusters, i.e., estimation of $K$. There exist plenty of approaches to estimate the number of the clusters. To name a few, Thorndike, (1953) proposed the Elbow method that aims to minimize the within-group sum of squares (WSS); Silhouette index (Rousseeuw,, 1987) is a measure of how similar an object is to its own cluster compared to other clusters, which takes values in $[-1,1]$; Tibshirani et al., (2001) proposed a gap statistic to estimate the number of clusters, etc. These methods either lack theoretical guarantees or have some restrictions in computation or settings. Hence, here we propose a theoretical guarantee and easily implemented approach to estimate the number of clusters. Notice that under Model (6), the number of the spiked eigenvalues of $\mathbf{S}_{n}$ or $\tilde{\mathbf{S}}_{n}$ is the same as the number of clusters if the means in terms of the different clusters are not linearly correlated. The estimation of the number of spikes in different models has been discussed in multiple literatures, and mostly are based on the setting of $\bm{\Sigma}=\mathbf{I}$; see e.g., Bai et al., (2018).

Motivated by the work of Bai et al., (2018) and Theorem 2, we propose two criteria to estimate the number of clusters. Without loss of generality, we assume $0<c<1$. Let

$$\begin{split}&\text{EDA}_{k}=-n(\lambda_{1}-\lambda_{k+1})+n(p-k-1)\log\widetilde{\theta}_{p,k}+2pk,\\ &\text{EDB}_{k}=-n\log(p)\cdot(\lambda_{1}-\lambda_{k+1})+n(p-k-1)\log\widetilde{\theta}_{p,k}+(\log n)pk,\end{split}$$

(21)

where $\widetilde{\theta}_{p,k}=\frac{1}{p-k-1}\sum_{i=k+1}^{p-1}\theta_{i}^{2}$, and $\theta_{k}=\exp\{\lambda_{n,k}-\lambda_{n,k+1}\},k=1,2,\ldots,p-1$.

We estimate the number of clusters by

	$\displaystyle\hat{K}_{\text{EDA}}$	$\displaystyle=$	$\displaystyle\arg\min\limits_{k=1,\ldots,w}\frac{1}{n}\text{EDA}_{k},$		(22)
	$\displaystyle\hat{K}_{\text{EDB}}$	$\displaystyle=$	$\displaystyle\arg\min\limits_{k=1,\ldots,w}\frac{1}{n}\text{EDB}_{k},$		(23)

where $w$ is the prespecified number of clusters satisfying $w=o(p)$. Note that under conditions of Theorem 2, it follows that for $k=1,2,\ldots,K-1$,

$$\theta_{k}\overset{p}{\rightarrow}\exp\{\varphi(\gamma_{k})-\varphi(\gamma_{k+1})\},~{}~{}\theta_{K}\overset{p}{\rightarrow}\exp\{\varphi(\gamma_{K})-\mu_{1}\},$$

(24)

where function $\varphi$ and $\mu_{1}$ are defined in (18) and (20), respectively. For simplicity, denote the limit of $\theta_{k}$ by $\xi_{k}$ for $k=1,\ldots,K$. Define two sequences $\{a_{s}\}_{s=2}^{K}$ and $\{b_{s}\}_{s=2}^{K}$ as follows

$$\begin{split}&a_{s}=\xi_{s}^{2}+\log\xi_{s}-2c-1+a_{s+1}\text{ and }a_{K+1}=0\text{ for }s=2,\ldots,K,\\ &b_{s}=\xi_{s}^{2}+{\log p}\log\xi_{s}-c\log n-1+b_{s+1}\text{ and }b_{K+1}=0\text{ for }s=2,\ldots,K.\end{split}$$

(25)

We propose two gap conditions for EDA and EDB, respectively, i.e.,

	$\displaystyle\min\limits_{s=2,\ldots,K}a_{s}$	$\displaystyle>$	$\displaystyle 0,$		(26)
	$\displaystyle\min\limits_{s=2,\ldots,K}b_{s}$	$\displaystyle>$	$\displaystyle 0.$		(27)

Note that Theorem 2 and (24) are obtained when the leading eigenvalues are bounded. Here we also investigate the cases of the leading eigenvalues tending to infinity.

Based on Theorem 2 and Lemma 1, we derive the consistency of $\hat{K}_{\text{EDA}}$ as follows.

In Bai et al., (2018), BIC is consistent when $\lambda_{K}\rightarrow\infty$ at a rate faster than $\log n$, which makes BIC less capable of detecting signals. This is because BIC has a more strict penalty coefficient $\log n$ compared to the penalty coefficient “2” in AIC. For the EDB construction of selecting the number of clusters, we add the coefficient “$\log p$” to the first term so that the spikes do not need to be very large and only the corresponding gap condition for EDB is required. By the analogous proof strategy of Theorem 4, we obtain the consistency of EDB as follows.

Once the estimator of the number of clusters is available, we can conduct spectral clustering. Specifically, let the eigenvectors corresponding to the first $\hat{K}$ eigenvalues of $\tilde{\mathbf{S}}_{n}$ be $\widehat{\mathbf{U}}=[\hat{\mathbf{u}}_{1},\ldots,\hat{\mathbf{u}}_{\hat{K}}]\in\mathbb{R}^{n\times\hat{K}}$. We then apply the following K-means optimization to the $\widehat{\mathbf{U}}$, i.e.,

$\displaystyle\mathbf{U}^{*}=\arg\max_{U\in\mathcal{M}_{n,\hat{K}}}\|U-\widehat{\mathbf{U}}\|_{F}^{2},$

(28)

where $\mathcal{M}_{n,K}=\{U\in\mathbb{R}^{n\times K}:U\text{ has at most $K$ distinct rows}\}$. Then, we return $\hat{\mathcal{V}}_{1},\ldots,\hat{\mathcal{V}}_{\hat{K}}$ as the indices for each cluster. From (28), we see that the spectral clustering is conducted from the obtained $\widehat{\mathbf{U}}$, and hence we look into the properties of $\widehat{\mathbf{U}}$.

In this section, we first evaluate the performance of the proposed criteria in the estimation of the number of clusters discussed in Section 3. Denote the sets of under-estimated, exactly estimated and over-estimated models by $\mathcal{F}_{-},\mathcal{F}_{*}$ and $\mathcal{F}_{+}$, respectively, i.e.,

$$\mathcal{F}_{-}=\{1,\ldots,K-1\},~{}~{}\mathcal{F}_{*}=\{K\},~{}~{}\mathcal{F}_{+}=\{K+1,\ldots,w\}.$$

The selection percentages corresponding to $\mathcal{F}_{-},\mathcal{F}_{*}$ and $\mathcal{F}_{+}$ are computed by $1000$ repetitions. Suppose that the entries of $\mathbf{W}_{n}$ are i.i.d. with the following distributions:

• •

  Standard normal distribution: $w_{i,j}\sim\mathcal{N}(0,1)$.

• •

  Standardized $t$ distribution with 8 degrees of freedom: $w_{i,j}\sim t_{8}/\sqrt{\text{Var}(t_{8})}$.

• •

  Standardized Bernoulli distribution with probability $1/2$: $w_{i,j}\sim(\text{Bernoulli}(1,1/2)-1/2)/(1/2)$.

• •

  Standardized chi-square distribution with 3 degrees of freedom: $w_{i,j}\sim(\chi^{2}(3)-3)/\sqrt{\text{Var}(\chi^{2}(3))}=(\chi^{2}(3)-3)/\sqrt{6}$

For comparison, three different methods are also considered: Average Silhouette Index (Rousseeuw, (1987)), Gap Statistic (Tibshirani et al., (2001)) and BIC with degrees of freedom (David, (2020)), denoted by ASI, GS and BICdf, respectively. This section considers the situations with $0<c<1$, and the cases with $c>1$ are demonstrated in the supplementary material. Here we set $c=1/3,3/4$ and the largest number of possible clusters $w=\lfloor 6\cdot n^{0.1}\rfloor$. Different means in terms of different clusters and the covariance matrices are set as follows :

Case 1. Let $\bm{\mu}_{1}=(5,0,-4,0,0,\ldots,0)^{\top}\in\mathbb{R}^{p}$, $\bm{\mu}_{2}=(0,4,0,-6,0,\ldots,0)^{\top}\in\mathbb{R}^{p}$, $\bm{\mu}_{3}=(0,-5,-5,0,0,\ldots,0)^{\top}\in\mathbb{R}^{p}$, $\bm{\mu}_{4}=(-6,0,0,6,0,\ldots,0)^{\top}\in\mathbb{R}^{p}$, and $\bm{\Sigma}=(\sigma_{i,j})_{p\times p}$, where $\sigma_{i,j}=0.2^{|i-j|}$. Define

$$\mathbf{A}_{n}=\big{(}\underbrace{\bm{\mu}_{1},\ldots,\bm{\mu}_{1}}_{n_{1}},\underbrace{\bm{\mu}_{2},\ldots,\bm{\mu}_{2}}_{n_{2}},\underbrace{\bm{\mu}_{3},\ldots,\bm{\mu}_{3}}_{n_{3}},\underbrace{\bm{\mu}_{4},\ldots,\bm{\mu}_{4}}_{n_{4}}\big{)},$$

where $n_{1}=n_{3}=0.3n,~{}n_{2}=n_{4}=0.2n$. Therefore, the true number of clusters is $K=4$.

Case 2. Let $\bm{\mu}_{1}=(3,0,0,0,\ldots,0)^{\top}\in\mathbb{R}^{p}$, $\bm{\mu}_{2}=(0,3,0,0,\ldots,0)^{\top}\in\mathbb{R}^{p}$, $\bm{\mu}_{3}=(0,0,3,0,\ldots,0)^{\top}\in\mathbb{R}^{p}$, $\bm{\Sigma}=\mathbf{I}$, where $\mathbf{I}$ is the identity matrix of size $p$. Then,

$$\mathbf{A}_{n}=\big{(}\underbrace{\bm{\mu}_{1},\ldots,\bm{\mu}_{1}}_{n_{1}},\underbrace{\bm{\mu}_{2},\ldots,\bm{\mu}_{2}}_{n_{2}},\underbrace{\bm{\mu}_{3},\ldots,\bm{\mu}_{3}}_{n_{3}}\big{)},$$

where $n_{1}=n_{2}=0.3n,~{}n_{3}=0.4n$. Therefore, the true number of clusters is $K=3$.

Case 3. The same setting as in the above Case 2 with $\bm{\Sigma}=(\sigma_{i,j})_{p\times p}$ instead of $\mathbf{I}$, where $\sigma_{i,j}=0.2^{|i-j|}$.

The spikes in the above cases are bounded. We also consider a case of spikes with $\gamma_{K}\rightarrow\infty$ at a rate faster than $\log n$ and $\gamma_{1}=O(p)$.

Case 4. Let $\bm{\mu}_{1}=(2a,a,-a,a,1,\ldots,1)^{\top}\in\mathbb{R}^{p}$, $\bm{\mu}_{2}=(a,a,2a,-3a,1,\ldots,1)^{\top}\in\mathbb{R}^{p}$, $\bm{\mu}_{3}=(a,-2a,-a,a,1,\ldots,1)^{\top}\in\mathbb{R}^{p}$, $\bm{\mu}_{4}=(-2a,a,a,a,1,\ldots,1)^{\top}\in\mathbb{R}^{p}$, and the sample size of cluster corresponding to each center be $n_{1}=n_{3}=0.3n,~{}n_{2}=n_{4}=0.2n$, such that the true number of clusters $K=4$. Suppose $\bm{\Sigma}=(\sigma_{i,j})_{p\times p}$, where $a=\sqrt{p/10}$, $\sigma_{i,j}=0.2^{|i-j|}$.

Tables 1 to 4 report the percentages of under-estimated, exactly estimated and over-estimated under 1000 replications. From the reported results, we see the criteria based on EDA and EDB work better and better as $n,p$ become larger. When $c=1/3$, the probabilities of the under-estimated number of clusters are equal to $0$ and increase when $c$ is getting closer to $1$. From (25), it is shown that the larger $c$ is, the harder the gap conditions are to be satisfied. EDB generally outperforms EDA except the case of $c=3/4$, when $p,n$ are large. It can be seen that when $c=3/4$, as $n$ increases, the probability of $\mathcal{F}_{-}$ estimated by EDB becomes larger, and is uniformly greater than that by EDA. This is due to the fact that the coefficient in the penalty term of EDB criterion is “$\log n$” which is different from the coefficient ”$2$” in EDA, so that the gap condition of EDB is more stronger than of EDA, that is, (27) is more difficult to be satisfied than (26). The criteria based on EDA and EDB show the highest accuracy under Bernoulli distribution, followed by normal, $t_{8}$ and $\chi^{2}(3)$ with relatively heavy right tail which may be destructive to the results.