3.1 Spectral clustering for the latent class model

Our primary goal is to estimate the latent label vector ${\boldsymbol s}^\star $ . However, maximizing the joint likelihood (3) directly still suffers from computationally intractable for large N and J. We first consider spectral clustering, which is a computationally efficient method that leverages the approximate low-rank structure of the data matrix. Let $\mathbf {Z}=(Z_{i,j})\in \left \{0,1\right \}^{N\times K}$ be a latent class membership matrix with $Z_{i,k}=1$ if $s^\star _i=k$ and $Z_{i,l}=0$ otherwise. Under the fixed-effect LCM, we have $R_{i,j}\sim \text {Bern}\left (\sum _{k=1}^{K}Z_{i,k}\theta _{j,k}\right )$ and we can write $R_{i,j}=\sum _{k=1}^{K}Z_{i,k}\theta _{j,k}+E_{i,j}$ , where $\sum _{k=1}^{K}Z_{i,k}\theta _{j,k} = \mathbb E[R_{i,j}]$ and $E_{i,j}$ is mean-zero random noise. We have the following crucial decomposition of the data matrix $\mathbf {R}$ into two additive components using matrix notation:

(5) $$ \begin{align} \mathbf R = \underbrace{\mathbf{Z}}_{N\times K} \underbrace{{\boldsymbol{\Theta}}^\top}_{K\times J} + \underbrace{\mathbf{E}}_{\text{noise}}. \end{align} $$

Here, $\mathbf {Z}\boldsymbol {\Theta }^\top = \mathbb E[\mathbf {R}]$ is the “signal” part containing the crucial latent class information (and thus is low-rank as $K\ll \min \left \{N,J\right \}$ ), and $\mathbf {E}=(E_{i,j})\in \mathbb {R}^{N\times J}$ is a mean-zero “noise” matrix with $\mathbb E[\mathbf {E}] = \mathbb E[\mathbf {R} - \mathbb E[\mathbf {R}]] = \mathbf {0}_{N\times J}$ .

Given the above observation (5), we first briefly outline the rationale of spectral clustering under the oracle noiseless setting (i.e., when $\mathbf {E}=\mathbf {0}$ ). Let $\mathbf {Z}\boldsymbol {\Theta }^\top = \mathbf {U}\boldsymbol {\Sigma }\mathbf {V}^\top $ be the rank-K SVD of $\mathbf {Z}\boldsymbol {\Theta }^\top $ . Notably, the ith row vector in $\mathbf {U}\boldsymbol {\Sigma }$ can be written as

(6) $$ \begin{align} {\boldsymbol e}_i^\top\mathbf{U}\boldsymbol{\Sigma}={\boldsymbol e}_i^\top\mathbf{Z}\boldsymbol{\Theta}^\top\mathbf{V}=\sum_{k=1}^KZ_{i,k}\boldsymbol{\Theta}_{:,k}^\top\mathbf{V}=\boldsymbol{\Theta}_{:,s_i}^\top\mathbf{V},\quad \forall i\in[N]. \end{align} $$

Thus, data points in the same latent class have identical row vectors in $\mathbf {U}\boldsymbol {\Sigma }$ , namely, one of the K distinct vectors $\left \{\boldsymbol {\Theta }_{:,k}^\top \mathbf {V}\right \}_{k=1}^K$ . Hence, the true latent class label information can be readily read off from rows of $\mathbf {U}\boldsymbol {\Sigma }$ . In the realistic noisy data setting, denote by $\widehat {\mathbf {U}}\widehat {\boldsymbol {\Sigma }}\widehat {\mathbf {V}}^\top $ the top-K SVD of the data matrix $\mathbf {R}$ . Then intuitively, $\mathbf {U}\boldsymbol {\Sigma }$ would be well-approximated by $\widehat {\mathbf {U}}\widehat {\boldsymbol {\Sigma }}$ provided that $\mathbf {E}$ is relatively small in magnitude, and the row vectors of $\widehat {\mathbf {U}}\widehat {\boldsymbol {\Sigma }}$ are noisy point clouds distributed around the K points. Therefore, K-means clustering can deliver reasonable estimates of class labels. We illustrate this point using an numerical example in Figure 1, where the row vectors of $\mathbf {U}\boldsymbol {\Sigma }$ (left) and $\widehat {\mathbf {U}}\widehat {\boldsymbol {\Sigma }}$ (right) are plotted as points in $\mathbb R^2$ , respectively.

Figure 1

An illustration for spectral clustering: row vectors of $\mathbf {U}\boldsymbol {\Sigma }$ (left) and $\widehat {\mathbf {U}}\widehat {\boldsymbol {\Sigma }}$ (right). Setting: $N=500$ , $J=250$ , and $K=3$ .

The above high-level idea can be formalized as the spectral clustering algorithm (Algorithm 1), investigated in Zhang and Zhou (Reference Zhang and Zhou2024) for sub-Gaussian mixture models. Algorithm 1 first obtains the rank-K SVD of $\mathbf {R}$ denoted by $\widehat {\mathbf {U}}\widehat {\boldsymbol {\Sigma }}\widehat {\mathbf {V}}^\top $ , and then performs K-means clustering on the rows of $\widehat {\mathbf {U}}\widehat {\boldsymbol {\Sigma }}$ , the left singular vectors weighted by the singular values.

Algorithm 1:

3.2 One-step likelihood refinement with spectral initialization

Despite computationally efficient, spectral clustering itself ignores the likelihood information in the LCM that might improve statistical accuracy. This motivates our proposed Algorithm 2 named spectral clustering with one-step likelihood refinement algorithm, or SOLA for short.

In particular, the first stage is to perform spectral clustering to obtain an initial class label estimate $\widetilde {\boldsymbol s}$ . This stage utilizes the inherent approximate low-rank structure of data but does not use the likelihood information. In the second stage, recall the joint likelihood function defined in (3) as

