{'fraction_non_alphanumeric': 0.043760869565217395, 'fraction_numerical': 0.015195652173913044, 'mean_word_length': 4.100454595853199, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Tensor-EM Method for Large-Scale Latent Class Analysis with Clustering Consistency

Latent class models are powerful statistical modeling tools widely used in psychological, behavioral, and social sciences. In the modern era of data science, researchers often have access to response data collected from large-scale surveys or assessments, featuring many items (large J) and many subjects (large N). This is in contrary to the traditional regime with fixed J and large N . To analyze such large-scale data, it is important to develop methods that are both computationally efficient and theoretically valid. In terms of computation, the conventional EM algorithm for latent class models tends to have a slow algorithmic convergence rate for large-scale data and may converge to some local optima instead of the maximum likelihood estimator (MLE). Motivated by this, we introduce the tensor decomposition perspective into latent class analysis. Methodologically, we propose to use a moment-based tensor power method in the first step, and then use the obtained estimators as initialization for the EM algorithm in the second step. Theoretically, we establish the clustering consistency of the MLE in assigning subjects into latent classes when N and J both go to infinity. Simulation studies suggest that the proposed tensor-EM pipeline enjoys both good accuracy and computational efficiency for large-scale data. We also apply the proposed method to a personality dataset as an illustration.


Introduction
Latent class models (LCMs) (Lazarsfeld and Henry, 1968;Goodman, 1974) are powerful statistical modeling tools widely used in psychological, behavioral, and social sciences. LCMs use a categorical latent variable to model the unobserved heterogeneity of multivariate categorical data and identify meaningful latent subgroups of subjects. LCMs have seen broad applications in a variety of scientific fields, including psychology and psychiatry (Bucholz et al., 2000;Keel et al., 2004), sociology and organizational research (Wang and Hanges, 2011), and biomedical and epidemiological studies (Bandeen-Roche et al., 1997;Dean and Raftery, 2010). There are also various extensions and generalizations building upon LCMs, including LCMs with covariates or distal outcomes (Vermunt, 2010;Lanza and Rhoades, 2013), longitudinal LCMs and latent transition analysis (Dunn et al., 2006;Collins and Lanza, 2009), and also restricted LCMs known as diagnostic classification models that involve multiple categorical latent variables (Rupp and Templin, 2008;Xu, 2017;von Davier and Lee, 2019). For general introductions and applications of LCMs, see Hagenaars and McCutcheon (2002) and Collins and Lanza (2009).
In this paper, we focus on LCMs with binary responses for large-scale data, which are typically collected in modern educational assessments (correct/wrong responses) and psychological or social science surveys (yes/no responses). Such data are characterized by many test takers with large N and also many items with large J. This is in contrary to the traditional regime with large N and fixed J, a relatively well understood setting. Such a large scope of data poses challenges to classical statistical analysis methods and calls for new developments for LCMs. In the following, we summarize the two main questions motivating our study.
The first question of large-scale latent class analysis is how to perform computation efficiently.
A conventional estimation method is the Expectation-Maximization (EM; Dempster et al., 1977) algorithm to maximize the marginal likelihood. EM algorithms have two potential drawbacks, the slow algorithmic convergence rate in high dimensions and a tendency to converge to some local optima when the initial values are poorly chosen (Balakrishnan et al., 2017). In practice, researchers often run EM with many random initializations and select the one that gives the largest log-likelihood value. This procedure can be very time-consuming, especially for large-scale data. It is thus desirable to develop more efficient computational tools for large-scale latent class analysis.
Motivated by this, we introduce the tensor decomposition perspective into latent class analysis.
There has been active research on tensor decompositions since they were introduced in Hitchcock (1927). The concept of tensors appeared in the literature of psychometrics dating back to 1960s-1970s (Tucker, 1964(Tucker, , 1966Harshman, 1970;Kruskal, 1976). The interest of tensor decompositions has expanded to other areas, including chemometrics (Smilde et al., 2005), signal processing (De Lathauwer and De Moor, 1998), and data mining (McCullagh, 2018). In particular, tensor methods have also been used in learning latent variable models (Anandkumar et al., 2012;Hsu and Kakade, 2013;Anandkumar et al., 2014). Anandkumar et al. (2014) summarized the common structures in several different latent variable models and used tensor power method to learn the parameters under a unified framework. Generally speaking, these tensor methods are all based on lower-order moments of observed variables rather than the entire likelihood function. Notably, such tensor decompositions based on the lower-order moments are convex optimization problems, while maximizing the log-likelihood with latent variables is highly non-concave with potentially many local optima. As a result, an advantage of using moment-based tensor decomposition algorithms for learning latent variable models is the provable guarantee of converging to the global optimum of the objective function; see Anandkumar et al. (2014) for more details.
Besides the computational challenge, the second question of large-scale latent class analysis is how to ensure the estimators in the large-N and large-J regime are theoretically valid and meaningful. Traditionally, the subjects' latent class indicators in an LCM are often treated as random variables and marginalized out to obtain the marginal likelihood; we call the resulting model a random-effect LCM. On the other hand, an alternative approach is to treat the subjects' latent class indicators as fixed unknown parameters and directly incorporate them into the likelihood; we call the resulting model a fixed-effect LCM. In the classical scenario with sample size N going to infinity and the number of items J held fixed, the fixed-effect LCMs are known to be inconsistent for estimating the subject-level latent class indicators (e.g., see Neyman and Scott, 1948). However, for data featuring large N and large J, with an increasing amount of information collected per subject, an interesting theoretical question is whether we can obtain consistency in clustering the subjects into their corresponding true latent class in fixed-effect LCMs?
In this paper, in the regime where both N and J go to infinity, we propose an efficient computational pipeline and develop the theory of clustering consistency for LCMs with binary responses. On the computational side, we propose an efficient two-step estimation pipeline integrating the moment-based tensor decomposition method and the EM algorithm. In the first step, we apply the tensor power method for LCMs to quickly and reliably find roughly accurate parameter estimates.
In the second step, we use the tensor estimates as initialization for the EM algorithm to refine the parameter estimation. It is known that the method of moments can be viewed as good complementary to the maximum likelihood approach (Chaganty and Liang, 2013;Zhang et al., 2014;Balakrishnan et al., 2017). With good initialization, EM algorithms typically converge in very few iterations. Therefore, such an estimation pipeline combines the advantages of both the tensor decomposition algorithm and the EM algorithm for latent class analysis. Our extensive simulation studies empirically show that such an estimation pipeline enjoys both computational efficiency and estimation accuracy. Further, on the theoretical side, we prove the clustering consistency of the joint maximum likelihood estimator (joint MLE) for fixed-effect LCMs. That is, we prove that the joint MLE is consistent in estimating the subjects' latent class memberships under certain mild assumptions when N and J both go to infinity. We also derive a bound on the rate of convergence of the joint MLE's clustering performance.
The rest of this paper is organized as follows. The setups of random-effect and fixed-effect LCMs are introduced in Section 2. The proposed estimation procedures of large-scale LCM are presented in Section 3. Some preliminaries about tensor are also provided in Section 3 to make this section self-contained. Section 4 presents our theoretical results on clustering consistency of joint MLE.
Section 5 presents simulation studies that evaluate the proposed estimation procedures and assess the empirical behavior of clustering consistency. A real data example is shown in Section 6, and we conclude this paper with some discussion in Section 7. The proof of the clustering consistency and additional simulation results are presented in Supplementary Material.

Latent Class Models with Binary Responses
In this section we introduce two perspectives of LCM. In random-effect LCMs, the latent class indicators are random variables; while in fixed-effect LCM, the latent class indicators are fixed and treated as unknown parameters. These two models share common assumptions on how the observed variables depend on the latent ones.

Latent Variables as Random Effects
We first introduce random-effect LCM. Consider a binary-outcome latent class model with J items and L classes. There are two types of individual-specific variables in the model, that is a binary response vector R ∈ {0, 1} J and a latent variable z ∈ [L]. Here [L] = {1, 2, . . . , L} is the set of positive integers smaller than or equal to L. The response vector R i = (R i,1 , . . . , R i,J ) contains the observed responses to the J items of ith subject. The jth component of R i will be 1 if this subject gives a positive response to the jth item and will be 0 otherwise. For instance, in a test with J items, if a student answers the jth item correctly, then R i,j , the jth component of R i , will be 1.
If the student fails to give a right answer then R i,j = 0. The latent variable z i is introduced to categorize different observations and explain the dependence among items.
The generative process for a response vector R i of an observation is as follows: first the class of this observation z i is drawn from a discrete distribution specified by the probability vector p = (p 1 , p 2 , . . . , p L ), where p k ≥ 0 and L k=1 p k = 1. So we have where p is the proportion of subjects belonging to th class in the population. Then given the latent class z i = , the responses to J items are drawn conditionally independently from a Bernoulli distribution with parameter θ j, for each item j. That is So θ j, measures the ability of subjects from th class to give a positive response on item j and is also known as item parameters. Like many other latent variables, local independence is assumed here, implying the dependence of item responses is fully explained by the latent classes. We collect all the item parameters for the L classes in the matrix θ = (θ j, ) ∈ [0, 1] J×L whose rows are indexed by the J items and columns indexed by the L classes. All the response vectors are collected in a N × J matrix R, and the corresponding log-likelihood function under the random-effect LCM is with (p, θ) the parameters to be estimated.

Latent Variables as Fixed Effects
Another way to model latent classes is to view latent class assignment as fixed unknown parameters.
For a fixed-effect LCM, denote the ith subject's latent class membership by a vector of binary entries Z i,· = (Z i,1 , . . . , Z i,L ), with Z i, = 1 if subject i belongs to the latent class . We also introduce another notation for the latent class membership z i ∈ {1, 2, . . . , L} and z i = corresponds to Z i,l = 1. Given a sample of size N , collect all the Z 1,· , . . . , Z N,· in a N × L matrix Z, then each row of Z contains only one entry of "1" and the remaining entries are zeros. We will use the two equivalent notations Z and z = {z 1 , . . . , z N } interchangeably. The components of response vector R i are independent Bernoulli variables with parameters specified by θ. So we have P (R i,j = 1) = θ j,z i .
The log-likelihood for (Z, θ) takes the following form The parameters to estimate are (Z, θ). The above display is also called the complete data likelihood in the literature.
In fixed-effect LCM, if we imagine that z is sampled from a discrete distribution which is separated from the sampling process of R, the data generating process of fixed-effect LCM and random-effect LCM coincides. So the estimation procedures developed in the following sections for random-effect LCM also apply to fixed-effect LCM.

Estimation Procedures
The EM algorithm is a popular method to maximize likelihood and estimate parameters in LCM by iterating between E-step and M-step. In E-step the probability of each subject belonging to each class is updated by current estimates of item parameters, and in M-step item parameters are updated given the probabilities of each subject's latent class membership. However, the likelihood function under LCM is nonconcave due to the mixture model formulation. Hence, EM algorithm may suffer from convergence to local optima and slow convergence rate under poor initializations.
Good initial values are critical to the success of the EM algorithm. In this section we introduce tensor method in Anandkumar et al. (2014) to find good initializations and hence improve the performance of EM algorithm. We will follow the discussions of Anandkumar et al. (2014) and be succinct and self-contained.

Preliminaries about Tensor
First we introduce some notations borrowed from Anandkumar et al. (2014). A real p-th order tensor T ∈ ⊗ p i=1 R n i is a p-way array of real numbers where [T ] i 1 ,...,ip is the (i 1 , . . . , i p )-th entry in the array.
We will mostly consider the case where n i = n for all i ∈ [p]. Vectors and matrices are special cases of tensors where p = 1 and p = 2, respectively. Another view of tensor is that it is a multilinear map from a set of matrices In this paper we will mainly consider three-way tensor and third-order case of this multilinear map. For a third-order tensor T ∈ ⊗ 3 R d and a vector u ∈ R d , we will make use of the following vector-valued map in the iteration of tensor power methods where e 1 , . . . , e d are the canonical basis vectors of R d ; that is, each e k is a d-dimensional vector with only the kth entry being one and the other entries being zero. And we will also use the following map in the iteration These maps are all special cases of (2).
Most tensors we consider in this paper are symmetric tensors, which means that an element of a tensor is invariant to permutations of its coordinates. If T ∈ ⊗ p R d is a symmetric tensor, then we have [T ] i 1 ,...,ip = [T ] i π(1) ,...,i π(p) for all permutations π on [p]. This concept is a generalization of symmetric matrices.
A simple case of a tensor is called rank-one tensor. A rank-one tensor T ∈ ⊗ p R d can be expressed as tensor product of p vectors:  Anandkumar et al. (2014) showed that for some latent variable models, their low-order moments can be expressed as a sum of rank-one tensors. Once this structure of cross moments is obtained for a particular model, one can apply the orthogonal tensor decomposition to learn the parameters of the model. For random-effect LCM, we show that there is also useful tensor structure in low-order moments by examining Theorem 3.6 in Anandkumar et al. (2014), which studies the multi-view models.

Tensor Structure in Random-Effect LCM
Recall that we use R to generally denote a response vector of length J. We consider R on a population level and divide the items into three disjoint parts (so we assume J ≥ 3) R 1 , R 2 , R 3 with each R t ∈ R Jt and J 1 + J 2 + J 3 = J. The item paramters of R t are denoted by θ t ∈ R Jt×L , which is a sub-matrix of θ, with rows corresponding to rows in R t . We need the following assumption to derive the tensor structure.
Note that the partition of items can be arbitrary as long as the item parameters for each part satisfy Condition 1. So we can try different partitions to estimate the parameters and take average to obtain the final estimates.
We denote the ith column of θ t to be θ t,i . The following theorem restates Theorem 3.6 in Anandkumar et al. (2014) in our setting and characterizes the tensor structure in random-effect LCM.
Theorem 1. Assume that Condition 1 holds and p 1 , . . . , p L > 0. Define where A + denotes the Moore-Penrose pseudoinverse of matrix A. Then we have Proof. First we compute the cross moment. For t = t , R t and R t are conditionally independent, we have If we denote D = diag{p 1 , . . . , p L }, then we have E[R t ⊗ R t ] = θ t Dθ T t . Then we calculate the conditional mean According to the model setting E[R 2 |z = i] = θ 2 e i , then we have Similarly one can get the decomposition for M 3 in (6).
In applications we only have finite samples and the moments in Theorem 1 should be approximated by empirical moments. In particular, once we have samples R 1 , . . . , R N ∈ R J , we partition each sample and obtain R t i ∈ R Jt corresponding to the partition on population level. Then the transformed response and estimated moments can be computed by Due to the randomness of sample, it is possible that r has (r − L) extra non-zero singular values. These singular values will be small since they equal to 0 . In this case we should discard these extra singular values and only use first L sigular values when calculating ( E[R 2 ⊗ R 3 ]) + , otherwise one has to compute the inverse of these small singular values which will incur large error.
After learning θ 1 from data by the tensor power method to be introduced in Section 3.3, we can obtain θ 2 and θ 3 by setting This can be derived in a same way as Theorem 1. So the main problem is to estimate θ 1 from moments M 2 and M 3 .
Although this structure only holds for random-effect LCMs, as we mentioned above, in fixedeffect LCMs we can view z 1 , . . . , z N as random. For instance, they are sampled from some discrete distributions on [L] independently from R. Then the data generation process of random-effect and fixed-effect LCMs are the same and the estimation procedures for random-effect LCMs also apply to fixed-effect LCMs.

Tensor Method to Learn the Parameters
In this section we briefly describe the procedures in Anandkumar et al. (2014) to recover the parameters in (6). That is, given where µ i ∈ R d , we want to obtain the elements of decomposition (w i , µ i )'s from M 2 and M 3 .
Condition 1 now becomes {µ 1 , . . . , µ L } are linearly independent. First we introduce the orthogonal decomposition of a tensor. Then we see how we can use tensor power method to recover the orthogonal decomposition of a tensor and estimate the parameters.

Orthogonal Decomposition
Since the moments structures in (8) are about at most a third-order tensor, we only consider the case p = 3 (third-order tensor).
A symmetric tensor T ∈ ⊗ 3 R d has an orthogonal decomposition if there exists a collection of orthonormal unit vectors {v 1 , . . . v L } and positive scalers λ i > 0 such that Without loss of generality we assume λ i > 0 because for third-order tensor we have −λ . This is a generalization of spectral decomposition for a symmetric matrix. We can also generalize the concept of eigenvalue and eigenvectors.
Recall the definition of the multilinear map induced by a tensor in (2). A unit vector u ∈ R d is an eigenvector of T with corresponding eigenvalue λ ∈ R if T (I, u, u) = λu. For an orthogonally and operation (4) reduces to

By the orthogonality of
The eigenpairs of a tensor are more complicated than those of a matrix and there are some subtle points; see Anandkumar et al. (2014) for more discussion. Fortunately, Theorem 4.1 and 4.2 in Anandkumar et al. (2014) guarantee the uniqueness and identifiability of orthogonal decomposition in (9).

Whitening Process
Comparing the third-order moment structure M 3 in (8) with the orthogonal decomposition form (9), we find they have almost the same form except that the vectors µ i 's in (8) may not necessarily be orthogonal to each other. So we need to whiten the tensor M 3 to M 3 , which has an orthogonal decomposition. In the whitening process we will make use of M 2 in (8).
is the diagnoal matrix containing all positive eigenvalues of M 2 and U ∈ R d×L is the matrix of corresponding orthogonal eigenvectors of M 2 . D −1/2 is well-defined since we assume µ i 's are linearly independent and thus M 2 is of rank L.
So µ i s are orthonormal vectors. Define Since µ i s are orthonormal vectors, this is the orthogonal decomposition of M 3 . We can use tensor power method described in section 3.3.3 to obtain the eigenpairs (λ

Tensor Power Method
Now we show how to recover the parameters (λ i , v i )'s in (9) from a tensor T . In analogy to matrix power method, here we use the tensor power method of De et al. (2000) to obtain the eigenpairs (9). First suppose a third-order tensor has an exact orthogonal decomposition. Then we have the following result on the algorithmic convergence of tensor power method (Lemma 5.1 in Anandkumar et al. (2014)).
Theorem 2. Let T ∈ ⊗ 3 R d have an orthogonal decomposition as given in (9). For a vector θ 0 ∈ R d , suppose that the set of numbers So repeated iteration starting from θ 0 converges to v 1 . To obtain all the eigenpairs, we use "deflation" after getting an eigenpair (λ i , v i ). That is, to obtain the j-th eigenpair, we subtract the previous j − 1 rank-one structures from T and then execute the power method on However, in practice we plug-in the empirical estimate of M 2 and M 3 in (6) and the estimate of whitened tensor M 3 may not have an exact orthogonal decomposition. So we should consider the case where we only have an approximation T of T and need a more robust algorithm to use an orthogonal decomposition to approximate T . Following Anandkumar et al. (2014), we present the following Algorithm 1. Multiple starting points are used in Algorithm 1 to ensure approximate convergence at first stage. Intuitively, by restarting from different points we can start from a point from which the initial n iteration from it dominates the error T − T . Perturbation analysis (Theorem 5.1 in Anandkumar et al. (2014)) shows that we can find such a point with high probability with K = poly(L) trials. In our simulation settings we find tensor estimator is good enough as initializations for EM algorithm with K = 10 trials on starting points.
In conclusion, the procedures to learn (ω i , µ i ) from M 2 and M 3 in (8) are : (1) Use the information of M 2 to whiten M 3 and get M 3 ; (2) Apply tensor power method to M 3 and learn the orthogonal

Algorithm 1 Robust tensor power method
Input: symmetric tensor T ∈ R d×d×d , number of iterations K, n.
Output: estimates of one of eigenpairs; the deflated tensor 1: for τ = 1 to K do 2: Draw θ (τ ) 0 uniformly from unit sphere in R d .

3:
for t = 1 to n do 4: Compute power iteration and re-normalization to obtain θ, and set λ = T ( θ, θ, θ). 9: return the estimated eigenpair ( λ, θ); the deflated tensor T − λ θ ⊗ θ ⊗ θ. In application we can plug-in the estimates of M 2 and M 3 and learn the parameters from samples.
In our random-effect LCM, we can apply this procedure to (6) and learn p i 's and θ. When the sample size is large enough, this method alone can yield an estimator close to the true parameters.
However if we do not have so many samples, estimates based on tensor are not so accurate since in tensor method we only take advantage of low-order moments and ignore other information of sampling distributions. So after obtaining the tensor estimates, we use them as initial values for EM algorithms to improve the accuracy. For fixed-effect LCM, we first learn p i 's and θ from (6) as a initialization step and apply Classification-EM (CEM) algorithm proposed in Celeux and Govaert (1992) to obtain the final estimator. We call this two-step estimation procedure the tensor-EM method.
Empirically, we find that the accuracy of this method is comparable to estimates obtained with EM algorithm starting from true parameters, indicating that the tensor-EM method can give the MLE of latent class model when the model is correct. Moreover, it is more computationally efficient than EM algorithm with several random initial values, especially in large-scale data. See our simulation study for more details.

Selecting the Number of Classes
In the discussion above we assume the number of classes L is known. Next we discuss the model selection of L. There exists a rich literature in selecting number of classes. Nylund et al. (2007) performed a Monte Carlo simulation study on several commonly used methods and found BIC proposed in Schwarz (1978) and likelihood ratio test based on bootstrap in McLachlan and Peel (2004) have a better performance. They recommend BIC and likelihood ratio test based on bootstrap to select the number of classes. Since we focus on large-scale datasets containing many items and samples, we follow the discussion of Chen et al. (2017) and apply generalized information criterion proposed in Nishii (1984) to selecting the number of classes.
Specifically, for a candidate set L and any L ∈ L , we apply the tensor-EM algorithm to learning the parameters and compute the generalized information criterion for random-effect and fixed-effect LCM as follows: where dim(p, θ) is the dimension of parameters to estimate and measures the model complexity.
We have dim(p, θ) = JL + L − 1 in random-effect LCM and dim(Z, θ) = JL + N in fixed-effect LCM. Sample size dependent quantity a N measures the level of penalty on model complexity. Here we consider two choice of a N .
• GIC 1 : a N = log(N ). This case corresponds to BIC and it enjoys some consistent results shown in Nishii (1984) when the model has low complexity (i.e. the dimension of parameters is fixed).
• GIC 2 : a N = log[log(N )]logN . This choice is considered in Fan and Tang (2013)  can be large and should not be treated as fixed. For instance, in the simulation study a random-effect LCM we consider has ten classes and one hundred items. This model has dimension d = 1009 while the sample size is N = 1000. So it is appropriate to adopt this choice in this setting. See Fan and Tang (2013) for discussion about theoretical results of this choice of a N .
After calculating the GIC for different models, we select the number of classes to minimize GIC(L): Although our simulation results show that the tensor-EM method gives the MLE of the model when the number of classes L is correctly selected, it may converge to some local optima when L is incorrect. It is also likely that the tensor method yields inaccurate results when L is misspecified as L. For instance, some item parameters may be negative or over one, which happens when the sample size is small. In this case, we propose to revise the tensor estimates as follows: we set the negative components to be a small positive number (e.g. 0.001) and those over one to be a number close to 1 (e.g. 0.999). Such a procedure will help us select the number of classes because we know as long as the true number of classes is specified, the tensor-EM method can yield an ideal estimate (MLE) with a small GIC value. On the other hand, a poor estimate based on a misspecified model will give a larger GIC value, indicating misfit of the model. Hence, GIC values computed by tensor-EM method can provide useful information to select the number of classes. According to our simulation study, the proposed method can select the right model most of the time.
After proposing the computational methods to select number of classes and to find the MLE, we next examine the theoretical properties of MLE in the large-N and large-J scenario. For random-effect LCM with fixed number of items J, the MLE is known to be consistent. However, the joint MLE for fixed-effect LCM may not be consistent (Neyman and Scott, 1948) when J is fixed. Intuitively, one cannot hope to recover each subject's latent class membership accurately with only a finite number of items observed for each subject. So in the next section, we will consider the consistency of joint MLE when J also goes to infinity in fixed-effect LCMs.

Clustering Consistency of the Joint MLE
In this section we consider large-scale fixed-effect LCMs and characterize the behavior of latent class assignment estimator Z under suitable conditions, where ( θ, Z) is the joint maximum likelihood estimator (MLE). We use a similar proof technique as in Gu and Xu (2020) to establish the clustering consistency of the joint MLE for fixed-effect LCMs.
First we need to define some notations. Denote the true parameters by (θ 0 , Z 0 ). Define The M ∈ [0, 1] above measures the average positive response rate over all subjects and items. Denote where the expectation is taken with respect to the distribution of R.
Given arbitrary Z, denote ), Then under any realization of Z, the following holds for any latent class a ∈ [L], We consider the joint maximum likelihood estimator ( θ, Z) subject to fitting a L-class fixed-effect maximizes the profile likelihood (R; Z, θ) given a particular realization Z. In Section 3 we will discuss procedures to compute the joint MLE efficiently.
We impose the following assumptions on the true parameters.
Assumption 1. There exists a constant d > 0 such that Assumption 2. There exists a positive sequence {β J } such that where · denotes the 2 norm.
Assumption 1 guarantees that the components of θ are bounded away from 0 and 1 but allowed to become very close to 0 or 1 as J becomes larger. It is a quite mild technical assumption. Assumption 2 is an identification condition for latent classes and guarantees that the item parameters of different class are different enough. Note that we allow different classes to have same probability to answer a single item correctly (i.e. θ j,a = θ j,b for some a = b). But their average performance on the J items should be different.
In fitting the latent class model, we are interested in controlling the number of incorrect latent class assignments. Formally, after obtaining some estimator z, let C 1 , . . . , C L be clusters from our estimator z such that subjects sharing same estimated membership are in one cluster. Let m l = arg max a∈[L] i∈ C l Z 0 i,a be the majority of true class membership among subjects in C l . Since the latent classes are identified up to permutations of class index, m l should be viewed as the class index of C l corresponding to true class assignments z 0 . The number of incorrect class assignments under estimated latent class assignments N e ( z) is defined as: So every subject i ∈ {1, . . . , N } whose true class under z 0 is not in the majority within its estimated class under z is counted. For instance, if we have eight subjects from two classes and we end up with estimated clusters C 1 = {1, 1, 1, 1, 2} and C 2 = {1, 2, 2}, where "1" and "2" represents the true class index for each subject. Then m 1 = 1, m 2 = 2 and number of errors N e ( z) = 2.
We have the following main theorem on the clustering consistency of joint MLE for fixed-effect LCM, which characterizes the asymptotic behavior of error rate N ( z)/N .
Assigning each subject to latent class resembles the process of clustering and hence we name our results as "clustering consistency". In  We compare the performance of the proposed tensor-EM method with three other methods: (1) EM-true, which is the EM algorithm starting from the true parameters as initial values; (2) EM-random, where we randomly generate the initial values for the EM algorithm and repeat five times, and then we select the estimators corresponding to maximum log-likelihood; (3) the tensor power method alone. In this third competitor, we permute the items randomly and obtain a tensor estimate for each permutation. Let π be a permutation of [J] and R π be subject response vector corresponding to permutation π (i.e. [R π ] i = R π(i) ). We obtain tensor estimates θ π from cross moments of R π and set the tensor estimates of original item parameters as θ j,a = θ π π −1 (j),a . We then repeat this procedure five times and finally take average of them. This repetition can reduce the MSE of item parameters a little but will remain the same magnitude.
We emphasize that in the proposed tensor-EM method, we do not repeat the tensor power method or take any average. Empirically, just one implementation of the tensor power method gives good initial values for the EM algorithm.

Verification of Clustering Consistency
In this subsection, we empirically verify Theorem 3 in fixed-effect LCMs with diverging N and J.
Specifically, we consider fixed-effect LCM with L = 5 classes. We let J increase from 20 to 100 by 10 and set N = 10J in all the simulations. Item parameters are generated uniformly from either

Performance of GIC in Selecting Number of Classes
We consider the accuracy of GIC in selecting L, which needs to be estimated in practice. The settings we consider are the same as in Section 5.1. The performance of GIC to select the number of classes in the random-effect LCM and the fixed-effect LCM are reported in Table 1. For settings with the true number of classes L being five, we let the candidate set of L to choose from be {2, 3, 4, 5, 6, 7}; while for settings with the true L being ten, we let the candidate set of L be {7, 8, 9, 10, 11, 12}. We see that for all these settings GIC 1 can always select the correct number of classes. And for most of the settings GIC 2 can choose the right model. The only setting where GIC 2 performs not so well is random-effect LCM with N = 1000, J = 100, L = 10, θ j,a ∈ {0.2, 0.4, 0.6, 0.8}. In general, both GIC 1 and GIC 2 enjoy desirable performance in selecting the correct number of classes.

Real Data Analysis
In this section we apply the proposed method to IPIP HEXACO equivalent scales data. The raw data can be found on https://openpsychometrics.org/_rawdata/. in total and 40 items measuring each one of them. Since all items are rated on a seven point scale, we convert them to binary in the following way. For those items positively related to a personality (marked by a '+' on the website), 5 (slightly agree), 6 (agree) and 7 (strongly agree) will be regarded as positive response and the answers will be set to 1. Otherwise they are set to 0. For those items negatively correlated with a quality (marked by a '-' on the website), 1 (strongly disagree), 2 (disagree) and 3 (slightly disagree) will be regarded as positive response and the answers will be set to 1. Otherwise they are set to 0.
After transforming the data, we consider fitting a random-effect LCM for the data. We apply the GIC method to selecting the number of latent classes. According to Figure 8, the GIC criterion selects eight classes for the data. Note that the real data may not strictly follow the model setup and the tensor-EM estimates may vary, depending on the initial values of the robust tensor power method in Algorithm 1. We further run the proposed method multiple times and find that fitting a model with fewer classes yields a higher GIC than the model fitted with eight classes shown in the p being fairly small (smaller than 0.001); thus the corresponding latent class can be neglected for better interpretability and parsimony, and it is still reasonable to select eight classes. We consider the fitted eight-class model in Figure 8. After obtaining the item paramters θ ∈ R 240×8 , we pool the 40 item parameters corresponding to each of six personalities for each class. To be more concrete, for each l-th class and k-th personality, we calculate the mean of item parameters in the set { θ j,l , 40(k − 1) + 1 ≤ j ≤ 40k} and obtainθ k,l = 40k j=40(k−1)+1 θ j,l /40.
These personality-level item parameters are collected in a matrixθ ∈ R 6×8 . Then each entrȳ θ k,l is compared with certain threshold to decide whether it shows a high level possession of k-th personality for l-th class. The results are shown in Figure 9. The left panel is based on the absolute value and the right based on quantiles in θ. In the left panel we compareθ k,l with 0.4, 0.6 to decide the level of k-th personality of l-th class. Ifθ k,l < 0.4, it is of low level (marked with blue). If 0.4 ≤θ k,l < 0.6, then it is of medium level (marked with green). Otherwise it is said to be high level (marked with yellow). In the right panel we use some data-dependent threshold. We define and ifθ k,l < θ k,0.4 we think the k-th personality of l-th class is of low level (blue). If θ k,0.4 ≤θ k,l < θ k,0.6 , it is said to be of medium level (green). Otherwise it is said to be of high level (

Discussion
This paper investigates the computation and theory for large-scale latent class models. In terms of computation, commonly used likelihood-based methods (e.g. EM algorithm) suffer from slow convergence rate and potential convergence to local optima under poor initializations. Recent developments in tensor decomposition and its applications provide a computationally efficient moment-based method to estimate the parameters. However, such tensor method is based on low-order moments of the observed variables rather than the entire likelihood function. Hence it is not statistically efficient and generally requires a large number of samples to ensure the accuracy of estimates. In this work, we propose a two-stage tensor-EM estimation pipeline which combines these two methods. Simulation studies empirically show the proposed procedure is both computationally efficient and statistically accurate. Further, we also theoretically establish the clustering consistency We also note that it is also possible to find similar tensor structures in more complicated latent structure models, for example, the diagnostic classification models or cognitive diagnostic models (CDMs) (Rupp and Templin, 2008;von Davier and Lee, 2019). This is because the low-order moments constructed in and exploited by the tensor power method here (also see Anandkumar et al., 2014) are essentially resulting from the local independence assumption in LCMs; that is, the observed variables are conditionally independent given the latent variable. Since CDMs can be viewed as extensions of LCMs where the local independence assumption also holds, the same tensor power method used here may also be used to find rough estimates of parameters in a CDM. However, CDMs also have a unique feature, the Q-matrix (Tatsuoka, 1983), which induces complicated parameter Zhang, Y., Chen, X., Zhou, D., and Jordan, M. I. (2014). Spectral methods meet EM: A provably optimal algorithm for crowdsourcing. arXiv preprint arXiv:1406.3824.

Appendices
In this appendices, the proof of Theorem 3 (clustering consistency) is presented in Appendix 1. In Appendix 2, more simulation results are provided to show the good performance of the proposed EM-tensor method.
Appendix 1: Proof of Theorem 3 Outline of proof idea. The proof follows the following 8 steps.
Step 1: Express (R; Z) −¯ (Z) in terms of a n a j D( θ j,a θ j,a ) + X − E(X), where X is a random variable depending on R andθ (Z) under Z, and n a = N i=1 Z i,a .
Step 2: Bound the first term j a n j,a D( θ j,a θ j,a ) in the above display uniformly over all possible Z.
Step 3: Bound the second term X − E(X) using possibly Bernstein type inequality. Combine this and Step 2 to obtain a bound for sup Z | (R; Z) −¯ (Z)|.
Step 4: (Denote the true latent class memberships by Z 0 and the those maximizing the likelihood by Z.) Establish¯ (Z 0 ) ≥¯ (Z) for all Z. Use triangle inequality to upper-bound the non-negative Since in the above display the middle group of terms [ (R; Z 0 ) − (R; Z)] ≤ 0, we have 0 ≤ Step 5: Introduce the notion of partitions and generalize¯ (Z) to¯ (Π).
Step 8: Apply results in step 6 and step 7 to MLE z, we have Now we formally begin the proof of Theorem 3 in the above several steps. In derivations below we abbreviateθ (Z) asθ and θ (Z) as θ to simplify notations.
Step 1. Define D(p q) = p log(p/q)+(1−p) log((1−p)/(1−q)), the Kullback-Leibler divergence of a Bernoulli distribution with parameter p from that with parameter q. In this step we prove a lemma as follows.
Lemma 1. Let (R i,j ; 1 ≤ N, 1 ≤ J) denote independent Bernoulli trials with parameters (P i,j ; 1 ≤ N, 1 ≤ J). Under a general latent class model, given an arbitrary Z, there is is random variable depending on Z and Given a fixed Z, denote n (Z) The maximizing properties of θ j,a andθ j,a in 16 imply that n a θ j,a = i Z i,a R i,j , n aθj,a = i Z i,a P i,j .
Using (17), we have the following, Define the random variable then X depends on Z and the above display becomes the summation of L a=1 n a j D( θ j,a θ j,a ) and X − E[X]. This establishes (15) in Lemma 1. In the following, we bound the first term L a=1 n a j D( θ j,a θ j,a ) and the second term X − E[X] in the above display uniformly over all possible Z, respectively in Step 2 and Step 3.
Step 2. In this step we prove the following theorem.
Theorem 4. The following event happens with probability at least 1 − δ, Given any fixed latent class memberships Z, every θ j,a is an average of n a independent Bernoulli random variables R 1,j , . . . , R N,j with meanθ j,a . We apply the Chernoff-Hoeffding theorem to obtain P( θ j,a ≥θ j,a + t) ≤ e −naD(θ j,a +t θ j,a ) , P( θ j,a ≤θ j,a − t) ≤ e −naD(θ j,a −t θ j,a ) .
Note that given a fixed Z, each θ j,a can take values only in the finite set {0, 1/n a , 2/n a , . . . , 1} of cardinality n a + 1. We denote this range of θ j,a by Θ j,a . Then P( θ j,a = ϑ) ≤ exp{−n a D(ϑ θ j,a )} for any ϑ ∈ Θ j,a . Now consider the cardinality of the set Θ given Z. Since for each of the J × L entries in θ, θ j,a can independently take on n a + 1 different values, there is | Θ| = [ a (n a + 1)] J .
Considering the natural constraint L a=1 n a = N , we have (n a + 1) 36 Define Θ = { θ ∈ Θ : j a n a D( θ j,a θ j,a ) ≥ }, then Θ ⊆ Θ. Note that the components of θ depend on different components of {R i,j , i ∈ [N ], j ∈ [J]} and thus are independent. We have The above result holds for fixed Z, so applying a union bound over all the L N possible assignment Now take δ = L N N L + 1 JL e − , then = N log L + JL log( N L + 1) − log δ. Therefore the following event happens with probability at least 1 − δ, This concludes the proof of Theorem 4.
Step 3. In this step we bound |X−E[X]|, with X defined in (18).
ing the Bernstein's inequality to the sum of independent bounded random variables, we have the following holds for any fixed Z, .
We next prove the following theorem.
Theorem 5. Under the following scaling, we have for any η > 0.
In this step we bound | (R; Z) −¯ (Z)| uniformly over all the Z. Combining the results of Step 2 and Step 3, since that there are L N possible assignments of Z, we apply the union bound to obtain In order for the second term on the right hand side of the above display to go to zero, the following of δ N J would suffice, for a small positive constant η. Then the right hand side of (22) goes to zero as N, J go large and hence the scaling of J described in the theorem yields This proves Theorem 5.
Step 4. Denote the true class assignments by Z 0 . We first establish First note that θ 0 j,z 0 i = P i,j , and The difference¯ (Z 0 ) −¯ (Z) can be written as therefore establishing (25). Since the above holds for every Z, it also holds for the maximum likelihood estimator Z. We further upper bound¯ (Z 0 ) −¯ (Z) from above as follows, Step 5.
as the average over all i in the subset S a indexed by Π i = a. We then define generalization of¯ (Z) Note thatθ Π Z j,a =θ Z j,a and hence¯ (Π Z ) =¯ (Z) when the partition Π Z is induced by latent class assignment Z.
We will proceed as follows: in step 6 we show a refined partition increases¯ (·). We then construct a refined partition Π * for every partition Π Z induced by Z and prove¯ (Z 0 ) −¯ (Π * ) ≥ 1 2 N e (z)β J in step 7. Finally we apply the results to MLE Z and obtain the desired results in step 8.
The first equality is obtained by rewriting¯ (Π) in terms of subsets. The inequality follows from non-negativity of K-L distance. Then we combine terms in same class under Π and obtain the second equality.
Step 7. Now we prove a result on refinement.
Lemma 3. For any latent class assignment Z, there exists a partition Π * that refines Π Z and For a given Z, partition each latent class assigned by Z into sub-classes according to true assignments Z 0 of each sample. For each sample i 1 that is incorrectly assigned by Z (by definition this means its true class under Z 0 is not in the majority within its estimated class under Z), we find another sample i 2 assigned to same class under Z but i 1 and i 2 belong to different class under Z 0 and make these two samples (i 1 , i 2 ) a pair. We allow two misclassified samples to form a pair.
Note that since incorrectly assigned samples are not in the majority of that class, we can find a pair for each of them.
The refined partition Π * contains all such pairs and remaining correctly assigned samples in all classes assigned by Z. So for the above example, the refined subset for that class is . Let e(z) by the set of incorrectly assigned sample. Clearly Π * is a refinement and we have¯ where i and i are in different classes under Z 0 while in same subset under Π * by definition. Apply Pinsker's inequality we have Hence we have¯ Step 8. Apply Lemma 3 to MLE z, there exists a refinement of Π Z denoted as Π * such that which completes the proof.

Appendix 2: More simulation results
First we present simulation results for random-effect LCM.