Mixed Membership sub-Gaussian Models Huan Qinga,∗ a School of Economics and Finance, Chongqing University of Technology, Chongqing, 400054, China
arXiv:2604.22633v1 [stat.ML] 24 Apr 2026
Abstract The Gaussian mixture model is widely used in unsupervised learning, owing to its simplicity and interpretability. However, a fundamental limitation of the classical Gaussian mixture model is that it forces each observation to belong to exactly one component. In many practical applications, such as genetics, social network analysis, and text mining, an observation may naturally belong to multiple components or exhibit partial membership in several latent components. To overcome this limitation, we propose the mixed membership sub-Gaussian model, which extends the classical Gaussian mixture framework by allowing each observation to belong to multiple components. This model inherits the interpretability of the classical Gaussian mixture model while offering greater flexibility for capturing complex overlapping structures. We develop an efficient spectral algorithm to estimate the mixed membership of each individual observation, and under mild separation conditions on the component centres, we prove that the estimation error of the per-individual membership vector can be made arbitrarily small with high probability. To our knowledge, this is the first work to provide a computationally efficient estimator with such a vanishing-error guarantee for a mixed-membership extension of the Gaussian mixture model. Extensive experimental studies demonstrate that our method outperforms existing approaches that ignore mixed memberships. Keywords: Mixed membership, sub-Gaussian mixture model, spectral method
1. Introduction Clustering is a fundamental task in unsupervised learning. Its goal is to partition unlabelled data into meaningful groups, thereby revealing hidden structures that are not directly observable (Ng et al., 2001; Von Luxburg, 2007). The importance of clustering extends to many areas of modern data analysis. In genomics, clustering helps identify cancer subtypes from gene expression profiles (Jain, 2010). In marketing, it segments customers based on purchasing behaviour. In social network analysis, it uncovers communities of individuals with shared interests (Fortunato & Hric, 2016). In medical imaging, clustering aids in disease diagnosis and tissue segmentation. In recommender systems, ∗ Corresponding author.
Email address: [email protected] Preprint submitted to
&
[email protected] (Huan Qing) April 27, 2026
it groups users with similar preferences to improve personalised suggestions. In satellite image analysis, it classifies pixels into land cover categories. In text mining, it groups documents by topic. In speech processing, it separates speakers. These diverse applications show why clustering has remained a central focus in statistics and machine learning for decades. Among the many clustering models, the Gaussian mixture model (GMM) is a cornerstone of unsupervised learning. Its core assumption is that each observation is generated from one of several latent Gaussian components, with the component membership unknown. This simple yet powerful framework has a long history, dating back to the work of (Pearson, 1894) on mixture distributions. The GMM’s flexibility and interpretability have made it a standard tool for density estimation and cluster analysis. Classical estimation methods for GMMs include the expectation-maximization (EM) algorithm (Dempster et al., 1977) and the k-means algorithm (Lloyd, 1982; McQueen, 1967). In fact, k-means can be obtained as a limiting case of EM when the component covariances are taken to be equal and isotropic, and the common variance tends to zero. Building on these foundations, a rich body of literature has developed sharp statistical guarantees for clustering under GMMs. (Löffler et al., 2021) proved that spectral clustering achieves the optimal mislabelling rate, matching the information-theoretic lower bound. (Ndaoud, 2022) gave a complete characterization of the exact recovery threshold for the two-component Gaussian mixture, revealing a sharp phase transition. (Chen & Yang, 2021) extended this analysis to multiple components, pinpointing the critical signal-to-noise ratio for exact recovery. (Li, 2025) further studied the exact recovery problem for k-component Gaussian mixtures, providing both necessary and sufficient conditions. Algorithmic aspects have also been rigorously investigated. (Lu & Zhou, 2016) established statistical and computational guarantees for Lloyd’s algorithm, showing that with a suitable initialisation it attains the optimal rate. (Balakrishnan et al., 2017) analysed the EM algorithm from a population-to-sample-based perspective, elucidating its convergence behaviour. (Srivastava et al., 2023) proposed a robust spectral clustering algorithm for sub-Gaussian mixture models in the presence of outliers. (Chen & Zhang, 2024) achieved optimal clustering in Gaussian mixtures with anisotropic covariance structures. (Jana et al., 2025) established optimality guarantees for adversarially robust clustering. Together, these works have profoundly advanced our understanding of when and how GMMs can be consistently estimated. A fundamental limitation of the classical Gaussian mixture model is its single-membership assumption: each observation is assumed to arise from exactly one latent component. In many real-world applications, however, this assumption is overly restrictive. Consider a social network, where an individual may simultaneously belong to several communities, such as a family, a circle of colleagues, and a sports club. A document can naturally cover multiple topics, like politics and economics. A gene may participate in several biological pathways. In all these cases, a binary 2
assignment to a single group is inadequate. A more adequate representation should allow each observation to belong to multiple components at the same time. The classical GMM and the rich body of work built upon it, including the expectation-maximization algorithm (Dempster et al., 1977) and sharp recovery results (Ndaoud, 2022; Chen & Yang, 2021; Li, 2025), are all designed for the single-membership setting. Consequently, they cannot be used to estimate fractional memberships or to detect overlapping clusters. This limitation has long been recognised in network analysis. There, the mixed membership stochastic blockmodel (MMSB) (Airoldi et al., 2008) generalises the classic stochastic blockmodel (SBM) (Holland et al., 1983) by allowing each node to belong to multiple communities with fractional intensities. This seminal work has inspired a rich line of research on computationally efficient and provably consistent estimation methods. For instance, (Mao et al., 2021) developed an efficient spectral algorithm, which exploits the simplex geometry of eigenvectors and provides sharp row-wise deviation bounds. (Jin et al., 2024) proposed Mixed-SCORE, a spectral method that uses a carefully designed ratio matrix to remove degree heterogeneity. Motivated by these advances, this paper aims to construct a similarly principled and computationally efficient mixed membership extension for the Gaussian mixture model, adapted to continuous data. This paper introduces the mixed membership sub-Gaussian model (MMSG) to overcome the single-membership limitation of classical Gaussian mixture models. A formal definition of MMSG is given in Section 2. The MMSG extends the classical Gaussian mixture model by allowing each observation to belong to multiple components with fractional weights, and the noise is sub-Gaussian. The main contributions of this work are as follows. • We introduce the MMSG model, which naturally captures overlapping clusters and fractional memberships, making it suitable for a wide range of applications where observations can belong to multiple latent groups. • Under the MMSG model, we develop an efficient spectral estimator for the mixed membership vectors. The method is simple to implement, does not require iterative optimization, and works in high-dimensional settings where the number of features can be much larger than the number of samples. • We prove that, under mild separation conditions on the component centres, the individual-wise estimation error for each individual’s membership vector vanishes with high probability. Our analysis accommodates highdimensional settings and does not require the clusters to be balanced in size. • Extensive experimental studies demonstrate that our method substantially outperforms existing algorithms designed for classical Gaussian mixture models. The remainder of the paper is organized as follows. Section 2 introduces the model. Section 3 presents the 3
algorithm. Section 4 gives the theoretical guarantees. Section 5 provides numerical studies on synthetic data. Section 6 applies the method to real-world datasets. Section 7 concludes and discusses future work. Notations. For a vector v, let kvkq denote its ℓq -norm; we drop the subscript q when q = 2. For any matrix M, M⊤ is its transpose, kMk its spectral norm, kMkF its Frobenius norm, kMk2→∞ the maximum ℓ2 -norm of its rows, max(0, M) its entrywise maximum with zero, σi (M) its i-th largest singular value, λi (M) its i-th largest eigenvalue (ordered by magnitude), κ(M) = σ1 (M)/σK (M) its condition number, and rank(M) its rank. For any positive integer m, Im is the m × m identity matrix and [m] = {1, 2, . . . , m}. ei is the i-th standard basis vector. E[·] denotes the expectation operator.
2. Model Specification Let X = [x1 , . . . , xn ] ∈ R p×n be the observed data matrix, where p ≥ 2 is the dimension of each observation and n ≥ 2 is the number of samples (each sample is also referred to as an individual). The number of latent components is denoted by K and is assumed to be known, where K satisfies 1 ≤ K ≤ min{p, n}. To allow each individual to belong to multiple components with fractional weights, we introduce a mixed membership matrix Π ∈ [0, 1]n×K whose rows πi satisfy: πik ≥ 0,
K X
πik = 1
k=1
∀i ∈ [n], ∀k ∈ [K].
Thus, the i-th individual can have a distribution over the K latent groups, a flexible extension of the classical hard assignment. We say that an individual is pure if its membership vector πi is a standard basis vector (i.e., exactly one entry equals 1, and the remaining K − 1 entries are 0); otherwise, the individual is mixed. Pure individuals correspond to those who belong exclusively to a single component, while mixed individuals exhibit fractional membership across multiple components. The component centres are collected in Θ = [θ1 , θ2 , . . . , θK ] ∈ R p×K , with θk ∈ R p being the centre of the k-th component. The noise matrix E = [ǫ 1 , ǫ 2 , . . . , ǫ n ] ∈ R p×n has independent entries ǫi j that are sub-Gaussian: there exists a constant η > 0 such that
kǫ i j kψ2 := inf t > 0 : E[exp(ǫ 2i j /t2 )] ≤ 2 ≤ η, and consequently Var(ǫ i j ) ≤ C0 η2 for an absolute constant C0 (for example C0 = 2). With these ingredients, we now formally define the proposed model. The idea is simple: each observation is the sum of a weighted combination of the component centres (the weights given by the individual’s membership vector) plus a sub-Gaussian noise term. This leads to the following definition. 4
Definition 1 (Mixed Membership sub-Gaussian Model (MMSG)). Under the Mixed Membership sub-Gaussian Model, the observed data matrix X is generated in the following way:
X = ΘΠ⊤ + E,
E[X] =: P = ΘΠ⊤ .
(1)
The matrix P represents the signal part, and the goal is to recover the mixed membership matrix Π (up to a permutation of the components) from the observed data X alone. To the best of our knowledge, this model has not been studied before. It is the first to combine sub-Gaussian noise with a full mixed membership structure. It naturally captures overlapping clusters, a feature absent from classical mixture models. Two important special cases are worth noting. If every individual is pure (i.e., each row of Π is a standard basis vector), then the model reduces to the classical sub-Gaussian mixture model. If, in addition, the noise entries are independent and identically distributed Gaussian with variance σ2 , we recover the classical Gaussian mixture model (GMM), which has been extensively analyzed in the literature (Löffler et al., 2021; Ndaoud, 2022; Chen & Yang, 2021). Hence, the proposed MMSG model provides a unified framework that bridges hard clustering and overlapping clustering under sub-Gaussian noise.
Figure 1: Comparison between the classical GMM and the proposed MMSG on a synthetic toy example with K = 3 components, isotropic sub-Gaussian noise of magnitude 0.8, and n = 300 observations. Left panel: GMM with hard assignment. Each observation belongs to exactly one component; colours indicate the true component (Comp1, Comp2, Comp3). Dashed circles enclose the points assigned to each component, showing the disjoint compact clusters enforced by the single-membership constraint. Right panel: MMSG with mixed membership. Pure individuals (80% of the data, circles) are coloured by their unique component. Mixed individuals (20%, red squares) have fractional membership vectors πi and therefore lie in the overlap regions between the dashed circles. The component centres (black crosses) are identical in both panels. Unlike GMM, MMSG allows observations to occupy intermediate positions, reflecting a smooth transition between clusters that cannot be represented by any hard-assignment model.
Figure 1 illustrates the fundamental difference between the two models using the same component centres and 5
noise level. In the left panel, the classical GMM forces a strict partition: every point is assigned to exactly one of the three components, and the dashed circles (drawn to contain all points of each component) are well separated. In the right panel, the MMSG model relaxes this restriction. Here, 80% of the individuals are pure (circles) and serve as anchors for the three components. The remaining 20% are mixed (red squares). Their locations are generated as a weighted combination of the component centres according to their membership vectors πi plus sub-Gaussian noise. As a result, these mixed points naturally populate the regions between the dashed circles, creating a gradual transition from one component to another. This toy example demonstrates concretely that the MMSG model can represent overlapping cluster structures and fractional memberships, whereas the GMM cannot. Before proceeding, we discuss identifiability of the parameters (Θ, Π). A standard condition in mixed membership models is the existence of pure individuals (Mao et al., 2021; Jin et al., 2024; Chen & Gu, 2024). Specifically, we assume there exists an index set I ⊂ [n] with |I| = K such that, after a suitable permutation of the component labels, Π(I, :) = IK . That is, each component has at least one pure individual belonging exclusively to that component. Under this pure-individual condition, the following identifiability result holds. The proof follows the same reasoning as Theorem 2 of (Chen & Gu, 2024) for the grade-of-membership model, relying only on the low-rank factorization P = ΘΠ⊤ and the existence of pure individuals. Proposition 1. Consider the MMSG model defined in Equation (1) with the mixed membership matrix Π having nonnegative rows summing to one and each component having at least one pure individual. Then: (a) If rank(Θ) = K, the model is identifiable up to a permutation of the K components. (b) If rank(Θ) = K − 1 and no column of Θ is an affine combination of the other columns, the model is also identifiable up to a permutation. (c) In any other case, if there exists an individual i with πik > 0 for every k, then the model is not identifiable. By Proposition 1, our MMSG model is well-defined and identifiable under the conditions stated. In this paper, to facilitate theoretical analysis, we assume that rank(Θ) = K. Consequently, the model is identifiable. The difficulty of recovering the mixed membership matrix Π depends critically on two quantities: the minimum distance ∆ between component centres and the balancedness β of component sizes. In the classical GMM literature, these parameters determine the fundamental limits of exact recovery (Löffler et al., 2021; Ndaoud, 2022; Chen & Yang, 2021). A larger ∆ makes the components easier to separate, while a larger β prevents any single component from being too small and therefore hard to estimate. Specifically, the two key parameters of our MMSG model are defined as
∆ = min kθk − θℓ k, k,ℓ
6
β=
σ2K (Π) . n/K
We should emphasize that our theoretical analysis permits β to approach zero, thereby accommodating arbitrarily small clusters (including the extreme case where some components have vanishing proportions). This flexibility is essential for handling unbalanced data, a common scenario in real applications. Although our definition of β uses σ2K (Π) and thus appears different from that in (Löffler et al., 2021) and related works, this is because we are operating in the mixed membership setting. When the MMSG model reduces to the classical GMM (that is, when all individuals are pure), our β degenerates exactly to the balancedness parameter used in (Löffler et al., 2021).
3. Estimation Procedure Given the MMSG model X = ΘΠ⊤ + E, recovering the mixed membership matrix Π from the observed data X can be tackled from a spectral perspective. The signal matrix P = ΘΠ⊤ has rank K. Its right singular vectors, as we will show below, have a simplex geometry that directly reveals the pure individuals and, consequently, the membership vectors. The challenge lies in extracting this geometry from the noisy observations X without knowing P. To build intuition, suppose first that P itself were available. Write its compact singular value decomposition as P = VΣU⊤ , where V ∈ R p×K and U ∈ Rn×K have orthonormal columns, and Σ = diag(σ1 , σ2 , . . . , σK ) with σ1 ≥ σ2 ≥ · · · ≥ σK > 0. The following lemma, whose proof relies only on the existence of pure individuals (one per component), reveals a crucial geometric property. Lemma 1. Under the MMSG model with the pure individual condition, there exists an invertible matrix B ∈ RK×K such that
U = ΠB,
B = U(I, :).
Consequently, each row of U is a convex combination of the K rows B1,: , B2,: , . . . , BK,: . Thus, the rows of U live in a simplex in RK , and the vertices of that simplex are precisely the rows of U that correspond to pure individuals. This simplex structure is a recurring theme in mixed membership models: it appears in mixed membership stochastic blockmodels for network data (Mao et al., 2021; Jin et al., 2024), in grade of membership models for ordinal categorical data (Qing, 2024a,b), and in topic models (Ke & Wang, 2024). If we can locate those vertices, we can recover B and then solve for Π from U = ΠB. Because each row of Π sums to one and is nonnegative, the final step is a simple row-wise normalization. Thus, in the oracle setting where P is known, we would therefore: 1. Compute the top K right singular vectors U of P. 7
2. Apply a vertex hunting algorithm, such as the successive projection algorithm (SPA) (Gillis & Vavasis, 2013), to the rows of U to obtain the vertex indices I. 3. Set B = U(I, :) and compute Z = UB−1 . 4. For each i = 1, 2, . . . , n, set Πi,: = Zi,: /kZi,: k1 . This ideal procedure recovers Π exactly up to a label permutation since the SPA algorithm can exactly recover the true vertex indices I as long as the mixed membership matrix Π satisfies the definitions used in MMSG and the pure-individual condition (Mao et al., 2021). Of course, P is not observed in practice, and we only observe the data matrix X. A naive replacement of P by X fails because the diagonal entries of X⊤ X contain the squared noise magnitudes, which dominate the signal when the dimension p is large. Under heteroscedastic noise (which is allowed in our model), the leading eigenvectors can be strongly affected. To remove this bias, we consider the Gram matrix with its diagonal entries set to zero. Define the off-diagonal operator Poff-diag that sets all diagonal entries of a square matrix to zero while leaving off-diagonal entries unchanged. Then form G = Poff-diag (X⊤ X). For any two distinct individuals i , j, the entry Gi j = x⊤i x j expands as Gi j = p⊤i p j + p⊤i ǫ j + ǫ ⊤i p j + ǫ ⊤i ǫ j . Under the sub-Gaussian assumption, the cross terms p⊤i ǫ j and ǫ ⊤i p j have mean zero and are well concentrated, while ǫ ⊤i ǫ j also has mean zero for i , j. Therefore, for i , j, E[Gi j ] = p⊤i p j = (P⊤ P)i j . The diagonal entries of G are zero by construction, whereas the diagonal entries of P⊤ P are kpi k2 . Hence E[G] equals P⊤ P minus its diagonal part. Removing the diagonal introduces a perturbation that shifts the eigenvalues, but under mild signal strength conditions, this perturbation does not change the subspace spanned by the leading eigenvectors in a harmful way. The key point is that the off-diagonal Gram matrix gives an unbiased estimate of the off-diagonal entries of P⊤ P while leaving out the noise-inflated diagonal. This hollowing technique has been successfully employed in various high-dimensional settings, including subspace estimation (Cai et al., 2021), Gaussian mixture models (Ndaoud, 2022), and a unified ℓ p theory of principal component analysis and spectral clustering (Abbe et al., 2022). We therefore adopt G as a replacement for P⊤ P in the subsequent spectral analysis. Let Û ∈ Rn×K be the matrix of orthonormal eigenvectors corresponding to the K largest eigenvalues of G, so that G = ÛΛ̂Û⊤ with Û⊤ Û = IK and Λ̂ = diag(λ̂1 , λ̂2 , . . . , λ̂K ). We then apply SPA to the rows of Û to obtain an estimated vertex index set Î. By Lemma 1, the rows of U form a perfect simplex. The rows of Û form a noisy version of that 8
simplex. As long as the perturbation is not too large, SPA will still correctly identify the vertices (or at least produce indices close to them). With Î in hand, we form an estimate of the transformation matrix: Ẑ = max 0, Û Û(Î, :)−1 , where the entrywise maximum with zero enforces the nonnegativity of the membership weights. Finally, we normalize each row to unit ℓ1 norm to obtain the estimated mixed membership matrix:
Π̂i,: =
Ẑi,: kẐi,: k1
i = 1, 2, . . . , n.
,
The complete algorithm, which we call SPG (Sequential Projection for sub-Gaussian), is presented as Algorithm 1, where the underlying vertex hunting routine SPA is given in Algorithm 2. Algorithm 1 Sequential Projection for sub-Gaussian (SPG) 1: Input: Data matrix X ∈ R p×n , number of components K. 2: Output: Estimated mixed membership matrix Π̂ ∈ Rn×K . 3: Form G = Poff-diag (X⊤ X). 4: Compute the top K eigen-decomposition G = ÛΛ̂Û⊤ with Û ∈ Rn×K , Û⊤ Û = IK . 5: Apply SPA (Algorithm 2) to the rows of Û to obtain Î ⊂ [n], |Î| = K. 6: Compute Ẑ = max(0, Û Û(Î, :)−1 ). 7: for i = 1, 2, . . . , n do 8: Set Π̂i,: = Ẑi,: /kẐi,: k1 . 9: end for 10: return Π̂.
Algorithm 2 Successive projection algorithm (SPA) (Gillis & Vavasis, 2013) 1: Input: Y ∈ Rn×K , integer K. 2: Output: Vertex index set K ⊂ [n], |K| = K. 3: K = ∅, R = Y, t = 1. 4: while t ≤ K do 5: i∗ = argmaxi kRi,: k2 (break ties arbitrarily). 6: K = K ∪ {i∗ }. 7: u = Ri∗ ,: . uu⊤ 8: R = In − kuk 2 R. 2 9: t = t + 1. 10: end while 11: return K. The computational cost of SPG is dominated by two steps: forming the Gram matrix X⊤ X and computing its leading K eigenvectors after diagonal removal. Computing X⊤ X directly requires O(n2 p) operations and O(n2 ) storage. The eigen-decomposition of the n × n matrix G can be performed in O(n3 ) time via a full decomposition, or in O(n2 K) 9
total time when only the top K eigenvectors are needed. The storage cost of G is O(n2 ). The successive projection algorithm (SPA) runs in O(nK 2 ) time: each of its K iterations finds the row with maximum ℓ2 norm (O(nK)) and updates the residual matrix via a rank-1 projection (O(nK)). This cost is negligible when K ≪ n. In summary, the overall time complexity of SPG is O(n2 p + n2 K + nK 2 ), and the space complexity is O(n2 ). When n is very large, approximation techniques can reduce both measures, but a detailed discussion is beyond the present scope.
4. Theoretical Guarantees In this section, we establish the vanishing estimation error property of the spectral estimator Π̂ returned by Algorithm 1 under the proposed mixed membership sub-Gaussian model. The analysis shows how the key model parameters—the minimum centre distance ∆, the balancedness β, the sub-Gaussian noise level η (all defined in Section 2), the condition numbers of the centre and membership matrices, and the incoherence of the signal matrix—together determine the estimation difficulty. Under mild sufficient conditions expressed in terms of these quantities, we prove that the row-wise ℓ1 error vanishes asymptotically with high probability. We first define the condition numbers that appear in the theorem and its proof. Set
κ=
σ1 (Θ) , σK (Θ)
κΠ =
σ1 (Π) , σK (Π)
κP =
σ1 (P) . σK (P)
Recall that P = ΘΠ⊤ is the signal matrix and its singular value decomposition is P = VΣU⊤ with U ∈ Rn×K , V ∈ R p×K having orthonormal columns. The incoherence parameters related to the signal matrix P are defined as µ0 =
pn max j∈[p], i∈[n] |P j,i |2 kPk2F
,
µ1 =
n max kUi,: k2 , K i∈[n]
µ2 =
p max kV j,: k2 . K j∈[p]
Set d = max{n, p} and µ = max{µ0 , µ1 , µ2 }. These constants appear in the main theoretical results given below, and they are bounded by absolute constants in many practical scenarios, which leads to simplified scaling laws. We now present the main theorem. Theorem 1. Let the estimator Π̂ be produced by Algorithm 1. Suppose that the data dimensions satisfy 1
8 2 np ≫ µ2 κ8 κΠ K log4 d,
8
8 4 2 κ8 κΠ µ 3 κ 3 κΠ Kσ13 (Π), p≫ K log2 d, n ≫ 2 β β3
10
(2)
and the separation conditions s r p 7 5 21 2 14 ηκ 2 κΠ µ Kσ1 (Π) log d ηκ2 κΠ µ K p 1p ∆≫ , ∆≫ σ1 (Π) ( ) 4 K log d √ 1 1 n n n β2 β2
(3)
hold, then with probability at least 1 − O(d−10 ), we have max kΠ̂i,: − (ΠP)i,: k1 = o(1), i∈[n]
where P is a permutation matrix . Theorem 1 establishes that, under the stated dimension and separation conditions, the spectral estimator Π̂ returned by the proposed SPG approach achieves vanishing ℓ1 estimation error for the MMSG model. The explicit upper bound for maxi∈[n] kΠ̂i,: − (ΠP)i,:k1 is derived in the proof (see Equation (A.5) in the appendix), where it depends on all model parameters n, p, K, ∆, η, µ, κ, κΠ, β, σ1 (Π), and it is o(1) when conditions in this theorem hold. The sufficient conditions are grouped into two natural families: scaling of dimensions (see conditions in Equation (2)) and separation of component centres (see conditions in Equation (3)). The dimension scaling conditions ensure that the sample size n, the ambient dimension p, and their product are sufficiently large relative to key problem parameters. These conditions appear complex because we explicitly track the influence of all model parameters (incoherence, condition numbers, balancedness, dimension ratios, and logarithmic factors). The separation conditions are required to ensure vanishing ℓ1 estimation error of mixed memberships for each sample. They impose a lower bound on the minimal centre distance ∆ that grows when the estimation problem becomes harder. Based on the explicit form of the lower bounds in conditions in Equation (3), we can analyze how each parameter affects the required ∆. • When the balancedness β = σ2K (Π)/(n/K) becomes smaller, the membership matrix Π is more ill-conditioned because its smallest singular value σK (Π) is small. This makes the problem harder, and indeed both lower √ bounds on ∆ contain a factor 1/ β. Hence, a smaller β forces a larger ∆ for vanishing estimation error of mixed memberships. • When the number of components K increases, the problem becomes more difficult because there are more centres to separate. Therefore, a larger K requires a larger ∆ for vanishing estimation error of mixed memberships. • When the condition number κ = σ1 (Θ)/σK (Θ) of the centre matrix increases, the centres are harder to distinguish. The first lower bound has κ7/2 and the second has κ2 . Thus, a larger κ demands a larger ∆. • When the condition number κΠ = σ1 (Π)/σK (Π) of the membership matrix increases, the rows of Π are more 11
5 skewed, making the pure individual recovery more challenging. The first lower bound has κΠ and the second 2 has κΠ . Hence, a larger κΠ also requires a larger ∆.
• When the noise level η increases, the data are noisier. Both lower bounds are linear in η, so a larger η directly forces a larger ∆. In summary, all these parameters affect the required separation ∆ monotonically: larger difficulty (smaller β, larger K, larger κ, larger κΠ , larger η) leads to a larger necessary ∆. This is precisely reflected in the explicit inequalities of conditions in Equation (3). Under these sufficient conditions, with probability at least 1 − O(d−10 ) (which is very high for large n and p), the estimated membership matrix satisfies maxi∈[n] kΠ̂i,: − (ΠP)i,:k1 = o(1), which means that up to a global relabelling of the components, the row-wise ℓ1 estimation error vanishes for every individual under mild conditions on ∆. The sufficient conditions in Theorem 1 involve several model-specific parameters. In many practical scenarios, these parameters are naturally bounded by absolute constants, leading to a dramatically simplified yet still rigorous set of conditions. The following corollary makes this precise and reveals the essential scaling laws that guarantee the vanishing estimation error. Corollary 1. Assume that the model parameters satisfy the following boundedness conditions:
κ = O(1),
β = O(1),
κΠ = O(1),
η = O(1),
µ = O(1).
Then, the sufficient conditions in Theorem 1 reduce to
np ≫ K 2 log4 d,
p ≫ K log2 d, ( 1/4 ) p p . ∆ ≫ K log d max 1, n
n ≫ K,
(4) (5)
Under these conditions, the same conclusion holds: with probability at least 1 − O(d−10 ), K max kΠ̂i,: − (ΠP)i,: k1 . + i∈[n] n
p
K log d K log d + ∆ ∆2
r
p = o(1), n
where P is a permutation matrix. The boundedness assumptions in Corollary 1 are remarkably mild and hold in a wide range of applications. • Condition number κ = O(1): This means that the ratio is bounded above by an absolute constant that does not depend on the sample size n, the ambient dimension p, or the number of components K. A bounded κ 12
ensures that the columns of the centre matrix Θ are well-conditioned. Equivalently, the component centres are not nearly linearly dependent. The notion of linear dependence here is global: even if some centres are statistically correlated in the sense of having small pairwise angles, the entire set may still be well-conditioned as long as no centre can be approximated as a linear combination of the others. The requirement κ = O(1) excludes degenerate configurations where the centres almost lie in a proper subspace of R p , a situation that would compromise identifiability and fundamentally hinder the possibility of achieving vanishing estimation error. This is a standard assumption in the analysis of mixture models, and it is considerably weaker than requiring the centres to be orthogonal or far apart. • Balancedness β = O(1): Recall that β = σ2K (Π)/(n/K). Here β = O(1) means that β is bounded above and below by positive absolute constants. Equivalently, there exist constants 0 < c1 ≤ c2 < ∞ such that c1 ≤ β ≤ c2 . Hence √ √ √ σK (Π) satisfies c1 n/K ≤ σK (Π) ≤ c2 n/K, i.e., the smallest singular value of Π is exactly of order n/K. This ensures that the membership vectors are well spread across the probability simplex, and no component receives asymptotically vanishing total weight. Such a balancedness condition is natural for mixed membership models: each latent component must appear with a substantial overall contribution. When combined with the √ additional assumption κΠ = σ1 (Π)/σK (Π) = O(1), we also obtain σ1 (Π) = O( n/K). Therefore, all singular √ values of Π are of the same order n/K, meaning that Π is well-conditioned. • Noise level η = O(1): The sub-Gaussian norm of the noise entries is bounded by a constant. This is the typical setting in high-dimensional statistics: any fixed noise variance qualifies. • Incoherence µ = O(1). The parameter µ controls the maximum leverage of the row and column spaces of the signal matrix P = ΘΠ⊤ . A bounded µ follows from explicit and verifiable assumptions on both Π and Θ. For √ the membership matrix, we require that the smallest singular value satisfies σK (Π) ≍ n/K, which is equivalent to the balancedness parameter β = σ2K (Π)/(n/K) being bounded below and above by positive constants. Under this condition, Lemma 4 directly gives µ1 = O(1). For the centre matrix Θ, we impose two mild regularity conditions. First, the entries of Θ are uniformly bounded by an absolute constant, and each centre vector θk has √ Euclidean norm at least of order p (so that the signal does not vanish as the dimension grows). Second, the column space of Θ is incoherent: the projection of any standard basis vector onto this space has Euclidean norm p at most C K/p for some absolute constant C. Under these conditions, a standard calculation shows that both µ0 and µ2 are bounded by absolute constants. Consequently, µ = max{µ0 , µ1 , µ2 } = O(1). In many practical
scenarios, such as when the rows of Π are drawn independently from a Dirichlet distribution with balanced concentration parameters (ensuring β = O(1) with high probability) and the centres are taken as, for example, 13
random matrices with independent bounded entries (so that the norm condition and the incoherence condition hold with high probability), the above assumptions are satisfied. Hence, the boundedness of µ is a rigorous consequence of explicit and interpretable conditions, not an ad-hoc postulate. Under these boundedness assumptions, the original complicated sufficient conditions collapse to the transparent scaling laws in Equations (4) and (5). Let us interpret them: • The sample size n and the dimension p need only grow polynomially in the number of components K and logarithmically in d = max{n, p}. Specifically, n ≫ K ensures that we have enough observations to estimate the K pure individuals. The condition p ≫ K log2 d guarantees, with high probability, that the operator norm of the sub-Gaussian noise matrix is sufficiently small relative to the smallest singular value of the signal matrix P. This ensures that the spectral perturbation caused by the noise does not destroy the low-rank structure of the signal, thereby enabling reliable recovery of the mixed membership vectors with vanishing error. The product condition np ≫ K 2 log4 d guarantees that the spectral norm of the noise contribution to the off-diagonal Gram matrix is sufficiently small relative to the signal, thereby controlling the eigenvector perturbation and ensuring accurate subspace estimation. • The minimal Euclidean distance ∆ between any two distinct component centres must dominate
p
K log d times
a factor that depends on the aspect ratio p/n. When p and n are of the same order, this factor is constant, p so ∆ ≫ K log d. If p is much larger than n, the condition becomes slightly more stringent because the
extra dimensions amplify the noise. The factor (p/n)1/4 captures this effect. This requirement is remarkably mild. For the canonical two-component case with balanced sample sizes or fixed K case, our condition reduces p to ∆ ≫ log n (up to constants after fixing the noise level). This matches the sharp information-theoretic threshold for exact recovery established in the classical Gaussian mixture model literature (Löffler et al., 2021; Chen & Yang, 2021; Ndaoud, 2022).
5. Numerical Experiments We now describe the simulation design used to validate the theoretical vanishing-error property of the SPG estimator (Algorithm 1) under the mixed membership sub-Gaussian model (MMSG). The primary goal is to verify that the row-wise ℓ1 estimation error vanishes under the conditions of Theorem 1, and to compare SPG with the weighted spectral clustering (WSC) for classical Gaussian mixtures studied in (Löffler et al., 2021; Zhang & Zhou, 2024) (which assumes a single membership per observation). For each parameter configuration, we generate 200
14
independent data sets and report averages of the aligned row-wise ℓ1 error rate (Hamming error rate):
Error rate =
n X 1 min Π̂i,: − (ΠP)i,: 1 , n P∈SK i=1
where SK is the set of all K × K permutation matrices 5.1. Data generation For each simulation run, we first construct the component centre matrix Θ = [θ1 , θ2 , . . . , θK ] ∈ R p×K . To achieve equal pairwise distances ∆ between centres, we take θk = √∆2 uk , where u1 , u2 , . . . , uK are orthonormal vectors obtained from the QR decomposition of a random p × K matrix with independent standard normal entries. This guarantees kθk − θℓ k = ∆ for all k , ℓ. The separation ∆ is chosen to satisfy the theoretical lower bounds. In the simplified setting of Corollary 1 (where boundedness assumptions hold), we set
∆ = c∆ ·
p
n o K log d max 1, p/n 1/4 ,
d = max{n, p},
with a constant c∆ that we vary to explore the phase transition. For experiments that verify the vanishing estimation error property, we take c∆ = 10, which comfortably exceeds the required lower bound. Importantly, ∆ is recomputed for every combination of n, p, K according to this formula, so that the theoretical separation condition is always satisfied. The mixed membership matrix Π ∈ [0, 1]n×K has rows summing to one. We fix a total number of pure individuals npure and distribute them equally among the K components (so each component receives npure /K pure rows, rounding as needed). The remaining n − npure rows are mixed. For mixed rows we generate independent membership vectors πi from a Dirichlet distribution Dirichlet(α1K ) with concentration parameter α > 0. This construction allows us to control the balancedness parameter β = σ2K (Π)/(n/K) via α: smaller α yields more extreme memberships (closer to pure), which increases σK (Π) and hence β; larger α produces more uniform mixtures and reduces β. Noise is generated according to the following two scenarios: • Gaussian heteroscedastic (GHe): Each entry ǫi j is independent N(0, σ2i ), where σi is drawn uniformly from [0.5, η]. • Sub-Gaussian heteroscedastic (SHe): We set ǫi j = σi ri j with i.i.d. Rademacher variables ri j ∈ {+1, −1} equiprobably and σi ∼ Uniform(0.5, η). All experiments below are conducted for both noise scenarios (GHe and SHe) and both dimension regimes: the 15
low-dimensional regime where p ≪ n and the high-dimensional regime where p ≫ n. Specific parameter choices for each regime are given within each experiment. 5.1.1. Experiment 1: Varying sample size n
0.5
0.4 SPG (GHe) SPG (SHe) WSC (GHe) WSC (SHe)
0.3
0.2
0.1
0 500
1000 1500 2000 2500 3000 3500 4000 4500 5000
Figure 2: Numerical results of Experiment 1
We examine the effect of the sample size n while keeping the centre separation ∆ constant. Fix the feature dimension at p = 2000 and choose a sufficiently large separation to satisfy the theoretical condition for all considered p n. Specifically, set ∆ = 10 · K log dmax · max{1, (p/nmin)1/4 }, where K = 4, nmin = 500 and dmax = max{5000, 2000}.
This yields a constant ∆ that exceeds the required lower bound for every n in the range. Let n ∈ {500, 1000, . . ., 5000}. This range covers both the high-dimensional regime (n < p) and the low-dimensional regime (n > p). To keep the pure proportion constant, we set the number of pure individuals to npure = ⌊0.4n⌋ (so that 40% of the individuals are pure, equally distributed among the K = 4 components). The other parameters are α = 0.5 and η = 1. This design isolates the effect of sample size on the estimation error, because the separation strength ∆ does not vary with n. Figure 2 displays the results. The proposed SPG estimator achieves an error below 0.1 for all n, with no visible decline as n grows. In contrast, WSC yields an error between 0.4 and 0.5, also flat across n. These patterns are direct consequences of the model structure and the theoretical bounds. Under the conditions of Corollary 1, the estimation √ q K log d p K log d error of SPG is bounded (up to universal constants) by Kn + + ∆ n , where d = max{n, p}. The first ∆2
and third terms decrease as n increases. Because ∆ is held fixed, the middle term remains essentially constant and dominates the bound once n is moderately large. Hence, SPG’s error stabilises at a small non-zero level, exactly as predicted. The error would vanish only if ∆ were sufficiently large to satisfy the separation condition. Here ∆ is
fixed, so the error approaches a constant rather than zero. For WSC, the hard assignment assumption is fundamentally mismatched with mixed membership data. The many mixed individuals (60% of the sample) cannot be assigned to 16
any single component, so any hard clustering has a bias that does not vanish. This bias also does not decrease as n grows, which explains why the error of WSC remains near 0.5 throughout. 5.1.2. Experiment 2: Varying the separation strength c∆ Experiment 2: n=200, p=2000
Experiment 2: n=2000, p=20 0.5
0.5
0.4
0.4 SPG (GHe) SPG (SHe) WSC (GHe) WSC (SHe)
0.3
SPG (GHe) SPG (SHe) WSC (GHe) WSC (SHe)
0.3
0.2
0.2
0.1
0.1
0
0 0
20
40
60
80
100
0
20
40
60
80
100
Figure 3: Numerical results of Experiment 2.
To investigate the sharpness of the separation condition, we consider two representative settings: a high-dimensional setting with n = 200, p = 2000 (so n ≪ p) and a low-dimensional setting with n = 2000, p = 20 (so n ≫ p). In both settings, we fix K = 4, α = 0.5, η = 1, and maintain a constant pure proportion of 40% by setting npure = ⌊0.4n⌋. We let c∆ vary from 10 to 100 in steps of 10. This allows us to observe the phase transition predicted by Theorem 1: below a certain threshold, the error is large, while above it the error becomes small. Figure 3 plots the estimation error against the separation factor c∆ . For the SPG estimator, the error falls below 0.1 even at the smallest c∆ and continues to decline toward zero as c∆ increases. For WSC, its error stays almost constant between 0.4 and 0.5, showing almost no response to larger centre distances. As Theorem 1 predicts, a larger c∆ increases the minimal centre distance ∆, which reduces noise perturbation. The spectral estimate then becomes more accurate, and the recovered membership vectors converge to the truth. Hence, SPG’s error vanishes for sufficiently large ∆. In contrast, WSC assumes hard assignments. Even with infinite centre separation, a mixed individual incurs an ℓ1 error that depends only on its true fractional membership vector and has a positive lower bound that does not vanish. Because the distribution of mixed membership vectors is fixed, the average error from mixed individuals is constant. Increasing ∆ improves pure individual classification but cannot reduce this irreducible error from the mixed majority, which constitutes 60% of the sample. Consequently, the overall WSC’s error remains unchanged as c∆ increases. This behaviour shows that mixed membership is an intrinsic property of the data, not a lack of separation. Estimating 17
fractional memberships, therefore, requires a dedicated method such as our SPG. 5.1.3. Experiment 3: Changing the balancedness β via the Dirichlet concentration α Experiment 3: n=200, p=2000 0.9
Experiment 3: n=2000, p=200 0.9
SPG (GHe) SPG (SHe) WSC (GHe) WSC (SHe)
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0.4
0.45
0.5
0.55
0.6
0.65
SPG (GHe) SPG (SHe) WSC (GHe) WSC (SHe)
0 0.4
0.7
0.45
0.5
0.55
0.6
0.65
0.7
0.75
Figure 4: Numerical results of Experiment 3.
We vary the Dirichlet concentration parameter α ∈ {0.2, 0.5, 1.0, 2.0, 5.0} to control the balancedness β = σ2K (Π)/(n/K). Smaller α produces more extreme membership vectors (closer to pure), which increases σK (Π) and hence β; larger α yields more uniform mixtures and reduces β. We consider two representative settings: a high-dimensional setting with n = 200, p = 2000 and a low-dimensional setting with n = 2000, p = 200. In both settings, we fix K = 4, c∆ = 10, η = 1, and maintain a constant pure proportion of 40% by setting npure = ⌊0.4n⌋. For each α we compute the empirical β from the generated Π and record the average row-wise ℓ1 error. This experiment examines how the estimation error depends on the balancedness of the membership matrix. Figure 4 reports the estimation errors of SPG and WSC as the balancedness parameter β changes. Recall that β = σ2K (Π)/(n/K). A larger β indicates more extreme membership vectors, in which most individuals are nearly pure. Smaller β corresponds to more uniform mixtures, with fractional memberships spread across components. The figure reveals two clear patterns. The estimation error of SPG stays below 0.1 across all β values and exhibits no systematic trend. In contrast, the error of WSC is much larger, ranging between 0.3 and 0.9, yet it decreases steadily as β grows. These results follow directly from the model structure and the algorithms. SPG is designed to recover fractional membership vectors. As Corollary 1 shows, under the maintained boundedness assumptions (which hold for all configurations in this experiment), the theoretical error bound for SPG does not involve β. In this experiment, the separation constant c∆ = 10 is fixed and sufficiently large to satisfy the theory for every β. Hence, the noise perturbation is well controlled, the estimated eigenvectors stay close to the true simplex, and the recovery of Π is equally accurate for all β. This explains why SPG’s error is stable and uniformly low. 18
WSC assumes a single component per observation, but this assumption is incorrect for mixed-membership data. When the overall balancedness β is small, the membership vectors across the population are nearly uniform. In this regime, a typical mixed individual lies near the centre of the simplex, far from any pure centre. Forcing such an individual into a single cluster creates a large ℓ1 error. As β increases, the population shifts toward more extreme membership vectors. Many individuals become nearly pure, and even the mixed ones move closer to the vertices. Hence, WSC can correctly classify a larger fraction of pure individuals, and the cost of forced assignments decreases. WSC’s error, therefore, declines with β. Yet even at the largest β, the error remains above 0.3, far exceeding SPG’s error. The remaining gap comes from the mixed individuals that are always present, as the proportion of pure individuals is fixed at 40%. A hard assignment can never match a true fractional vector, so the error cannot vanish. 5.1.4. Experiment 4: Varying the proportion of pure individuals Experiment 4: n=200, p=2000
Experiment 4: n=2000, p=200 SPG (GHe) SPG (SHe) WSC (GHe) WSC (SHe)
0.7 0.6
SPG (GHe) SPG (SHe) WSC (GHe) WSC (SHe)
0.7 0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0 0
0.1
0.2
0.3
0.4
0.5
0
0.1
0.2
0.3
0.4
0.5
Figure 5: Numerical results of Experiment 4.
We investigate how the estimation error depends on the proportion of pure individuals. Fix the total sample size n and feature dimension p for two regimes: a low-dimensional regime with n = 2000, p = 20 and a high-dimensional regime with n = 200, p = 2000. In each regime, we vary the pure proportion cpure ∈ {0.05, 0.1, . . . , 0.5} and set the number of pure individuals to npure = ⌊cpure n⌋, with the constraint that npure ≥ K (so each component has at least one pure individual). Pure individuals are equally distributed among the K = 4 components. The remaining n − npure individuals are mixed, generated from a Dirichlet distribution with α = 0.5. The other parameters are K = 4, c∆ = 10, η = 1. Figure 5 plots the estimation error against the pure proportion cpure . The SPG estimator maintains an error below 0.1 across the entire range, showing almost no dependence on cpure . The WSC estimator, by contrast, starts with an error near 0.8 when only 5% of the individuals are pure and declines steadily to about 0.4 when half of the sample is 19
pure. Despite this decline, WSC’s error remains much larger than that of SPG for every value of cpure . SPG’s insensitivity to cpure follows from two facts. First, the estimation error of SPG is primarily controlled by the noise level and the centre separation ∆, both held constant in this experiment. Second, the vertex hunting step requires only one pure individual per component. Adding more pure individuals does not reduce the noise, increase ∆, or improve the simplex identification beyond what a single pure individual already provides. Hence, SPG’s error stays flat. For WSC, the decline is straightforward. Pure individuals are correctly classified when the centres are well separated, as they are here with c∆ = 10. Mixed individuals, generated with α = 0.5, contribute a roughly constant error that does not vanish. As cpure increases, the average error falls because pure individuals dominate the sample, yet the error never reaches the level of SPG because the mixed individuals remain. It should be emphasized that this experiment differs from Experiment 3. In Experiment 3, we varied α, which changes the distribution of membership vectors among the mixed individuals themselves. Here α is fixed at 0.5, so the mixed individuals have the same mixing behaviour throughout; only the number of pure individuals changes. The fact that SPG performs equally well in both settings, while WSC improves in both settings but always remains far behind, reinforces the conclusion that fractional memberships require a method specifically designed to estimate them.
6. Real Data Applications We now consider three real datasets that have long served as standard benchmarks for clustering: Iris, Wine, and Dermatology 1 . Each comes with a set of hard labels assigned by domain experts. This conventional perspective forces every observation into a single class. Such a hard assignment, however, may ignore the possibility that some individuals have mixed memberships. Our goal is to assess whether mixed membership patterns exist in these data. Under our model, each observation is equipped with a membership vector whose entries are non-negative and sum to one. This vector can be degenerate (a pure individual) or have multiple positive entries (a mixed individual). By re-examining these datasets, we can identify which individuals, if any, exhibit mixed membership characteristics. Table 1 summarises the basic information of the three datasets. The Iris data, from Fisher (1936), contains three iris species described by four morphological measurements. The Wine data (Aeberhard et al., 1994) consists of three wine cultivars with 13 chemical attributes such as alcohol, malic acid, and ash. The Dermatology data (Güvenir et al., 1998) covers six erythemato-squamous diseases with 34 clinical and histopathological features. Each dataset has a known ground truth number of latent classes. 1 The three real datasets are available for download at the UCI Machine Learning Repository: https://archive.ics.uci.edu.
20
Dataset Iris Wine Dermatology
Table 1: Basic information of the three real datasets. Sample meaning Feature meaning
Source Fisher (1936) Aeberhard et al. (1994) Güvenir et al. (1998)
Iris flower Wine sample Skin lesion
Sepal and petal dimensions Chemical concentrations (alcohol, malic acid, ash, etc.) Clinical and histopathological scores
n
p
K
150 178 358
4 13 34
3 3 6
We apply the SPG algorithm to each dataset. The output is an estimated membership matrix Π̂ whose rows sum to one. For each individual, we define its home base component as ĉi = arg maxk Π̂ik . We call an individual highly pure if maxk Π̂ik ≥ 0.9 and highly mixed if maxk Π̂ik ≤ 0.6. Let τpure and τmixed denote the proportions of such individuals. We also compute the condition number κ(Π̂) = σ1 (Π̂)/σK (Π̂). This quantity measures how balanced the estimated membership vectors are across the components. A value close to one indicates that the membership weights are distributed relatively evenly, while larger values suggest that some components dominate the memberships for most individuals. Table 2: Estimated mixing characteristics for the three real datasets.
Dataset
τpure
τmixed
κ(Π̂)
Iris Wine Dermatology
0.2467 0.5506 0.3324
0.0600 0.1685 0.2737
7.6167 1.2813 2.1764
The estimated mixing characteristics are reported in Table 2. Several observations follow: • For the Iris dataset, the proportion of highly pure individuals is 0.2467 and the proportion of highly mixed individuals is 0.0600. The condition number κ(Π̂) = 7.6167 is large. This indicates that the estimated membership vectors are strongly unbalanced. Only a small fraction of observations serve as nearly pure anchors, while most individuals have moderate weights across the three components. The low mixing proportion suggests that few Iris flowers lie on the boundaries between species. • For the Wine dataset, the pure proportion is 0.5506 and the mixed proportion is 0.1685. The condition number is 1.2813, which is close to one and suggests that the membership vectors are well balanced across the three cultivars. The high pure proportion indicates that most wines have a clear dominant cultivar, while the mixed proportion captures those samples that lie on the boundaries between cultivars. • For the Dermatology dataset, the pure proportion is 0.3324, the mixed proportion is 0.2737, and the condition number is 2.1764. This dataset exhibits the highest mixing proportion among the three. The pattern is consistent with the clinical reality that several erythemato-squamous diseases share overlapping features. About one-third of the patients are highly pure, and about one-quarter are highly mixed. The remaining individuals (nearly 40%) 21
fall into neither category, indicating that many patients have membership patterns that are moderately spread across several disease classes. To visualise the estimated memberships, we draw ternary plots for the two datasets with three components, namely Iris and Wine. Figure 6 presents the results. Each point represents an individual, and different symbols indicate three categories: highly pure (maximum membership ≥ 0.9), moderate (0.6 < maximum membership < 0.9), and highly mixed (maximum membership ≤ 0.6). The vertices of the triangle correspond to the purest samples. Highly mixed points cluster near the centre, moderate points appear closer to the edges or vertices, and highly pure points are located at the vertices. For the Iris data, the number of highly mixed points is very small, consistent with τmixed = 0.0600 in Table 2. Most moderate points concentrate along the edge between Components 1 and 3, with a visible bias toward Component 3. For the Wine data, highly mixed points are far more numerous (τmixed = 0.1685) and fill the central region of the triangle. Moderate points are widely spread across the simplex, showing no strong preference for any particular edge or vertex. These visual patterns confirm that the mixed membership structure differs substantially between the two datasets: Wine exhibits a larger and more evenly distributed fraction of mixed and moderate individuals.
Figure 6: Ternary plots of estimated membership vectors for the Iris data (left) and the Wine data (right). Colours indicate the home base component.
7. Conclusion In this paper, we have proposed a novel and interpretable statistical model, the mixed membership sub-Gaussian model (MMSG). Unlike the classical Gaussian mixture model, which forces each observation to belong to exactly one component, the MMSG model allows an observation to belong to multiple components simultaneously, with a membership vector that sums to one. This flexibility makes the model suitable for a wide range of applications where mixed membership is natural, such as in genetics, social sciences, and text analysis. For this model, we developed an efficient spectral algorithm, called SPG, to estimate the mixed memberships of individuals. The algorithm constructs 22
the off-diagonal Gram matrix, extracts its top eigenvectors, applies the successive projection algorithm to identify pure individuals, and then recovers the mixed membership matrix via row normalization. Under mild separation conditions on the component centres, we proved that the estimation error of the per-individual mixed membership vector can be made arbitrarily small with high probability. Extensive numerical experiments demonstrate that our method significantly outperforms existing algorithms that do not account for mixed membership. To the best of our knowledge, this is the first work to provide a computationally efficient estimator with such a vanishing-error guarantee for a mixed-membership extension of the Gaussian mixture model. The present work opens up several interesting and challenging directions for future research. First, the number of components K is assumed to be known throughout this paper. Developing methods with theoretical guarantees to estimate K from data under the MMSG model is a fundamental and challenging problem that remains largely open. Second, model selection between the classical GMM and the MMSG model is an important practical question. One may design information criteria, likelihood ratio tests, or cross-validation procedures to decide whether a hard clustering or a mixed membership structure better explains a given data set. Third, it would be interesting to develop a test for whether two individuals have exactly the same mixed membership vector, which is more challenging than testing for pure membership alone. Fourth, tensor methods that exploit higher-order moments could potentially improve statistical efficiency. Fifth, robust extensions that tolerate heavy-tailed noise or adversarial outliers are worth investigating, as real data may deviate from sub-Gaussian assumptions. Sixth, scaling the SPG algorithm to very large sample sizes using randomised sketching techniques is a promising direction for big data applications. Just as the classical mixed membership stochastic blockmodel has inspired a large body of research in network analysis, we expect that the MMSG model will serve as a foundation for many future studies in mixed membership modelling for continuous data, both in theory and in practice, and this paper is not an end but rather the beginning of a rich line of research.
CRediT authorship contribution statement Huan Qing is the sole author of this article.
Declaration of competing interest The author declares no competing interests.
Data availability Data and code will be made available on request. 23
Appendix A. Technical Proofs Appendix A.1. Proof of Proposition 1 Proof. The MMSG model shares exactly the same low-rank structure P = ΘΠ⊤ with the mixed membership matrix Π satisfying the pure-individual condition as the grade-of-membership (GoM) model studied in (Chen & Gu, 2024). Theorem 2 of (Chen & Gu, 2024) establishes the above identifiability statements for the GoM model. The proof of that theorem uses only the factorization P = ΘΠ⊤ , the row-sum-to-one and non-negativity constraints on Π, and the existence of pure individuals. It never requires the entries of Θ to lie in [0, 1]. Therefore, the same conclusions hold for the MMSG model without any modification. We omit the repetitive details and refer the reader to the original proof in (Chen & Gu, 2024) for a similar argument. Appendix A.2. Proof of Lemma 1 Proof. Without loss of generality, suppose that there exists a set I ⊂ [n] of pure individuals such that Π(I, :) = IK after a suitable permutation of the components. Take the compact singular value decomposition P = VΣU⊤ , where U⊤ U = IK , V⊤ V = IK , and Σ = diag(σ1 , σ2 , . . . , σK ) with σ1 ≥ σ2 ≥ · · · ≥ σK > 0. Then U = P⊤ VΣ−1 . Since P⊤ = ΠΘ⊤ , we obtain U = ΠΘ⊤ VΣ−1 .
Define B = Θ⊤ VΣ−1 . Because rank(Θ) = K and V shares the same column space as Θ (since col(P) = col(Θ)), there exists an invertible matrix R such that Θ = VR. Consequently, Θ⊤ V = R⊤ is invertible. Hence B is invertible, and we have U = ΠB. For the pure individuals, using Π(I, :) = IK gives U(I, :) = Π(I, :)B = B. Thus B = U(I, :). Finally, for any i ∈ [n], we have Ui,: = π⊤i B =
K X
πik Bk,: ,
k=1
where πik ≥ 0 and
PK
k=1 πik = 1, so each row of U is a convex combination of the vertex rows B1,: , B2,: , . . . , BK,: .
Appendix A.3. Proof of Theorem 1 Proof. When the following conditions hold:
np ≫ µ2 κP8 K 2 log4 d,
(A.1)
p ≫ µ1 κP8 K log2 d,
(A.2) 24
1 1 η ≪ min , , p p √ κP 4 np log d κ3 n log d σK (P)
(A.3)
P
n ≫ µ1 κP4 K,
(A.4)
Theorem 1 in (Cai et al., 2021) with sampling rate psamp = 1 guarantees that there exists an orthogonal matrix O ∈ −10 RK×K such that with probability at least 1 − O(dmax ),
kÛO − Uk2,∞ κP2
r
µK Egeneral , n
where
Egeneral =
µ1 κP2 K η2 √ ηκP p + n log d + 2 np log d. n σK (P) σK (P)
By Lemmas 2, 3, and 4, we know that σK (P) ≥ √∆2κ σK (Π), κP ≤ κκΠ , and µ1 ≤ β1 . Therefore, as long as the dimensions n, p, and the separation parameter ∆ satisfy
8 2 np ≫ µ2 κ8 κΠ K log4 d, p ≫
8 4 ηκ4 κ3 p κ8 κΠ κ4 κΠ ηκ2 κΠ p 1 p K log2 d, n ≫ K, ∆ ≫ √ ( ) 4 K log d, ∆ ≫ √ Π K log d, β β β n β
the conditions in Equations (A.1)-(A.4) hold naturally. Since Û′ Û = IK , U′ U = IK , by basic algebra, we have kÛÛ′ − UU′ k2→∞ ≤ 2kÛO − Uk2→∞ . Set ̟ = kÛÛ′ − UU′ k2→∞ . Since both SPG and the SPACL algorithm without the prune step of (Mao et al., 2021) apply the SPA algorithm to an eigenvector matrix to hunt for pure individuals and then estimate the mixed memberships, the proof of SPG’s estimation error is the same as SPACL. By Lemma Equation (3) in Theorem 3.2 of (Mao et al., 2021), there exists a K × K permutation matrix P such that r ! ! p µK 2 2 ′ Egeneral . max kΠ̂i,: − (ΠP)i,: k1 = O ̟κ(Π Π) λ1 (Π Π) = O κP κΠ σ1 (Π) i∈[n] n ′
Since σK (P) ≥ √∆2κ σK (Π), κP ≤ κκΠ , and µ1 ≤ β1 , we get Egeneral = ≤
µ1 κP2 K η2 √ ηκP p n log d + 2 np log d + n σK (P) σK (P)
2 K κ2 κΠ ηκκΠ p 2κ2 η2 √ + 2κn log d + 2 2 np log d βn ∆σK (Π) ∆ σK (Π)
25
=
2 κ2 κΠ K ηκκΠ p 2κ2 η2 K 2κK log d + + √ βn ∆2 β ∆ β
r
p log d, n
which gives
4 σ1 (Π) max kΠ̂i,: − (ΠP)i,: k1 = O κ2 κΠ i∈[n]
1
When n ≫
8
4 µ 3 κ 3 κΠ 2 β3
7
2 3
Kσ1 (Π), ∆ ≫
r
2 K ηκκΠ p µK κ2 κΠ 2κ2 η2 K 2κK log d + + √ n βn ∆2 β ∆ β
1
√
1
2 4 5 2 Kσ (Π) log d ηκ2 κΠ µ ηκ 2 κΠ µ 1 √ , and ∆ ≫ 1 1 n β2 β2
maxi∈[n] kΠ̂i,: − (ΠP)i,: k1 = o(1) in Equation (A.5).
r
σ1 (Π)
q
r
!! p log d . n
(A.5)
p K p 41 K log d, we have n (n)
Appendix A.4. Proof of Corollary 1 Proof. We start from the sufficient conditions in Theorem 1 and substitute the boundedness assumptions κ = O(1), β = O(1), κΠ = O(1), η = O(1), µ = O(1). Recall that β = σ2K (Π)/(n/K), so β = O(1) implies σ2K (Π) = O(n/K). Because κΠ = σ1 (Π)/σK (Π) = O(1), there exists an absolute constant c > 0 such that
σ1 (Π) ≤ c σK (Π) = O
p
n/K .
First, the three conditions in Equation (2) become:
8 2 np ≫ µ2 κ8 κΠ K log4 d = O(1) · K 2 log4 d, 8 κ8 κΠ K log2 d = O(1) · K log2 d, β 4 µ1/3 κ8/3 κΠ K σ2/3 n≫ 1 (Π). β2/3
p≫
2/3 (n/K)1/3, the right-hand side is bounded above by an absolute constant times In the last line, using σ2/3 1 (Π) ≤ c
K · (n/K)1/3 = K 2/3 n1/3 . Hence the condition n ≫ K 2/3 n1/3 is equivalent to n2/3 ≫ K 2/3 , i.e. n ≫ K. Absorbing all universal constants into the ≫ notation yields np ≫ K 2 log4 d,
p ≫ K log2 d,
Second, the two conditions in Equation ((3)) are: p 5 1/2 ηκ7/2 κΠ µ Kσ1 (Π) log d , ∆≫ √ β1/2 n 26
n ≫ K.
∆≫
2 1/4 ηκ2 κΠ µ 1/2 β
q p p σ1 (Π) K/n p/n 1/4 K log d.
√ Plugging the boundedness assumptions and σ1 (Π) ≤ c n/K into the first inequality gives ∆ ≫ O(1) · so ∆ ≫
p
K·
p √ p n/K · log d = O(1) · K log d, √ n
K log d. For the second inequality, note that p p p σ1 (Π) K/n ≤ c n/K · K/n = c = O(1),
hence the square root factor
q
√ σ1 (Π) K/n is also O(1). Consequently, we have p ∆ ≫ O(1) · (p/n)1/4 K log d.
Combining the two lower bounds, we obtain the more stringent requirement
∆≫
p
o n K log d max 1, (p/n)1/4 .
All omitted constants are absolute and independent of n, p, K, d, η, ∆. Therefore, under the stated boundedness assumptions, the conditions of Theorem 1 imply the simplified conditions (4)–(5), and the conclusion of the theorem remains valid with probability at least 1 − O(d−10 ). This completes the proof. Appendix B. Technical Lemmas Lemma 2. Under MMSG, we have ∆ σK (P) ≥ √ σK (Π). 2κ Proof of Lemma 2. Because Π⊤ is row-full-rank, we can apply the multiplicative singular value inequality for matrices with full column/row rank: σK (ΘΠ⊤ ) ≥ σK (Θ) σK (Π⊤ ).
27
Now we bound σK (Θ) from below. For any distinct k, ℓ ∈ [K], we have ∆ ≤ kθk − θℓ k = kΘ(ek − eℓ )k ≤ σ1 (Θ) kek − eℓ k =
√ 2 σ1 (Θ),
√ √ hence σ1 (Θ) ≥ ∆/ 2. By definition κ = σ1 (Θ)/σK (Θ), we have σK (Θ) = σ1 (Θ)/κ ≥ ∆/( 2κ). Substituting this into the previous inequality gives ∆ σK (P) ≥ √ σK (Π), 2κ which completes the proof. Lemma 3. Under MMSG, let P = ΘΠ⊤ . Define κ = σ1 (Θ)/σK (Θ), κP = σ1 (P)/σK (P), and κΠ = σ1 (Π)/σK (Π). Then, we have
κP ≤ κκΠ .
(B.1)
Proof of Lemma 3. The pure-individual condition guarantees that Π contains a K × K identity submatrix; since K ≤ n by the model definition, Π has full column rank K. For any matrices A ∈ R p×K and B ∈ Rn×K with B having full column rank, the singular values satisfy σ1 (AB⊤ ) ≤ σ1 (A) σ1 (B),
(B.2)
σK (AB⊤ ) ≥ σK (A) σK (B).
(B.3)
Applying Equations (B.2) and (B.3) with A = Θ and B = Π yields
σ1 (P) ≤ σ1 (Θ) σ1 (Π),
(B.4)
σK (P) ≥ σK (Θ) σK (Π).
(B.5)
Dividing Equation (B.4) by Equation (B.5) gives σ1 (P) σ1 (Θ) σ1 (Π) ≤ · = κκΠ , σK (P) σK (Θ) σK (Π) which is exactly Equation (B.1). This completes the proof.
28
(B.6)
Lemma 4. Under MMSG, let U be the right singular vectors of P = ΘΠ⊤ , i.e., P = VΣU⊤ with U⊤ U = IK . Then max kUi,: k2 ≤ i∈[n]
1 , λK (Π⊤ Π)
where λK (·) denotes the smallest eigenvalue. Consequently, we have
µ1 :=
1 n max kUi,: k2 ≤ . K i∈[n] β
Proof of Lemma 4. The row vectors of U satisfy the same simplex structure as in the mixed membership stochastic blockmodel (MMSB) (Mao et al., 2021). Specifically, by Lemma 1, there exists an invertible matrix B such that U = ΠB and B = U(I, :), where I indexes one pure node per component. Then, exactly the same geometric argument as in Lemma 3.1 of (Mao et al., 2021) shows that
max kUi,: k2 = max(D−1 )kk ≤ λmax (D−1 ) = i∈[n]
k∈[K]
1 , λK (D)
where D = Π⊤ Π. Hence, the claimed bound holds. Substituting it into the definition of µ1 gives µ1 ≤ 1/β because β = σ2K (Π)/(n/K) = λK (Π⊤ Π) · (K/n). This completes the proof. References Abbe, E., Fan, c., & Wang, K. (2022). An lp theory of PCA and spectral clustering. Annals of Statistics, 50, 2359–2385. Aeberhard, S., Coomans, D., & De Vel, O. (1994). Comparative analysis of statistical pattern recognition methods in high dimensional settings. Pattern Recognition, 27, 1065–1077. Airoldi, E. M., Blei, D. M., Fienberg, S. E., & Xing, E. P. (2008). Mixed Membership Stochastic Blockmodel. Journal of Machine Learning Research, 9, 1981–2014. Balakrishnan, S., Wainwright, M. J., & Yu, B. (2017). Statistical guarantees for the EM algorithm: From population to sample-based analysis. Annals of Eugenics, 45, 77–120. Cai, C., Li, G., Chi, Y., Poor, H. V., & Chen, Y. (2021). Subspace estimation from unbalanced and incomplete data matrices: ℓ2,∞ statistical guarantees. Annals of Statistics, 49, 944 – 967. Chen, L., & Gu, Y. (2024). A spectral method for identifiable grade of membership analysis with binary responses. psychometrika, 89, 626–657. Chen, X., & Yang, Y. (2021). Cutoff for exact recovery of gaussian mixture models. IEEE Transactions on Information Theory, 67, 4223–4238. Chen, X., & Zhang, A. Y. (2024). Achieving optimal clustering in gaussian mixture models with anisotropic covariance structures. Advances in Neural Information Processing Systems, 37, 113698–113741. Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39, 1–22. Fisher, R. A. (1936). The use of multiple measurements in taxonomic problems. Annals of Eugenics, 7, 179–188. Fortunato, S., & Hric, D. (2016). Community detection in networks: A user guide. Physics Reports, 659, 1–44.
29
Gillis, N., & Vavasis, S. A. (2013). Fast and robust recursive algorithmsfor separable nonnegative matrix factorization. IEEE transactions on pattern analysis and machine intelligence, 36, 698–714. Güvenir, H. A., Demiröz, G., & Ilter, N. (1998). Learning differential diagnosis of erythemato-squamous diseases using voting feature intervals. Artificial Intelligence in Medicine, 13, 147–165. Holland, P. W., Laskey, K. B., & Leinhardt, S. (1983). Stochastic blockmodels: First steps. Social Networks, 5, 109–137. Jain, A. K. (2010). Data clustering: 50 years beyond K-means. Pattern Recognition Letters, 31, 651–666. Jana, S., Yang, K., & Kulkarni, S. (2025). Adversarially robust clustering with optimality guarantees. IEEE Transactions on Information Theory, . Jin, J., Ke, Z. T., & Luo, S. (2024). Mixed membership estimation for social networks. Journal of Econometrics, 239, 105369. Ke, Z. T., & Wang, M. (2024). Using SVD for topic modeling. Journal of the American Statistical Association, 119, 434–449. Li, Z. (2025). Exact recovery of community detection in k-community gaussian mixture models. European Journal of Applied Mathematics, 36, 491–523. Lloyd, S. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory, 28, 129–137. Löffler, M., Zhang, A. Y., & Zhou, H. H. (2021). Optimality of spectral clustering in the Gaussian mixture model. Annals of Statistics, 49, 2506–2530. Lu, Y., & Zhou, H. H. (2016). Statistical and computational guarantees of lloyd’s algorithm and its variants. arXiv preprint arXiv:1612.02099, . Mao, X., Sarkar, P., & Chakrabarti, D. (2021). Estimating mixed memberships with sharp eigenvector deviations. Journal of the American Statistical Association, 116, 1928–1940. McQueen, J. B. (1967). Some methods of classification and analysis of multivariate observations. In Proc. of 5th Berkeley Symposium on Math. Stat. and Prob. (pp. 281–297). Ndaoud, M. (2022). Sharp optimal recovery in the two component Gaussian mixture model. Annals of Statistics, 50, 2096–2126. Ng, A., Jordan, M., & Weiss, Y. (2001). On spectral clustering: Analysis and an algorithm. Advances in Neural Information Processing Systems, 14. Pearson, K. (1894). III. Contributions to the mathematical theory of evolution. Proceedings of the Royal Society of London, 54, 329–333. Qing, H. (2024a). Finding mixed memberships in categorical data. Information Sciences, 676, 120785. Qing, H. (2024b). Grade of Membership Analysis for Multi-Layer Ordinal Categorical Data. Statistica Sinica, 38. Srivastava, P. R., Sarkar, P., & Hanasusanto, G. A. (2023). A robust spectral clustering algorithm for sub-Gaussian mixture models with outliers. Operations Research, 71, 224–244. Von Luxburg, U. (2007). A tutorial on spectral clustering. Statistics and Computing, 17, 395–416. Zhang, A. Y., & Zhou, H. Y. (2024). Leave-one-out singular subspace perturbation analysis for spectral clustering. Annals of Statistics, 52, 2004–2033.
30