$$ \begin{align*} L\left({\boldsymbol s},{\boldsymbol{\Theta}}\mid \mathbf R\right) = \sum_{i\in[N]} \sum_{j\in[J]} \Big(R_{i,j} \log\left(\theta_{j,s_i}\right) + (1-R_{i,j}) \log\left(1-\theta_{j,s_i}\right)\Big). \end{align*} $$

We then perform the following two steps:

1. 1. (Maximization) For fixed $\widetilde {\boldsymbol s}$ , find the $\boldsymbol {\Theta }$ maximizing the objective $ L\left (\widetilde {\boldsymbol s},{\boldsymbol{\Theta}} \mid \mathbf R\right )$ , which turns out to have a explicit solution as

   $$ \begin{align*}\widehat\theta_{j,k}= \frac{\sum_{i\in[N]} R_{i,j} \mathbb{I}\left(\widetilde s_i=k\right)}{\sum_{i\in[N]} \mathbb{I}\left(\widetilde s_i=k\right)},\quad \forall j\in[J],~ k\in[K]. \end{align*} $$

   In other words, this is the sample average of the j-th column of $\mathbf {R}$ based on estimated class labels $\left \{i\in [N]:\widetilde s_i=k\right \}$ .

2. 2. (Assignment) Given $\widehat {\boldsymbol {\Theta }}$ , the objective function $ L({\boldsymbol s},\widehat {\boldsymbol {\Theta }}\mid {\mathbf R})$ can be decoupled into N terms, each corresponding to an individual data point, and it is equivalent to assign label $k\in [K]$ to sample i which maximize the log-likelihood independently for each $i\in [N]$ :

   $$ \begin{align*} \widehat s_i = \operatorname*{\mbox{arg max}}_{k\in[K]} \sum_{j\in[J]} \Big[R_{i,j} \log\Big(\widehat\theta_{j,k}\Big) + (1-R_{i,j}) \log\Big(1-\widehat\theta_{j,k}\Big)\Big],\quad i\in[N]. \end{align*} $$

It is readily seen that these two steps increase the joint likelihood monotonically. However, convergence to true ${\boldsymbol s}^\star $ depends on the performance of the initial estimator $\widetilde {\boldsymbol s}$ in the first stage, which we obtained using spectral clustering. The above details are summarized in Algorithm 2, and we provide its algorithmic guarantee in Section 4.

Algorithm 2:

Interestingly, our two-stage procedure resonates with the spirit of the “initialize-and-refine” scheme that has been used in network analysis for stochastic block models (SBMs) (Chen et al., Reference Chen, Liu and Ma2022; Gao et al., Reference Gao, Ma, Zhang and Zhou2017; Lyu et al., Reference Lyu, Li and Xia2023), while the model considered here differs fundamentally in that LCMs handle high-dimensional response data without relying on pairwise interaction structures inherent to SBMs.

In practice, we observe that performing multiple likelihood-based refinement steps may lead to improved convergence and performance compared to the one-step refinement procedure. Thus, we also consider a multiple-step refinement variant of Algorithm 2, which we denote by SOLA+, in the later simulation studies. We set the number of “Maximization $+$ Assignment” refinement steps for SOLA+ to be $10$ in all our simulations. In principle, one could also adapt the algorithm by monitoring the change in the joint log-likelihood values and stopping the algorithm once the absolute difference between successive iterations falls below a small threshold (Di Mari et al., Reference Di Mari, Bakk, Oser and Kuha2023).

3.3 Extension: Classification EM refinement

Our proposed framework can be readily extended to one-step CEM refinement, where CEM is another widely used approach to estimate LCMs (Celeux & Govaert, Reference Celeux and Govaert1992; Zeng et al., Reference Zeng, Gu and Xu2023). It can be regarded as a variant of the standard EM algorithm designed specifically for clustering. Unlike traditional EM, which calculates “soft” assignments by estimating the probability that each observation belongs to every latent class, CEM uses “hard” assignments: each observation is directly assigned to the class with the highest posterior probability. By explicitly modeling the mixture proportion parameters for the latent classes, CEM adjusts the likelihood to account for varying cluster sizes. Concretely, in our setting, one can adapt the one-step refinement in Algorithm 2 to a CEM step by including class proportions. Besides updating $\widehat {\theta }_{j,k}$ as in (7), we also estimate the proportion parameters $\widehat {{\boldsymbol p}} = (\widehat p_1,\ldots ,\widehat p_K)\in [0,1]^K$ by

(9) $$ \begin{align} \widehat{p}_k \;=\; \frac{1}{N} \sum_{i=1}^N \mathbb{I}\{\widetilde{s}_i = k\}, \quad k \,\in\, [K]. \end{align} $$

Then, for each individual i, we update their latent class label based on both the item parameters and the estimated class proportions:

(10) $$ \begin{align} \widehat s_i \;=\; \arg\max_{k \in [K]} \Biggl[\log \widehat p_k \;+\; \sum_{j=1}^J \Bigl(\,R_{i,j}\,\log\widehat{\theta}_{j,k} \;+\; (1 - R_{i,j})\,\log\bigl(1 - \widehat{\theta}_{j,k}\bigr)\Bigr) \Biggr]. \end{align} $$

This one-step CEM update replaces the likelihood-refinement stage of our two-stage estimator: start with spectral clustering, then perform the CEM step above. We provide its theoretical guarantee in Section 4. In practice, one could easily generalize the above CEM procedure to a multiple-step refinement variant, which is analogous to SOLA+.

4.1 Theoretical guarantee of spectral clustering

In this section, we introduce the theoretical guarantee of spectral clustering (Algorithm 1). To that end, we define some necessary notations. The difficulty of recovering $\mathbf {s}^\star $ depends on how distinct the latent classes are, which, in the context of spectral clustering, is characterized by the following quantity:

(11) $$ \begin{align} \Delta:=\min_{k_1\ne k_2\in[K]}\left\|\boldsymbol{\Theta}_{:,k_1}-\boldsymbol{\Theta}_{:,k_2}\right\| = \min_{k_1\ne k_2\in[K]}\left[\sum_{j=1}^J (\theta_{j,k_1} - \theta_{j,k_2})^2\right]^{1/2}. \end{align} $$

This quantity measures the minimal Euclidean distance between any two class-specific item parameter vectors. Intuitively, a larger $\Delta $ yields a larger gap between the columns of $\boldsymbol {\Theta }$ and a bigger distinction among the latent classes, and therefore produces clearer separation in the spectral embeddings (referring to (6)), making recovery of the true class labels easier. In addition, let $\lambda _1\ge \cdots \ge \lambda _K>0$ denote the singular values of $\mathbb {E}[\mathbf {R}]=\mathbf {Z}\boldsymbol {\Theta }^\top $ in decreasing order. The following result characterizes the clustering error of Algorithm 1.

Proposition 1 (Adapted from Theorem 3.1 in Zhang & Zhou, Reference Zhang and Zhou2024).

Define $\sigma ^2_{\theta }:=\left (1-2\theta \right )/\left [2\log \left (\left (1-\theta \right )\theta ^{-1}\right )\right ]$ for $\theta \in (0,1)\backslash \left \{1/2\right \}$ and $\sigma ^2_{\theta }:=1/4$ for $\theta =1/2$ , and let $\bar \sigma :=\max _{j,k}\sigma _{\theta _{j,k}}$ . Assume $\text {rank}\left (\boldsymbol {\Theta }\right )=K$ , $\alpha N/K^2 \ge 10$ and

(12) $$ \begin{align} \min\left\{\frac{\Delta}{\alpha^{-1/2}K\left(1+\sqrt{\frac{J}{N}}\right)\bar\sigma},\frac{\lambda_K}{\left(\sqrt{N}+\sqrt{J}\right)\bar\sigma}\right\}\rightarrow\infty. \end{align} $$

Then, for the spectral clustering estimator $\widetilde {\boldsymbol s}$ in Algorithm 1, we have

(13) $$ \begin{align} \mathbb{E} \ell\left(\widetilde {\boldsymbol s},{\boldsymbol s}^\star\right)\le \exp\left(-\frac{\Delta^2}{8\bar\sigma^2}\left(1-o(1)\right)\right)+\exp\left(-\frac{N}{2}\right), \end{align} $$

where $\ell (\cdot ,\cdot )$ is defined in (4). Moreover, if further assume that $\Delta ^2\ge 8\left (1+c\right )\bar \sigma ^2\log N$ for any constant $c>0$ , then we have $\ell \left (\widetilde {\boldsymbol s},{\boldsymbol s}^\star \right )=0$ with probability $1-\exp \left (-N/2\right )-\exp \left (-\Delta /\bar \sigma \right )$ .

Proposition 1 delivers that the simple and computationally efficient spectral clustering can essentially lead to exponentially small error rate with respect to the latent class separation $\Delta $ . Notably, we can have exact recovery of all latent class labels (i.e., $\ell \left (\widetilde {\boldsymbol s},{\boldsymbol s}^\star \right )$ with $\widetilde s_i=\pi (s_i)$ for all $i=1,\ldots ,N$ ) with high probability as long as $\Delta ^2/\bar \sigma ^2\ge 8\left (1+c\right )\log N$ . In the LCM we study, we have $\bar \sigma \le 1/4$ by definition, and hence it suffices to have $\Delta ^2\ge 2(1+c)\log N$ to achieve exact recovery. This is not hard to achieve even for slowly growing number of items, because if we have constant separation of item parameters for each item with $\delta = \min _{k_1\neq k_2\in [K],j\in [J]} |\theta _{j,k_1} - \theta _{j,k_2}|$ , then $\Delta ^2 = J\delta ^2$ . In this case, having $J> 2(1+c)\delta ^{-2} \log N$ for exact recovery corresponds to a moderately growing dimension. In sharp contrast, the JML estimator analyzed in Zeng et al. (Reference Zeng, Gu and Xu2023) only achieves a polynomial convergence rate of order $O (J^{-1/2})$ . Moreover, they cannot achieve exact recovery according to the technical condition and conclusion therein. Therefore, the new spectral estimator $\widetilde {\boldsymbol s}$ represents a substantial improvement with exponential error decay in (13).

We also note that the condition (12) on $\Delta $ and $\lambda _K$ is mild and can be easily satisfied in the modern high-dimensional data setting. We next give a concrete example. Let us consider a common generative setting for item parameters $\boldsymbol {\Theta }$ , where its entries are i.i.d. generated from the Beta distribution.

Proposition 2. Consider $\theta _{j,k} \stackrel {i.i.d.}{\sim } \text {Beta}(a,b)$ for all $j \in [J]$ and $k \in [K]$ with $a,b>0$ . Define

$$ \begin{align*}B:=\frac{a\left(a+b+a\left(a+b\right)^{-1}\right)}{\left(a+b\right)\left(a+b+1\right)},\end{align*} $$

then, under Assumption 2, we have $\lambda _K\gtrsim \sqrt {\alpha B NJ}$ and $\Delta \gtrsim \sqrt {BJK}$ with probability at least $1-K(e^{-c_1B^2J}+e^{-c_2 BJ})$ for some universal constants $c_1,c_2\in (0,1)$ .

Proposition 2 entails that the condition on $\Delta $ and $\lambda _K$ for the success of spectral clustering in Proposition 1 can be easily satisfied when $\boldsymbol {\Theta }$ consists of Beta distributed entries, a common assumption in Bayesian inference for the LCM. Specifically, combined with Assumption 1, a sufficient condition on B for spectral clustering to have the desirable error control (13) would be $B\gg K(\frac {1}{N}+\frac {1}{J})$ , which is easily satisfied if the Beta parameters a and b are of the constant order.

Despite the above nice theoretical guarantee for spectral clustering, the clustering error’s reliance on the Euclidean metric $\Delta $ in Proposition 1 is inherent from the general framework analyzed in Zhang and Zhou (Reference Zhang and Zhou2024), and it turns out to be not statistically optimal due to ignorance of the likelihood information. We will discuss and address this issue next.

4.2 Fundamental statistical limit for LCMs

In this section, we provide the understanding of the information-theoretic limit (fundamental statistical limit) of clustering in LCMs. To compare different statistical procedures, we consider the minimax framework that had been widely considered in both statistical decision theory and information theory in the past decades (Le Cam, Reference Le Cam2012; Takezawa, Reference Takezawa2005). The minimax risk is defined as the infimum, over all possible estimators, of the maximum loss (here, the mis-clustering error) taken over the entire parameter space. In other words, it is the rate of the best estimator one can have among all statistical procedures in the worst-case scenario.

Formally, we establish the minimax lower bound for mis-clustering error under the fixed-effect LCM in the following result, whose proof is adapted from Lyu et al. (Reference Lyu, Chen and Gu2025) and can be found in Section B.1.

Theorem 1 (Minimax lower bound).

Consider the following parameter space for $\mathbb {E}[\mathbf {R}]$ :

$$ \begin{align*} &{\mathcal{P}}_K\left({\boldsymbol s},\boldsymbol{\Theta}\right):=\left\{\bar{\mathbf{R}}:\bar R_{i,j}=\theta_{j,s_i}, {\boldsymbol s}\in[K]^{N}, \theta_{j,k}\in[0,1] , \forall i\in[N], j\in[J], k\in[K]\right\}. \end{align*} $$

Define $I^\star :=\min _{k_1\ne k_2\in [K]}\sum _{j\in [J]}I (\theta _{j,k_1},\theta _{j,k_2} )$ with

$$ \begin{align*} I(p,q):=-2\log\left(\sqrt{pq}+\sqrt{\left(1-p\right)\left(1-q\right)}\right), \quad \forall p,q\in(0,1), \end{align*} $$

then we have

(14) $$ \begin{align} \inf_{\widehat{\boldsymbol s} }\sup_{{\mathcal{P}}_K\left({\boldsymbol s},\boldsymbol{\Theta}\right)}\mathbb{E}\ell\left(\widehat {\boldsymbol s},{\boldsymbol s}^\star\right)\ge \exp\left(-\frac{I^\star}{2}\left(1+o(1)\right)\right), \end{align} $$

where the infimum is taken over all possible estimators of the latent class labels.

Roughly speaking, Theorem 1 informs us that under the fixed-effect LCM, the error rate $\exp (-I^\star /2 )$ is the lowest possible clustering error rate that no estimator can uniformly surpass. Notably, $I(p,q)$ is exactly the Renyi divergence (Rényi, Reference Rényi1961) of order $1/2$ between two Bernoulli distributions, Bern $(p)$ and Bern $(q)$ , and hence $\sum _{j\in [J]}I(\theta _{j,k_1},\theta _{j,k_2})$ can be interpreted as the Renyi divergence between two Bernoulli random vectors with parameters $\boldsymbol {\Theta }_{:,k_1}$ and $\boldsymbol {\Theta }_{:,k_2}$ , respectively. $I^\star $ can be regarded as the overall signal-to-noise ratio (SNR) which further takes account into all possible pairs of $(k_1,k_2)\in [K]^2$ with $k_1\ne k_2$ . Although the form of $I^\star $ might seem complicated at its first glance, it actually exactly reflects the likelihood information in the LCMs. By definition of the minimax lower bound, we know that $\exp (-I^\star /2 )$ cannot be larger than the spectral clustering error rate $\exp (-\Delta ^2/(8\bar \sigma ^2) )$ obtained in Proposition 1. We will discuss this in more detail in Section 4.4.

In the following, we take a closer look at why the quantity $I^\star $ shows up as the fundamental statistical limit. We analyze an idealized scenario where the true item parameters $\boldsymbol {\Theta }$ are known, which we term as the oracle setting. For simplicity, we consider $K=2$ latent classes and $I^\star =\sum _{j\in [J]}I(\theta _{j,1},\theta _{j,2})$ . In this regime, we can relate the clustering problem for the i-th sample to the following binary hypothesis testing problem:

$$ \begin{align*} H_0:s_i^\star=1,\quad \text{v.s.} \quad H_0:s_i^\star=2. \end{align*} $$

By the Neyman–Pearson lemma, the test that gives the optimal Type-I plus Type-II error is the likelihood ratio test. Consequently, the oracle classifier assigns an observation to the class that maximizes the log-likelihood:

(15) $$ \begin{align} \bar s_i = \operatorname*{\mbox{arg max}}_{k\in[2]} \sum_{j\in[J]} \Big[R_{i,j} \log\theta_{j,k} + (1-R_{i,j}) \log\left(1-\theta_{j,k}\right)\Big]. \end{align} $$

Assuming without loss of generality that $s^\star _i=1$ , one can show that the mis-clustering probability satisfies

$$ \begin{align*} &\mathbb{P}\left(\bar s_i\ne s^\star_i\right) =\mathbb{P}\left(\sum_{j\in[J]}R_{i,j}\log\frac{\theta_{j,2}\left(1-\theta_{j,1}\right)}{\theta_{j,1}\left(1-\theta_{j,2}\right)}\ge \sum_{j\in[J]}\log\frac{1-\theta_{j,1}}{1-\theta_{j,2}}\right)\le \exp\left(-\frac{\sum_{j\in[J]}I(\theta_{j,1},\theta_{j,2})}{2}\right). \end{align*} $$

In other words, the above error rate in the oracle setting directly informs the information-theoretic limit of any clustering procedure in the LCM with binary responses, as indicated by Theorem 1. Note that the oracle estimator (15) that achieves this optimal rate is constructed using the true $\boldsymbol {\Theta }$ and cannot be applied in practice. Later, we will construct a computationally feasible estimator (i.e., our SOLA estimator) that achieves this optimal rate without relying on prior knowledge of $\boldsymbol {\Theta }$ .

4.3 Theoretical guarantee of SOLA

In this section, we give the theoretical guarantee of SOLA. To facilitate the theoretical analysis, we consider a sample-splitting variant of SOLA (Algorithm 3), where the sample data points are randomly partitioned into two subsets ${\mathcal{S}}_1,{\mathcal{S}}_2$ of size $\lceil N/2\rceil $ and $N-\lceil N/2\rceil $ . The key idea is to estimate $\boldsymbol {\Theta }$ using one subset of samples, and then update labels for the remaining ones. The procedure can be repeated by switching the roles of the two subsets, and we can obtain labels for all samples. Finally, an alignment step is needed due to the permutation ambiguity of clustering in the two subsets. We have the following result, whose proof is deferred to Section B.2.

Algorithm 3:

Theorem 2. Suppose Assumptions 1 and 2 and the assumptions of Proposition 1 hold. In addition, assume that $N\gtrsim K\log N$ , $JK\lesssim \min \{\exp (-c\Delta ^2/\bar \sigma ^2 ),N^{\gamma _0}\}$ , for some universal constants $0<c<1/8$ and $\gamma _0>0$ . If the SNR condition satisfies

$$ \begin{align*} \frac{\Delta/\bar\sigma}{\sqrt{J}K^{3/4}\left(\log N\right)^{1/4}/N^{1/4}}\rightarrow\infty, \end{align*} $$

then, for $\widehat {\boldsymbol s}$ from Algorithm 3, we have

$$ \begin{align*} \mathbb{E}\ell \left(\widehat {\boldsymbol s}, {\boldsymbol s}^\star\right)\le \exp\left(-\frac{I^\star}{2}\left(1-o(1)\right)\right). \end{align*} $$

Combining Theorems 1 and 2, we can see that the estimator of Algorithm 3 achieves the fundamental statistical limit and is optimal for estimating the latent class labels. We have several remarks on Theorem 2. First, it is interesting to consider the regime when $J=o (\sqrt {\frac {N\log N}{K}} )$ , since otherwise the SNR condition would already implies exact recovery (perfect estimation of latent class labels with high probability) for spectral clustering, as shown by Theorem 2 and Proposition 1. Second, the sample splitting step is crucial in the current proof of Theorem 2 to decouple the dependence between $R_{i,j}$ and $\widehat {\boldsymbol {\Theta }}$ in (8). This step serves purely as a theoretical device to facilitate the theoretical analysis of clustering error rate, as commonly adopted in theoretical statistics (Abbe, Reference Abbe2018; Lei & Zhu, Reference Lei and Zhu2017; Rinaldo et al., Reference Rinaldo, Wasserman and G’Sell2019; Vu, Reference Vu2018). In practice, halving the sample only increases the estimation error of $\widehat {\boldsymbol {\Theta }}$ by a constant factor, which is negligible in theory and immaterial in practice when the sample size N is large. Alternatively, it may be possible in the future to show that the same theoretical guarantee holds for the original version of SOLA, that is, Algorithm 2, by developing a leave-one-out analysis similar to that in Zhang and Zhou (Reference Zhang and Zhou2024). However, the analysis will be highly technical and much more involved than that in Zhang and Zhou (Reference Zhang and Zhou2024) and beyond the scope of our current focus. Last, although sample splitting simplifies theoretical analysis, it may not be necessary in practice; we use the original SOLA (Algorithm 2) for all numerical experiments and observe satisfactory empirical results.

As an extension of Theorem 2, we have a similar theoretical guarantee for the one-step CEM refinement proposed in Section 3.3, whose proof is given in Section B.4.

Corollary 1. Suppose the conditions in Corollary 2 hold. Consider the output of Algorithm 3 with (16) replaced by (9) and sample splitting version of (10), we have

$$ \begin{align*} \mathbb{E}\ell \left(\widehat {\boldsymbol s}, {\boldsymbol s}^\star\right)\le \exp\left(-\frac{I^\star}{2}\left(1-o(1)\right)\right). \end{align*} $$

Theorem 1 indicates that one-step CEM refinement achieves the same exponential error rate $\exp (-I^*/2)$ as Algorithm 3. There might be a gap between these two methods when Assumption 2 does not hold, which is an interesting future direction to consider.

4.4 Comparison of spectral clustering and SOLA

In this section, we compare the convergence rate of spectral clustering and SOLA under LCM. While spectral clustering is widely employed due to its computational efficiency and solid theoretical guarantees in general clustering settings, we demonstrate that SOLA attains minimax-optimal clustering error rates in LCMs that spectral clustering fails to achieve.

In view of Theorem 1, the minimax-optimal clustering error rate for LCM is of order $\exp (-I^\star /2 )$ . On the other hand, spectral clustering is shown to achieve the error rate scaling as $\exp (-\Delta ^2/(8\bar \sigma ^2) )$ in Proposition 1. To understand the difference between these two quantities, we have the following result, whose proof can be found in Section B.5.

Proposition 3. Assume that (a) $\min _{k_1\ne k_2\in [K]}\left |\Omega _0(k_1,k_2)\right |\gtrsim J$ , where $\Omega _0(k_1,k_2):= \{j\in [J]:\left |\theta _{j,k_1}-\theta _{j,k_2}\right |=o(1) \}$ and (b) $\min _{k\in [K]}\left |\left \{j\in [J]:\left |\theta _{j,k}-1/2\right |>c_0\right \}\right |\gtrsim J$ for some universal constant $c_0>0$ . Then we have

$$ \begin{align*} \exp\left(-\frac{I^\star}{2}\right) &\lesssim\exp\left(-\frac{\Delta^2}{8\bar\sigma^2}\left(1+c\right)\right), \end{align*} $$

for some universal constant $c>0$ .

A few comments are in order regarding Proposition 3. First, under the conditions (a) and (b), the rate of spectral clustering becomes sub-optimal (in order) due to a loose constant in the exponent as the latent class separation $\Delta \rightarrow \infty $ . Second, it is important to note that the mis-clustering error $\ell (\widehat {\boldsymbol s}, {\boldsymbol s}^*)$ takes values in the finite set $\left \{0,1/N,\ldots ,1\right \}$ according to the definition (4). Given this, it is only meaningful to compare error rates when the exponent $\Delta ^2/\bar \sigma ^2$ (or $I^\star $ ) does not exceed a constant factor of $\log N$ , since otherwise exact clustering of all subjects is achieved with high probability. To see this, we use Markov inequality to get

$$ \begin{align*}\mathbb{P}(\ell(\widehat{\boldsymbol s}, {\boldsymbol s}^*) \neq 0) = \mathbb{P}(\ell(\widehat{\boldsymbol s}, {\boldsymbol s}^*) \ge 1/N) \leq N \mathbb{E}\ell(\widehat{\boldsymbol s}, {\boldsymbol s}^*) \leq N \exp\left(-C\log N\right) \to 0, \end{align*} $$

for a constant $C>1$ . In view of this, the condition (a) helps us exclude the trivial case where $I^\star \asymp \Delta ^2\asymp J$ , which would only allow $J\lesssim \log N$ . Third, the essential order difference occurs when most of $\theta _{j,k}$ ’s are bounded away from $1/2$ . To see this, consider the scenario when $K=2$ , $\theta _{j,1}=0.5$ , and $\theta _{j,2}=0.5-\delta $ for some $\delta =o(1)$ . Then we have $\Delta =\frac {1}{2}J\delta ^2$ and $I^\star =-J\log \left (\sqrt {\frac {1}{2}\left (\frac {1}{2}-\delta \right )}+\sqrt {\frac {1}{2}\left (\frac {1}{2}+\delta \right )}\right )=\frac {1}{2}J\delta ^2\left (1+o(1)\right )$ . This implies that without condition (b), the exponents $\Delta $ and $I^\star $ are of the same order, resulting in at most a constant-factor difference in clustering error for spectral clustering and SOLA. We also conduct numerical experiments in Section 5 to verify the necessity of (b).

4.5 Estimation of the number of latent classes K

In practice, the true number of latent classes K may be unknown. To see how we can estimate K from data, recall that $\mathbf {R}=\mathbf {Z}\boldsymbol {\Theta }^\top +\mathbf {E}$ . By Weyl’s inequality and classical random-matrix theory (Vershynin, Reference Vershynin2018), the singular values of $\mathbf {R}$ satisfy

$$ \begin{align*}\left|\lambda_i\left(\mathbf{R}\right)-\lambda_i\left(\mathbf{Z}\boldsymbol{\Theta}^\top\right)\right|\le \left\|\mathbf{E}\right\|=O_p\left(\sqrt{J}+\sqrt{N}\right),\quad i\in[N\wedge J].\end{align*} $$

Since $\mathbf {Z}\boldsymbol {\Theta }^\top $ has exactly K nonzero singular values $\lambda _1\ge \cdots \ge \lambda _K>0$ , and Proposition 2 shows that $\lambda _K\gtrsim \sqrt {\alpha BNJ}$ under a common generative setting, the top K singular values of $\mathbf {R}$ lie well above the noise level. Accordingly, we estimate K by counting how many singular values of $\mathbf {R}$ exceed a threshold slightly above the noise level:

(17) $$ \begin{align} \widehat K=\left|\left\{i\in[N\wedge J]:\lambda_i\left(\mathbf{R}\right)>2.01 \left(\sqrt{J}+\sqrt{N}\right)\right\}\right|. \end{align} $$

Here, the factor $2.01$ ensures we count only those singular values contributed by the low-rank signal. We have the following lemma certifying the consistency of $\widehat K$ , whose proof is deferred to Section B.6.

Lemma 1. Instate the assumptions of Proposition 1. For $\widehat K$ defined in (17), we have

$$ \begin{align*}\mathbb{P}\left(\widehat K=K\right)=1-o(1),\qquad \text{as~}N,J\rightarrow\infty.\end{align*} $$

This thresholding rule provides a formal analog of the classical scree-plot method (Cattell, Reference Cattell1966; Zhang et al., Reference Zhang, Chen and Li2020), and its proof shares a similar spirit of Theorem 2 in Zhang et al. (Reference Zhang, Chen and Li2020) but with a sharper characterization of the constant factor in the noise magnitude.

5.1 Simulation studies

We compare our proposed methods SOLA (Algorithm 2) and $\text {SOLA}+$ (Algorithm 2 with multiple steps refinement), with spectral method (Zhang & Zhou, Reference Zhang and Zhou2024), EM (Linzer & Lewis, Reference Linzer and Lewis2011), and Tensor-CEM (Zeng et al., Reference Zeng, Gu and Xu2023). For the EM algorithm, we employ the popular polca package (Linzer & Lewis, Reference Linzer and Lewis2011) in R, and for Tensor-CEM, we use the original MATLAB code (Zeng et al., Reference Zeng, Gu and Xu2023) for tensor initialization combined with our R-based implementation of the CEM refinement. Our evaluation focuses on three aspects: (a) clustering accuracy, (b) stability, and (c) computational efficiency. Throughout the simulations, we generate the true latent label $\left \{s_i^\star \right \}_{i=1}^N$ independently uniformly from $[K]$ and set $K=3$ , $N=2J$ . The parameters $\theta _{j,k}$ are independently sampled from Beta $(\beta _1,\beta _2)$ . We generate $200$ independent replicates and report the average mis-clustering error and computation time, excluding any replicates in which the algorithm fails (Tensor-CEM and EM). Notably, our algorithm is numerically stable without any failures in all simulations.

Remark 3 (EM implementations: polca vs. LatentGOLD).

Beyond polca, we also evaluated LatentGOLD (Vermunt & Magidson, Reference Vermunt and Magidson2005), which implements a robust EM routine with safeguards against numerical underflow. Table 1 reports a representative comparison across three regimes ( $N=2J$ , $N=J$ , and $N=J/2$ ), showing that LatentGOLD is substantially more stable and accurate than polca in these settings. Our overall message remains unchanged: EM-type methods can be competitive in the classical regime with relatively large N and smaller J, whereas SOLA/SOLA+ become increasingly advantageous as J grows and the problem enters the high-dimensional regime (with J comparable to or larger than N). While retaining polca for comparison purposes in our large-scale simulation studies, we emphasize that when analyzing data with a small-to-moderate J, it is recommended to use LatentGOLD for their better numerical stability and performance.

Table 1

Clustering error across different methods and settings under $100$ replicates

Simulation 1: Latent class estimation accuracy for dense \ $\boldsymbol {\Theta }$ close to $1/2$ . In the first simulation, we set $(\beta _1,\beta _2)=(5,5)$ . Under these parameters, the entries of $\boldsymbol {\Theta }$ tend to be close to $0.5$ . As shown in Figure 2, our proposed methods achieve lower mis-clustering errors compared to the alternative approaches. In this setting, the improvement of SOLA over the plain spectral clustering method is negligible, which is expected since $\theta _{j,k}$ ’s values near $0.5$ reduce the additional gain from incorporating the likelihood information. This observation supports the necessity of having $\theta _{j,k}$ ’s well separated from $0.5$ (see Proposition A.1) to have likelihood-based refinement improve upon spectral clustering. It is also worth noting that Tensor-CEM performs poorly when both the number of items J and the sample size N are small, likely due to instability introduced by its tensor-based initialization step.

Figure 2

Simulation 1: Mis-clustering proportions versus number of items J. Entries of \ $\boldsymbol {\Theta }$ are independently generated from $\text {Beta} (5,5)$ .

Simulation 2: Latent class estimation accuracy for sparse item parameters $ \ \boldsymbol {\Theta }$ . In the second simulation, we choose $(\beta _1,\beta _2)=(1,8)$ , a sparse scenario, where the $\theta _{j,k}$ ’s tend to be close to $0$ . Figure 3 illustrates that our proposed methods exhibit superior clustering accuracy when N and J are small, while maintaining competitive performance as these dimensions increase. In this sparse setting, the traditional EM algorithm performs the worst because $\theta _{j,k}$ ’s are near the boundary (close to $0$ ), making it difficult for random initializations to converge to optimal solutions. Furthermore, the gap between spectral clustering and SOLA becomes more pronounced when the $\theta _{j,k}$ ’s are bounded away from $0.5$ .

Figure 3

Simulation 2: Mis-clustering proportions versus number of items J. Entries of \ $\boldsymbol {\Theta }$ are independently generated from $\text {Beta} (1,8)$ .

Simulation 3: Stability analysis. We next assess the stability of the different methods under the sparse setting $(\beta _1,\beta _2)=(1,8)$ . In particular, we record a “failure” point if either (a) the estimated number of classes degenerates from the pre-specified K to a smaller value or (b) the algorithm produces an error during execution. As shown in Figure 4, our proposed methods are robust across a wide range of N and J without any failures in any simulation trial. In particular, the failure rate for Tensor-CEM exceeds 15% for small N and J, likely due to its sensitive initialization, while the EM algorithm becomes increasingly unstable as the dimensions grow. These results underscore the robustness of our SOLA approaches.

Figure 4

Simulation 3: Failure rate versus the number of items J. Entries of $\boldsymbol {\Theta }$ are independently generated from $\text {Beta} (1,8)$ .

Simulation 4: Computational efficiency. Finally, we compare the computation time required by each method. In our experiments, the SVD step and our SOLA and $\text {SOLA}+$ algorithms are implemented in R. The total computation time for Tensor-CEM is the sum of the MATLAB initialization and the R refinement stages. As reported in Table 2 and illustrated in Figure 5, our proposed methods are significantly faster than both the EM and Tensor-CEM approaches, with computation time comparable to that of SVD.

Table 2

Simulation 4: Running time (seconds) of different methods

Figure 5

Simulation 4: Running time (seconds) of different methods.

Simulation 5: Latent class estimation accuracy under different regimes. In the fifth simulation study, we also employ Spec-EM (EM algorithm initialized by spectral clustering), motivated by the fact that when the number of items is moderately large, spectral clustering provides an effective initialization for the classical EM algorithm for LCMs. In particular, we consider the following regimes: (i) $N=J$ (Figure 6), (ii) $N=J/2$ (Figure 7), (iii) fix $N=100$ and vary J (Figure 8), (iv) fix $J=100$ and vary N (Figure 9), (v) fix $J=30$ and vary N (Figure 10), and (vi) $N=J$ with $\left \{s_i^\star \right \}_{i=1}^N$ being independently generated from $[3]$ with probability $(p_1,p_2,p_3)=(1/6,1/3,1/2)$ (Figure 11). In all these settings, we set $(\beta _1,\beta _2)=(1,8)$ . In addition, we consider (vii) $N=J$ under the fixed design of $\boldsymbol {\Theta }$ (Figure 12), where the fixed item parameter matrix $\boldsymbol \Theta $ is constructed by vertically concatenating copies of a $5\times 3$ block matrix $\boldsymbol \Theta _0$ as follows:

$$\begin{align*}\boldsymbol\Theta=\begin{bmatrix} \boldsymbol\Theta_0\\ \vdots\\ \boldsymbol\Theta_0 \end{bmatrix}\Bigg\}~b~\text{copies},\qquad \boldsymbol\Theta_0:=\begin{bmatrix} 0.3 & 0.7 & 0.7\\ 0.3 & 0.7 & 0.3\\ 0.7 & 0.3 & 0.7\\ 0.7 & 0.3 & 0.3\\ 0.7 & 0.7 & 0.3 \end{bmatrix}. \end{align*}$$

Here, the block matrix $\boldsymbol \Theta _0$ is vertically stacked b times to yield $J=5b$ items with $K=3$ latent classes. It is clear that SOLA and SOLA+ perform significantly better in most scenarios except the most traditional large-N, small-J regime, where Spec-EM performs better. This underpins both the validity of SOLA and spectral clustering.

Figure 6

Simulation 5-1: Mis-clustering proportions versus number of items J under $N=J$ . Entries of $\boldsymbol {\Theta }$ are independently generated from $\text {Beta}(1,8)$ .

Figure 7

Simulation 5-2: Mis-clustering proportions versus number of items J under $N=0.5J$ . Entries of $\boldsymbol {\Theta }$ are independently generated from $\text {Beta}(1,8)$ .

Figure 8

Simulation 5-3: Mis-clustering proportions versus number of items J under $N=100$ . Entries of $\boldsymbol {\Theta }$ are independently generated from $\text {Beta}(1,8)$ .

Figure 9

Simulation 5-4: Mis-clustering proportions versus sample size N under $J=100$ . Entries of $\boldsymbol {\Theta }$ are independently generated from $\text {Beta}(1,8)$ .

Figure 10

Simulation 5-5: Mis-clustering proportions versus sample size N under $J=30$ . Entries of $\boldsymbol {\Theta }$ are independently generated from $\text {Beta}(1,8)$ .

Figure 11

Simulation 5-6: Mis-clustering proportions versus number of items J under imbalanced latent classes with $N=J$ . Entries of $\boldsymbol {\Theta }$ are independently generated from $\text {Beta}(1,8)$ .

Figure 12

Simulation 5-7: Mis-clustering proportions versus number of items J under $N=J$ . The $\boldsymbol {\Theta }$ matrix is set to a fixed design matrix.

Simulation 6: Estimation of number of latent classes. Our theoretical guarantee for the selection criterion of K relies on random matrix theory under the asymptotic regime $N,J \to \infty $ . This means relatively large N and J are needed for the theoretical guarantee to hold, which is consistent with our focus on high-dimensional data. The constant “2” in the threshold reflects a worst-case bound and can indeed be conservative, especially when entries of $\boldsymbol {\Theta }$ are mostly concentrated away from $1/2$ . To illustrate its effectiveness in practice, we have performed simulation studies with $N=2J$ and $K=3$ , under several generative mechanisms for $\boldsymbol {\Theta }$ (i.i.d. Beta $(1,1)$ , Beta $(1,4)$ , and Beta $(1,5)$ , which reflect different sparsity levels, and also the fixed design as described in Simulation 5). As shown in Table 3, the criterion successfully recovers K in high-dimensional regimes. Our real-data analysis in Section 5.2 also demonstrates that the method works well when the data are sufficiently large. For smaller N and J, we recommend that practitioners use classical information criteria, such as BIC, as a more reliable and interpretable tool for selecting K.

Table 3

Proportion of successfully selecting $\widehat{K}$ to be the true K based on $200$ simulation replicates, under different generative mechanisms of the item parameters $\boldsymbol \Theta =(\theta _{j,k})_{J\times K}$

5.2 Real data example

We further evaluate the performance of our methods on the United States 112th Senate Roll Call Votes data, which is publicly available at https://legacy.voteview.com/senate112.htm. The dataset contains voting records for $J=486$ roll calls by $102$ U.S. senators. Following the preprocessing procedure in Lyu et al. (Reference Lyu, Chen and Gu2025), we convert the response matrix to binary and remove senators who are neither Democrats nor Republicans, as well as those with more than $10\%$ missing votes. This results in a final sample of $N=94$ senators. For the remaining missing entries, we impute values by randomly assigning $0$ or $1$ with probability equal to the positive response rate observed from each senator’s non-missing votes.

We compare the performance of the methods considered in our simulation studies. Table 4 summarizes the mis-clustering error and computation time (in seconds) for each method. Interestingly, the SVD approach yields an error comparable to those obtained with SOLA and SOLA+, suggesting that the refinement step may be unnecessary for this dataset. This observation can be partially explained by our remark following Theorem 2: when the number of items J is substantially larger than the number of individuals N, spectral clustering alone can perform very well. In contrast, the traditional EM method exhibits significantly higher mis-clustering error, and EM and Tensor-CEM suffer from high computational cost in this large-J regime.

Table 4

Mis-clustering error and running time (seconds) of different methods on Senate Roll Call Votes data

We also validate the performance of our method for estimating K on this real dataset. In particular, the threshold defined in (17) is $2.01\times (\sqrt {486}+\sqrt {94})\approx 63.8$ . The first three largest singular values of $\mathbf {R}$ are $\{148.1,\,64.4,\,16.6\}$ , which leads to an estimated number of clusters $\widehat K=2$ . Notably, this estimated value exactly matches the known number of clusters in the dataset (i.e., $K=2$ ).