A Latent Semiparametric Mixture Model for Robust Heterogeneity Learning in High Dimensions

Introduction

High-dimensional data heterogeneity has received growing attention in fields such as genomics, precision medicine, and social science. Cluster analysis, a widely used unsupervised learning approach, aims to uncover such heterogeneity by partitioning populations into subgroups. Popular methods include $K$-means [20], which minimizes within-cluster variation; hierarchical clustering, which recursively splits data based on group-wise dissimilarity; and spectral clustering (Lu et al., 2016), which uses the graph Laplacian and its eigen-decomposition to detect heterogeneity. These approaches are largely algorithmic. In contrast, model-based clustering [11] adopts a probabilistic framework, modeling data as drawn from a finite mixture of distributions, with each component representing a cluster. A common assumption is that the components follow Gaussian or some specific parametric distributions [33], which simplifies computation and analysis but can be overly restrictive in practice.

A key motivation for learning data heterogeneity, involving complex distributions not well captured by simple parametric models, is clustering analysis of single-cell RNA sequencing (scRNA-seq) data, where the primary goal is to uncover biological heterogeneity by identifying unknown cell types [2]. Accurate clustering is crucial for reliable downstream analyses, such as detecting highly or spatially variable genes [39, 49]. In this paper, we introduce a novel method, Clustering with LatEnt Semiparametric MixturE (CESME), which shows strong empirical performance closely aligned with its theoretical guarantees. As a demonstration, Figure 1 presents a two-dimensional Uniform Manifold Approximation and Projection (UMAP; [32]) plot of clustering results on a benchmark scRNA-seq dataset ([5]). Specifically, we compare results from CESME, Seurat ([18]), and a well-established sparse $K$-means method (SK-means; [47]), using the true cell labels as reference. CESME clearly outperforms the others, particularly Seurat, a widely used scRNA-seq analysis pipeline, highlighting its effectiveness in clustering high-dimensional data with complex distributions and feature dependencies.

In Figure 2, we showcase the key idea behind CESME. Panel (a) displays an observed dataset with $n = 90$ observations, where the clustering boundary is difficult to discern. We then apply a coordinate-wise cubic-root transformation to $X_1$ and $X_2$ to obtain the transformed features $X'_1$ and $X'_2$, and plot the latent data based on these new features. As shown in Panel (b), the latent data become more clearly separated, which in turn makes the clustering boundary much easier to visualize in Panel (c). This example illustrates the key idea in CESME: we assume that observed data exhibit clear cluster patterns in the latent space via a coordinate-wise transformation; so a seemingly difficult clustering problem in the observation space is actually much easier in the latent space.

Vector (PDF)

Formulation and contributions

To formally introduce CESME, we consider the following data-generating model for the observed data and the clustering problem. First, let $\mathbf{z}^{*} = (z_{1}^{*}, \ldots, z_{n}^{*})$ denote the cluster membership variables, where for $i = 1, \ldots, n$, $z_{i}^{*}$'s are independent and identically distributed with $\mathbb{P}(z_{i}^{*} = k) = \lambda_{k}^{*}$, and $\sum_{k} \lambda_{k}^{*} = 1$. Given $z_{i}^{*} = k$, a latent variable $\mathbf{X}_{i}^{\prime}$ is generated from a multivariate normal distribution, $\mathbf{X}_{i}^{\prime} \mid z_{i}^{*} = k \sim\mathcal{N}(\boldsymbol{\mu}_{k}^{*}, \boldsymbol{\Sigma}^{*})$.

Let $\mathbf{h} := (h_{1}, \ldots, h_{p})$ denote a collection of $p$ strictly increasing functions. The observed variable $\mathbf{X}_{i}$ is obtained via applying $\mathbf{h}$ coordinate-wise to latent $\mathbf{X}_{i}^{\prime}$, that is, $\mathbf{X}_{i} = \mathbf{h}(\mathbf{X}_{i}^{\prime}) = (h_{1}(X_{i1}^{\prime}), \ldots, h_{p}(X_{ip}^{\prime}))^{\mathsf{T}}$. For $\ell= 1, \ldots, p$, let $g_{\ell}^{*}$ denote the inverse of $h_{\ell}$, and write $\mathbf{g}^{*} := (g_{1}^{*}, \ldots, g_{p}^{*})$. Then, $X_{i\ell}^{\prime} = g_{\ell}^{*}(X_{i\ell})$ for $\ell= 1, \ldots, p$. Joint identifiability of the latent transformations $\mathbf{g}^{*}$ and mixture parameters $\boldsymbol{\beta}^{*}$ remains open for general $K$ and covariance structures. Section B of the supplementary material discusses partial identifiability and establishes identifiability in several structured settings. The clustering guarantee in Theorem 3.1 is formulated relative to a fixed data-generating representation: under its stated assumptions, the theorem bounds the excess clustering risk of CESME relative to the corresponding oracle Bayes rule. This guarantee does not require global identifiability and should not be interpreted as establishing uniqueness of the latent parameterization. We also assume that, for some $\delta> 0$, $\max_{1 \le\ell\le p} \mathbb{E}|X_{i\ell}|^{4+\delta} \le M_{4+\delta} < \infty$, which provides sufficient condition on moment and tail behavior of the observed data $\mathbf{X}_i$, and enables analyzing the empirical $cdf$ of $\mathbf{X}$.

Note the optimal clustering (Bayes rule) maximizes posterior probability [33]. In our generative model, we see that $\mathbb{P}(z_i^* \mid\mathbf{X}_i) = \mathbb{P}(z_i^* \mid\mathbf{X}')$, which means that the Bayes rule can be calculated in the latent space. Based on the Gaussian assumption, the posterior maximization reduces to minimizing the quadratic loss $[\{\mathbf{g}^*(\mathbf{X}_i)-\boldsymbol{\mu}_k^*\}^{\mathsf{T}}(\boldsymbol{\Sigma}^*)^{-1}\{\mathbf{g}^*(\mathbf{X}_i)-\boldsymbol{\mu}_k^*\}/2-\log(\lambda_k^*)]$. Note that we assume a common covariance matrix, which enables a simplified form of the optimal decision rule

where the discriminant vectors are $\boldsymbol{\gamma}_k^*=(\boldsymbol{\Sigma}^*)^{-1}(\boldsymbol{\mu}_1^*-\boldsymbol{\mu}_k^*)$ for $k=2,\ldots,K$, and $\boldsymbol{\gamma}_1^*$ is defined to be zero. The linearity of the optimal rule permits an efficient and optimal clustering in high-dimensional spaces when we further assume that each $\boldsymbol{\gamma}_k^*$ is $s$-sparse, which has been well demonstrated in Cai et al. [7]. Another important observation is that we assume the covariance matrices are equal across clusters only in the transformed space. With nonlinear transformations, which is the case of primary interest here, the covariance matrices in the original space can actually be cluster-dependent. This also distinguishes our approach from that of Cai et al. [7] (CHIME). It would be of interest to further allow cluster-dependent covariance after transformation, adding more flexibility to our framework. When $p$ is small, this modeling choice does not substantially increase the complexity of the learning algorithm. In high dimensions, however, one would need efficient algorithms to learn the resulting quadratic Bayes rule. Under additional structural assumptions on the covariances, this may be achieved by incorporating large covariance learning algorithms into our procedure, such as the graphical lasso under a sparse precision matrix assumption [12].

One may wonder why not use cluster-dependent transformations to generate the observed data. We prefer cluster-independent transformations for the sake of interpretability, because in (1), $X_{i\ell}$ contributes to the optimal clustering if and only if $g_\ell^*(X_{i\ell})$ has contributions via the nonzero elements of the discriminant vector. If we choose cluster-ependent transformations dependent transformations, then the optimal clustering compares $\left[\{\mathbf{g}_k^*(X_i)-\boldsymbol{\mu}_k^*\}^{\mathsf{T}}(\boldsymbol{\Sigma}^*)^{-1}\{\mathbf{g}_k^*(X_i)-\boldsymbol{\mu}_k^*\}/2-\log(\lambda_k^*)\right]$. Thus, we lose the clear interpretability and the inference in high dimensions becomes much more complicated to handle.

Based on the above discussions, we prefer cluster-independent covariance matrices and transformations to strike an excellent balance between model flexibility, interpretability, and learning feasibility which is particularly important in high dimensions. Our modeling strategy has been proven to be practically useful in the benchmark data study in § 5. In fact, our method delivers the best clustering results across six scRNA-seq datasets. Our method is expected to be useful in other applications such as imaging, flow cytometry, metabolomics, and financial data, where original features have highly skewed or heavy-tailed distributions.

Our model is also related to the nonparametric mixture models studied by Hall and Zhou [16]. When $\boldsymbol{\Sigma}^*=\mathbf{I}$, the transformed coordinates $g_1^*(X_{i1}),\ldots,g_p^*(X_{ip})$ are conditionally independent given $z_i^*=k$. Thus, the conditional density factorizes as

Hence, our model includes the conditionally independent nonparametric mixture models in Hall and Zhou [16] as a special case. However, we believe that it is important to account for within-cluster dependence among features for the success of our method. We therefore do not impose a pre-specified $\boldsymbol{\Sigma}^*$ and rely on a data-driven $\boldsymbol{\Sigma}^*$. The numerical experiments with nonparametric mixture models in § 4.1 also support our view.

Under the proposed data-generative model, we aim to approximate the optimal rule via the CESME algorithm. We adopt an alternating procedure that iteratively estimates the latent transformations and the parameters of the mixture of Gaussians. Given previous estimates, it alternates between updating mixture parameters via a modified expectation–maximization (EM) and estimating transformations using weighted empirical conditional cumulative distribution functions (CDFs). This strategy is computationally efficient in high dimensions.

To establish theoretical guarantees for CESME, we identify the population solution path, i.e., the population-level alternating algorithm, in parallel with the empirical path using data. The estimation error is then decomposed into the numerical error, representing the gap between the population path and true parameters, and the statistical error, capturing the deviation between the empirical and population paths. While our analysis is motivated by EM analysis for the Gaussian mixtures [3], the presence of latent transformations makes it challenging to establish gradient stability and restricted concavity of the auxiliary $Q$-function with respect to all model parameters, rendering standard EM analysis for the population path inapplicable. To address this, we identify the appropriate parameter spaces for both the mixture parameters and latent transformations, and carefully analyze the Taylor expansion and its remainders over these spaces.

For parameters in traditional Gaussian mixtures, the estimators are Lipschitz, allowing Rademacher-based analysis to directly yield finite sample error bounds between empirical and population paths. In contrast, the estimators of latent transformations in our model are composite functions of empirical CDFs and EM-estimated mixture parameters, making Rademacher complexity analysis intractable. Instead, we develop a fundamentally different framework by bounding the complexity via the Vapnik–Chervonenkis (VC) dimensions of the corresponding function classes. We show that for clustering $n$ subjects with $p$ features, after sufficiently many yet manageable iterations, the estimation error is of order $O_P\left[\left\{s\log(np)/n\right\}^{1/2}\right]$, with sparsity $s$ defined later. This establishes both convergence and estimation consistency of our algorithm. Moreover, we benchmark CESME’s clustering performance against the optimal Bayes rule, and show that they are nearly risk-equivalent, with excess risk of order $O(s\log(np)/n)$. If transformations were known, the excess-risk bound would be $O(s\log(p)/n)$ according to Cai et al. [7]. The bound with unknown transformations is $O(s\log(np)/n)$, which has the same order when $\log n = O(\log p)$, including the regime $p \gg n$. Thus, in this regime, estimating the marginal monotone transformations incurs no additional cost in the order of the excess-risk bound.

Related works

This paper aims to develop a practical, off-the-shelf model-based clustering method for high-dimensional data with complex distributions varying across dimensions and with feature dependencies. While the classical Gaussian mixture model is widely used [33, 7, 30], extensions include mixtures of parametric heavy-tailed distributions [27, 26] and nonparametric mixtures [16]. See Gormley et al. [13] for a comprehensive review. Parametric mixtures, like Gaussian ones, are often task-specific and not broadly applicable. Nonparametric mixtures offer greater flexibility, but often assume restrictive conditions, such as within-cluster feature independence, which can degrade empirical performance. These models also tend to struggle in high dimensions.

Several methods have been proposed for high-dimensional clustering. Although spectral clustering is popular for low-dimensional data, Zhang and Zhou [50] showed that it can be suboptimal for clustering under sub-Gaussianity. Wang et al. [46] extended it to high dimensions using $L_1$-regularized principal components, but, as shown in § 5, its performance is limited when real data depart substantially from Gaussianity. Pan and Shen [35] proposed a modified model-based approach with diagonal covariance and sparsity-inducing penalties for dimension reduction; however, by ignoring within-cluster dependence, it is less effective in more general settings. As noted earlier, the widely used sparse $K$-means imposes an $L_1$ penalty on the cluster means and has been extended to sparse convex clustering [44], although theoretical guarantees for both methods remain incomplete. Jin and Wang [23] proposed IF-PCA, which selects informative features using Kolmogorov–Smirnov tests and, in high dimensions, higher criticism, before applying PCA and $K$-means on the principal components. While effective for high-dimensional clustering, IF-PCA’s guarantees rely on Gaussianity or sub-Gaussianity. Cai et al. [7] introduced CHIME, which exploits sparsity in the optimal clustering rule for Gaussian mixtures and offers theoretical guarantees. However, it is closely tied to the Gaussian assumption and may perform poorly when the data are far from Gaussian.

Organization of the paper

In § 2 we present the details of CESME algorithm, with theoretical analysis given in § 3. § 4 provides simulations supporting the theoretical results and demonstrating CESME’s empirical performance. In § 5 we apply CESME to six scRNA-seq datasets and compare it with benchmark pipelines. Technical proofs and additional results appear in the online supplementary file.

An alternating algorithm and the model-based clustering

Throughout the following algorithm, we treat the number of clusters $K$ as known. We use nested loops to estimate the transformations $\mathbf{g}^{*}$ and mixture parameters $\boldsymbol{\beta}^{*} = (\boldsymbol{\lambda}^{*}, \boldsymbol{\mu}^{*}, \boldsymbol{\Sigma}^{*})$, where $\boldsymbol{\lambda}^{*} = (\lambda_{1}^{*},\ldots,\lambda_{K}^{*})^{\mathsf{T}}$ and $\boldsymbol{\mu}^{*} = ((\boldsymbol{\mu}_{1}^{*})^{\mathsf{T},\ldots,(\boldsymbol{\mu}_{K}^{*})^{\mathsf{T}})^{\mathsf{T}}}$. At iteration $t$, the outer loop alternates between updating $\widehat{\mathbf{g}}^{(t)}$ using $(\widehat{\boldsymbol{\beta}}^{(t-1)},\widehat{\mathbf{g}}^{(t-1)})$ and updating $\widehat{\boldsymbol{\beta}}^{(t)}$ via an EM algorithm using $(\widehat{\mathbf{g}}^{(t)},\widehat{\boldsymbol{\beta}}^{(t-1)})$, which is the inner loop.

Latent transformations via conditional CDF inversion

To fix a convenient coordinate-wise location-scale convention for the latent representation, we impose that $\sum_{k=1}^{K}\lambda_{k}^{*}\mu_{k\ell}^{*}=0$, and $(\sigma_{\ell}^{*})^{2}+\sum_{k=1}^{K}\lambda_{k}^{*}(\mu_{k\ell}^{*})^{2}=1$ for each $\ell=1,\ldots,p$. Thus, the marginal means and variances of $\mathbf{g}^{*}(\mathbf{X}_{i})$ are zero and one, respectively. Otherwise, we can impose a shift-scale transformation on top of the $\mathbf{g}^{*}$ functions to achieve that.

Remark 1. The above centering and scaling conditions specify a coordinate-wise location-scale convention, achievable by reparameterization, but they do not establish global identifiability. Identifiability in the general case remains open, but is not required for the clustering guarantee in Theorem 3.1, interpreted as discussed in § 1.1. See Section B of the supplementary material for detailed identifiability results.

By the monotonicity of $g_{\ell}^{*}$, we have $\mathbb{P}(X_{i\ell}\leq x\mid z_{i}^{*}=k)=\mathbb{P}\{g_{\ell}^{*}(X_{i\ell})\leq g_{\ell}^{*}(x)\mid z_{i}^{*}=k\}=$ $\Phi[(\sigma_\ell^*)^{-1}\{g_\ell^*(x)-\mu_{k\ell}^*\}]$ for any $x$ in the support of $X_{i\ell}$, with the standard normal CDF $\Phi$. Using $\sum_{k=1}^K \lambda_{k\ell}^*\mu_{k\ell}^*=0$, this yields the following natural population representation of $g_\ell^*$:

The above equation motivates estimating the latent transformation as a weighted combination of cluster-specific inverse normal CDF transforms, where each conditional CDF is replaced by its empirical estimate.

Given the update $\widehat{\boldsymbol{\beta}}^{(t-1)}$ from the $(t-1)$th iteration of the alternating maximization algorithm, we estimate $g_\ell^*$ by replacing $\lambda_{k\ell}^*$ and $\sigma_\ell^*$ in (2) with $\widehat{\lambda}_{k\ell}^{(t-1)}$ and $\widehat{\sigma}_\ell^{(t-1)}$ respectively, where $\widehat{\sigma}_\ell^{(t-1)}>0$ is the square root of the $\ell$th diagonal entry of $\widehat{\boldsymbol{\Sigma}}^{(t-1)}$. Then, we can estimate $\mathbb{P}(X_{i\ell}\le x\mid z_i^*=k)$ by $\widehat{\psi}_{\ell k}(x)=\{\sum_{i=1}^n\mathbb{I}(\widehat{z}_i^{(t-1)}=k)\}^{-1}\sum_{i=1}^n\mathbb{I}(X_{i\ell}\le x,\widehat{z}_i^{(t-1)}=k)$, where $\widehat{z}_i^{(t-1)}$ is the cluster assignment from the $(t-1)$th iteration, with its explicit form given in § 2.2. To sum up, we have the update

Since $\widehat{g}^{(t)}$ depends on $\widehat{z}^{(t-1)}$ and $\widehat{\boldsymbol{\beta}}^{(t-1)}$, and $\widehat{z}^{(t-1)}$ in turn depends on both $\widehat{\boldsymbol{\beta}}^{(t-1)}$ and $\widehat{g}^{(t-1)}$, it follows that $\widehat{g}^{(t)}$ is ultimately a function of $\widehat{\boldsymbol{\beta}}^{(t-1)}$ and $\widehat{g}^{(t-1)}$.

In practice, as the empirical CDF can be zero or one which should be avoided in the inverse normal CDF, we consider truncating $\widehat{\psi}_{\ell k}(x)$ as $\widetilde{\psi}_{\ell k}(x)=\widehat{\psi}_{\ell k}(x)\mathbb{I}\{\delta_n\le\widehat{\psi}_{\ell k}(x)\le1-\delta_n\}+(1-\delta_n)\mathbb{I}\{\widehat{\psi}_{\ell k}(x)>1-\delta_n\}+\delta_n\mathbb{I}\{\widehat{\psi}_{\ell k}(x)<\delta_n\text{ or }\sum_{i=1}^n\mathbb{I}(\widehat{z}_i^{(t-1)}=k)=0\}$ and replacing $\widehat{\psi}_{\ell k}(x)$ by $\widetilde{\psi}_{\ell k}(x)$ in (3), where $\delta_n=1/(4n^{1/4}\sqrt{\pi\log n})$ is similarly set as in Lafferty et al. [25] for estimating the nonparanormal graphical model.

Modified EM iterations

Given $\widehat{g}^{(t)}$ and $\widehat{\boldsymbol{\beta}}^{(t-1)}$, the update $\widehat{\boldsymbol{\beta}}^{(t)}$ is obtained via an EM algorithm, which serves as the inner loop of our algorithm. Let $\widehat{\boldsymbol{\beta}}^{(t,r)}$ denote the update at the $r$th iteration of the inner loop, with $r=1,2,\ldots$, and $\widehat{\boldsymbol{\beta}}^{(t,0)}=\widehat{\boldsymbol{\beta}}^{(t-1)}$. In the $E$-step, since $\gamma_k^*$’s are assumed to be $s$-sparse, we estimate them by

for $k=2,\ldots,K$, and set $\widehat{\gamma}_1^{(t,r)}(\rho)=0$. We use the lasso penalty because it has been used in sparse $K$-means [47] and CHIME [7]. This step differs from the canonical EM algorithm for updating the discriminant vectors.

Given $\widehat{\gamma}_k^{(t,r)}(\rho)$, the cluster assignment for the $i$th observation is updated as $\widehat{z}_i^{(t,r)}(\rho)=\arg\min_k\left[\left\{\widehat{\mathbf{g}}^{(t)}(\mathbf{X}_i)-\left(\widehat{\mu}_k^{(t,r-1)}+\widehat{\mu}_1^{(t,r-1)}\right)/2\right\}^{\top}\widehat{\gamma}_k^{(t,r)}(\rho)-\log\left(\widehat{\lambda}_k^{(t,r-1)}\right)\right]$ by following the Bayes rule. For each $k$, the posterior probability is updated as

Then, in the $M$-step, the conditional expectation of the log-likelihood with respect to the distribution of $z_i|\mathbf{X}_i$, with current estimates $\widehat{\mathbf{g}}^{(t)}$ and $\widehat{\beta}^{(t,r-1)}$, is

Maximizing $Q_{\widehat{\mathbf{g}}^{(t)}}(\beta|\widehat{\beta}^{(t,r-1)})$ with respect to $\beta$ yields the update $\widehat{\beta}^{(t,r)}=\arg\max_\beta Q_{\widehat{\mathbf{g}}^{(t)}}(\beta|\widehat{\beta}^{(t,r-1)})$, where

For a pre-specified threshold $\varepsilon_0$ and small $\epsilon>0$, we define the relative discrepancy

and stopping criterion

In the $t$th iteration of the outer loop, the EM algorithm alternates between the $E$- and $M$-steps, and terminates at iteration $r$ if $\operatorname{RD}(\widehat{\beta}^{(t,r)},\widehat{\beta}^{(t,r-1)})<\varepsilon_0$. Once the criterion is met, we set $\widehat{\beta}^{(t)}=\widehat{\beta}^{(t,r-1)}$ and $\widehat{\mathbf{z}}^{(t)}=\widehat{\mathbf{z}}^{(t,r-1)}$. The outer loop terminates at iteration $t$ if $\operatorname{RD}(\widehat{\boldsymbol{\beta}}^{(t)},\widehat{\boldsymbol{\beta}}^{(t-1)}) < \varepsilon_0$. Our procedure is summarized in Algorithm 1.

Algorithm 1 High-dimensional CESME

Input: $p$ by $n$ data matrix $\mathbf{X}=(\mathbf{X}_1,\ldots,\mathbf{X}_n)$ and prespecified number of clusters $K$.

1. Initialize: $\widehat{\boldsymbol{\beta}}^{(0)}=(\widehat{\boldsymbol{\lambda}}^{(0)},\widehat{\boldsymbol{\mu}}^{(0)},\widehat{\boldsymbol{\Sigma}}^{(0)})$ and $\widehat{\mathbf{g}}^{(0)}$ (see Section 3.4 for a provable initialization).

2. Compute $\widehat{\gamma}_k^{(0)}$ and $\widehat{\mathbf{z}}^{(0)}$.

3. for $t=1,2,\ldots$ do {Outer Loop}

4. Compute $\widehat{\mathbf{g}}^{(t)}$ based on $\widehat{\mathbf{z}}^{(t-1)}$ and $\widehat{\boldsymbol{\beta}}^{(t-1)}$.

5. Let $\widehat{\boldsymbol{\beta}}^{(t,0)}=\widehat{\boldsymbol{\beta}}^{(t-1)}$.

6. for $r=1,2,\ldots$ do {Inner Loop}

7. The $E$-step: compute $\widehat{\gamma}_k^{(t,r)}$, $\widehat{w}_k^{(t,r)}$, and $\widehat{\mathbf{z}}^{(t,r)}$ based on $\widehat{\mathbf{g}}^{(t)}$ and $\widehat{\boldsymbol{\beta}}^{(t,r-1)}$.

8. The $M$-step: compute $\widehat{\boldsymbol{\beta}}^{(t,r)}$ based on $\widehat{\boldsymbol{\gamma}}^{(t,r)}$, $\widehat{\mathbf{g}}^{(t)}$, and $\widehat{\boldsymbol{\beta}}^{(t,r-1)}$.

9. if condition (2.4) is satisfied for $\widehat{\boldsymbol{\beta}}^{(t,r)}$ and $\widehat{\boldsymbol{\beta}}^{(t,r-1)}$ then

10. Terminate the inner loop and output $\widehat{\boldsymbol{\beta}}^{(t)}=\widehat{\boldsymbol{\beta}}^{(t,r-1)}$ and $\widehat{\mathbf{z}}^{(t)}=\widehat{\mathbf{z}}^{(t,r-1)}$.

11. if condition (2.4) is satisfied for $\widehat{\boldsymbol{\beta}}^{(t)}$ and $\widehat{\boldsymbol{\beta}}^{(t-1)}$ then

12. Terminate the outer loop and output $\widehat{\mathbf{z}}=\widehat{\mathbf{z}}^{(t)}$.

Output: Clustering assignments $\widehat{\mathbf{z}}$.

Whenever the stopping criterion (2.4) is met for $\widehat{\boldsymbol{\beta}}^{(T)}$ and $\widehat{\boldsymbol{\beta}}^{(T-1)}$ at the $T$th outer iteration, our algorithm terminates with $\widehat{\mathbf{z}}=\widehat{\mathbf{z}}^{(T)}$ as the estimated cluster assignment, and $\widehat{\boldsymbol{\theta}}=(\widehat{\boldsymbol{\beta}},\widehat{\mathbf{g}})=(\widehat{\boldsymbol{\beta}}^{(T)},\widehat{\mathbf{g}}^{(T)})$ as the final estimates of the parameters $\boldsymbol{\beta}$ of the mixture and transformations $\mathbf{g}$, respectively.

As a practical alternative to the nested two-loop scheme, one may use a one-step alternating schedule: perform a single EM update for $\boldsymbol{\beta}$ followed by one transform update for $\mathbf{g}$, and repeat. This one-step variant is straightforward to implement and often reduces per-iteration cost, though it typically requires more outer iterations to reach a comparable fit. For numerical comparisons between the one-step EM variant and Algorithm 1, see Section F.1 in the supplementary file.

Implementation details

The theoretical analysis uses a penalty parameter $\rho$ of order $\sqrt{\log(np)/n}$; see the proof of Theorem 3.1 in the supplementary file. For an estimate $\widehat{\boldsymbol{\theta}}(\rho)$ obtained with penalty $\rho$, we define the data-driven BIC-type criterion $\operatorname{Cr}(\rho)=2l_K(\widehat{\boldsymbol{\theta}}(\rho);\mathbf{X})-\{b_K+s(\rho)\}\log(p)$, where $l_K(\boldsymbol{\theta};\mathbf{X})=\sum_{i=1}^{n}\log\left(\sum_{k=1}^{K}\lambda_k\phi(\mathbf{g}(\mathbf{X}_i)\mid\boldsymbol{\mu}_k,\boldsymbol{\Sigma})\right)$, $b_K=K-1+p(K-1)+p(p-1)/2 is the baseline mixture parameter count, and $\tilde{s}(\rho)=\sum_{k=2}^{K}\sum_{j=1}^{p}\mathcal{I}(\tilde{\gamma}_{kj}(\rho)\ne0)$ is the effective sparsity.

At outer iteration $t\ge1$, we update the penalty as $\rho^{(t)}=\nu\rho^{(t-1)}+C_{\rho}^{(t)}\sqrt{\log(np)/n}$, where $\nu\in(0,1)$ is a fixed numerical decay parameter. At the first inner iteration, we select $C_{\rho}^{(t)}$ from a finite grid $\mathcal{C}$ of positive multipliers by maximizing $\operatorname{Cr}$ evaluated at the corresponding penalty values and estimates. The selected penalty $\rho^{(t)}$ remains fixed for the remaining inner iterations. § 4.1 specifies the initial penalty $\rho^{(0)}$, the value of $\nu$, and the candidate grid.

To enforce the zero mean and unit variance condition in practice, we standardize $\hat{g}(\mathbf{X}_1),\ldots,\hat{g}(\mathbf{X}_n)$. Specifically, let $\operatorname{sd}(\hat{g}_{\ell})=\left[(n-1)^{-1}\sum_{i=1}^{n}\{\hat{g}_{\ell}(X_{i\ell})-n^{-1}\sum_{i=1}^{n}\hat{g}_{\ell}(X_{i\ell})\}^{2}\right]^{1/2}$, and define the standardized transformation as $\tilde{g}_{\ell}(x)=\{\hat{g}_{\ell}(x)-n^{-1}\sum_{i=1}^{n}\hat{g}_{\ell}(X_{i\ell})\}/\{\operatorname{sd}(\hat{g}_{\ell})+\varepsilon_{n,p}\}$, where $\varepsilon_{n,p}=C\sqrt{\log(np)/n}$ with constant $C>0$ independent of $n$ and $p$. This small additive term ensures numerical stability when $\operatorname{sd}(\hat{g}_{\ell})$ is close to zero.

We have thus far assumed that the number of clusters $K$ is known. To select $K$ in practice, we run Algorithm 1 for each candidate $K$ and compute the corresponding $\operatorname{Cr}_{K}=\max_{\rho}\operatorname{Cr}(\rho)$ using estimates $\hat{\theta}(K)$. The selected number of clusters is the one that maximizes $\operatorname{Cr}_{K}$: $\hat{K}=\arg\max_{K}\operatorname{Cr}_{K}$. This approach is widely used in model-based clustering [13]. § 4 and § 5 demonstrate that CESME combined with the data-driven BIC-type criterion effectively selects the number of clusters, outperforming competing methods in both simulations and real data applications. We emphasize that this construction is heuristic and effective in practice; accordingly, we refer to $\operatorname{Cr}_{K}$ as a data-driven BIC-type criterion.

Statistical Analysis

Conditions and parameter space

To establish statistical guarantees, we first define $\|\theta\|_{s}^{2}=\|\beta\|_{2,s}^{2}+\|g\|_{\infty,s}^{2}$, where $\|\beta\|_{2,s}^{2}=\|\lambda\|_{2}^{2}+\sum_{k=1}^{K}\|\mu_{k}\|_{2,s}^{2}+\|\Sigma\|_{2,s}^{2}$. For a vector $\mathbf{a}$, matrix $\mathbf{M}$, and function $\mathbf{g}$, define $\|\mathbf{a}\|_{2,s}=\sup_{\mathbf{x}\in\Gamma(s),\|\mathbf{x}\|_{2}=1}|\mathbf{a}^{\mathsf{T}}\mathbf{x}|$, $\|\mathbf{M}\|_{2,s}=\sup_{\mathbf{x}\in\Gamma(s),\|\mathbf{x}\|_{2}=1}\|\mathbf{M}\mathbf{x}\|_{2}$. Here, $\Gamma(s)={\mathbf{a}\in\mathbb{R}^{p}:2|\mathbf{a}_{S^{c}}|_{1}\le $4\lVert a_S\rVert_1 + 3\sqrt{s}\lVert a\rVert_2$, for $S \subset\{1,\ldots,p\}$ and $|S|=s\}$ denotes a cone of sparse vectors, where $a_S$ is the subvector of $a$ supported on $S$. For $\mathbf{x}=(x_1,\ldots,x_p)^\mathsf{T}$, write $\mathbf{g}(\mathbf{x})=\{g_1(x_1),\ldots,g_p(x_p)\}^\mathsf{T}$. For a deterministic sequence $a_{n,p}\in(0,1/2)$, define $D_{n,p}=\{\mathbf{x}\in\mathbb{R}^p:a_{n,p}\leq\Phi((g_\ell^*(x_\ell)-\mu_{k\ell}^*)/\sigma_\ell^*)\leq1-a_{n,p}$ for each $k=1,\ldots,K$ and $\ell=1,\ldots,p\}$.

The sequence $a_{n,p}$ is chosen so that $D_{n,p}$ is nonempty; further requirements are specified in Lemma C.4 and the subsequent analysis. For functions $\mathbf{g}$ bounded on $D_{n,p}$, define the restricted supremum seminorms $\lVert\mathbf{g}\rVert_\infty=\max_{1\leq\ell\leq p}\sup_{\mathbf{x}\in D_{n,p}}|g_\ell(x_\ell)|$ and $\lVert\mathbf{g}\rVert_{\infty,s}=\sup_{\mathbf{a}\in\Gamma(s),\lVert\mathbf{a}\rVert_2=1}\sup_{\mathbf{x}\in D_{n,p}}|\mathbf{a}^\mathsf{T}\mathbf{g}(\mathbf{x})|$. The norm $\lVert\cdot\rVert_{2,s}$ and semi-norm $\lVert\cdot\rVert_{\infty,s}$ will provide refined measures of difference between estimated and true parameters, accounting for the sparsity. Finally, we define $R_{\min}=\min_{1\leq j<k\leq K}\lVert(\boldsymbol{\Sigma}^*)^{-1/2}(\boldsymbol{\mu}_k^*-\boldsymbol{\mu}_j^*)\rVert_2$ to quantify the separability between clusters, and impose the following conditions on the model.

Assumptions similar to Condition 3.1 have been considered for discriminant analysis [6, 31] and high-dimensional clustering [7]. The rate of $s$ relative to $n^{1/2}$ has also been used in high-dimensional clustering with sparse cluster means [17], ensuring the consistent estimation of $s$-sparse discriminant vectors [24].

Condition 3.2 requires a reasonable separation between clusters. Similar assumptions have appeared in [3] and [7]. Notably, there exist constants $\overline{M}$ and $\underline{m}$, independent of $n$, $p$, and $s$, such that $M\leq\overline{M}$ and $m\geq\underline{m}$. Under the mixing-weight and covariance bounds in Condition 3.2, its separation requirement holds whenever $R_{\min}>4\{\log_+(8\overline{M}/\underline{m})\}^{1/2}$. Since $\overline{M}$ and $\underline{m}$ are independent of $p$ and $s$, this gives a sufficient separation requirement independent of these dimensions; see the discussion following Lemma C.2 in Section A of the supplement. Since our goal is to benchmark clustering performance against the optimal Bayes rule, this requirement on $R_{\min}$ is milder than those in works that focus on recovering true labels [50, 29], where $R_{\min}$ is typically assumed to diverge with $n$ and $p$. Particularly, this condition holds when the non-zero entries of $(\boldsymbol{\Sigma}^{*})^{-1}(\boldsymbol{\mu}_{k}^{*}-\boldsymbol{\mu}_{j}^{*})$ are sufficiently large for $j \ne k$. Also, since our alternating maximization algorithm solves a non-convex problem, initialization is necessary as shown in Algorithm 1; we thus impose the following condition to ensure algorithmic convergence.

Condition 3.3 only requires the initial estimates to lie within a local basin of attraction whose radius is of order $R_{\min}$, rather than to be very close to the true parameters. Since Condition 3.2 places a lower bound on $R_{\min}$, this requirement does not shrink with $n$ and hence does not require the initialization to be consistent. This makes our initialization condition more attainable in practice than those imposed in prior works [7, 3, 17, 29, 53, 45]. For example, $\lVert\widehat{\boldsymbol{\beta}}^{(0)}-\boldsymbol{\beta}^{*}\rVert_{2,s} \lesssim s^{-1}$ is assumed in Cai et al. [7], which becomes more restrictive as $s$ diverges, and $\lVert\widehat{\boldsymbol{\mu}}_{k}^{(0)}-\boldsymbol{\mu}_{k}^{*}\rVert_{2} \le C\max_{k'\ne k}\lVert\widehat{\boldsymbol{\mu}}_{k'}^{(0)}-\boldsymbol{\mu}_{k'}^{*}\rVert_{2}$ for some constant $C>0$ is assumed in [53] and [45].

Guarantees on the estimations

The estimation guarantees of Algorithm 1, motivated by the traditional EM analysis, rely on decomposing the estimation error into statistical and numerical parts. By identifying the population path $\boldsymbol{\beta}^{(t,r)}$ and establishing its local concavity and gradient stability (Lemmas C.1 and C.2 in the supplementary file), we show that the numerical error decays geometrically over iterations. Since $\widehat{\boldsymbol{g}}^{(t)}$ in (3) is a composite function of empirical CDFs and earlier estimates $\widehat{\boldsymbol{\beta}}^{(t)}$, we control the statistical error by bounding the VC dimensions, instead of the intractable Rademacher complexity, of their function classes. In particular, the mixture parameter estimators are characterized by projections onto sparse sets with VC dimension $O(Ks)$, and the complexity of $\widehat{g}^{(t)}$ is controlled via Donsker preservation results [8]; see Lemma C.5 in the supplementary file. These arguments yield a statistical error of order $O(\sqrt{s\log(np)/n})$, while the numerical error decays geometrically and reaches the same order after $T_0$ iterations. Likewise, the number of inner-loop iterations is not fixed a priori, and $R_t$ EM updates in the $t$th outer loop suffice to reduce the optimization error below the statistical error. Consequently, the estimators satisfy

with high probability. See Theorems A.1 and A.2 in the supplementary file.

Clustering guarantees

To study the clustering guarantees of CESME, we first consider the optimal Bayes rule $\widehat{\boldsymbol{z}}_{\mathrm{opt}}=(\widehat{z}_{\mathrm{opt},1},\ldots,\widehat{z}_{\mathrm{opt},n})^{\mathsf{T}}$, where $\widehat{z}_{\mathrm{opt},i}=\arg\min_k\{(\widehat{\boldsymbol{\gamma}}_k^*)^{\mathsf{T}}[g^*(\mathbf{X}_i)-(\boldsymbol{\mu}_1^*+\boldsymbol{\mu}_k^*)/2]-\log(\lambda_k^*)\}$, which leverages oracle information. Following [7], we consider $\mathcal{L}(G)=\min_{\Pi\in\mathcal{S}_K}\mathbb{E}[\mathbb{I}\{\Pi(G(\mathbf{X}_{\mathrm{new}}))\ne z_{\mathrm{new}}^*\}]$, where $G$ is a generic clustering rule, $\mathcal{S}_K$ is the set of all label permutations of $\{1,\ldots,K\}$, and $z_{\mathrm{new}}^*$ is the true latent label of new observation $\mathbf{X}_{\mathrm{new}}$. As noted by [7], benchmarking against $\widehat{\boldsymbol{z}}_{\mathrm{opt}}$ allows to evaluate the CESME clustering performance objectively, so we compare $\mathcal{L}(\widehat{\boldsymbol{z}})$ and $\mathcal{L}(\widehat{\boldsymbol{z}}_{\mathrm{opt}})$. We prove that CESME is nearly risk-equivalent to the Bayes rule, as shown in Theorem 3.1.

Theorem 3.1 uses prescribed iteration numbers $T_0$ and $R_t$ to control finite-iteration errors. Since the requirements on $T_0$ and $R_t$ depend on unknown population quantities, we use the relative-change criterion (5) to terminate both loops in practice. Under Conditions 3.1-3.3, if $s\log(np)/n\to0$ and the iteration counts satisfy the stated requirements, the excess clustering risk of CESME relative to the oracle Bayes rule converges to zero, without requiring knowledge of the transformations $\mathbf{g}^{*}$. Theorem 3.1 is directly comparable to Theorem 3.2 in Cai et al. [7] where the bound for CHIME is $s\log(p)/n$.

Note that if we used the true transformations in CESME, then its error bound would be that of CHIME. For the $p\gg n$ case, our upper bound matches the Gaussian-mixture benchmark rate asymptotically because $\log(np)\leq2\log(p)$. Therefore, Theorem 3.1 shows that CESME attains the same rate as its oracle counterpart that has access to the true transformations. This theoretical conclusion is also supported by our simulation results.

Theorem 3.3 of Cai et al. [7] establishes a minimax lower bound of order $s\log(p)/n$ over the Gaussian-mixture parameter class considered there, and CHIME attains this rate. This provides a benchmark for our excess-risk bound. When $\log n=O(\log p)$, Theorem 3.1 gives an excess-risk bound for CESME of the same order, $O(s\log(p)/n)$.

Spectral-based initialization

Initialization plays an important role in fitting mixture models. For CHIME under Gaussian mixtures, Cai et al. [7] employed the HP algorithm (Hardt and Price, [19]) when $K=2$. For two-component Gaussian mixtures, Wu and Zhou (2021) showed that a random initialization is sufficient for estimation. However, such a strategy becomes unreliable when the number of components exceeds two. Indeed, Jin et al. [22] showed that, for Gaussian mixtures with $K\geq3$, random initialization can lead the EM algorithm to undesirable local optima with non-negligible probability. Alternatively, the ordereddistance method [21] applies to more general Gaussian mixtures, but under a fixed dimension $p$.

Our numerical results suggest that CESME is not overly sensitive to initialization, provided that the initializer is reasonable; see § 4.3 and Figure 6. Some initialization schemes considered there do not satisfy Condition 3.3, yet still perform satisfactorily in practice. This is compatible with our theory, since Condition 3.3 is sufficient rather than necessary. Nevertheless, from the practical perspective, it is desired to construct an initialization procedure that satisfies Condition 3.3 with high probability.

Motivated in part by spectral initialization methods for Plackett-Luce mixtures [34], we consider a spectral-based initialization in Algorithm 2. The key observation is that, in the latent space, $\mathbf{g}^{*}(\mathbf{X})$ follows a Gaussian mixture, so the Gram matrix $\mathbf{g}^{*}(\mathbf{X})^\top\mathbf{g}^{*}(\mathbf{X})$ admits a low-rank signal-plus-noise decomposition whose leading eigenspace reflects the cluster structure. Although the observed data $\mathbf{X}$ are obtained through unknown coordinate-wise monotone transformations, this structure is still approximately preserved when the transformations are sufficiently regular near the component centers. In turn, the leading eigenspace of $\mathbf{X}^\top\mathbf{X}$ continues to capture the dominant cluster geometry well enough to yield a suitable initializer, provided that the latent cluster signal is not masked by the effect of the unknown coordinate-wise transformations on the Gram matrix, together with the sampling noise. Since these transformations need not affect all spectral directions equally, we introduce a regularization parameter $\alpha_{1} > 0$ to down-weight weaker and less stable directions, thereby stabilizing the resulting spectral representation.

Proposition D.1 in Section D.2 of the supplementary file states the initialization guarantee for Algorithm 2. Under its stated assumptions, the initial estimators satisfy Condition 3.3 with high probability, up to a permutation of the component labels. The proof is provided by the argument for Proposition 3.1 of Zhang et al. [51], which uses random matrix theory [52, 1] and leave-one-out perturbation arguments [10] to control estimation of the leading eigenspace of $\mathbf{X}^\top\mathbf{X}$ and the subsequent $K$-means assignments and parameter estimates. Section D.2 also discusses Algorithm 2 A Spectral-Based Initialization

Input: $p$-dimensional data $\mathbf{X} = (\mathbf{X}_1,\ldots,\mathbf{X}_n)$, preconceived number of clusters $K$, and the normalization parameter $\alpha_1 > 0$.

1. Perform eigen-decomposition on $\mathbf{A} = \mathbf{X}^{\mathsf{T}}\mathbf{X}$ to obtain the rank-$K$ approximation $\mathbf{V}\mathbf{D}_0^2\mathbf{V}^{\mathsf{T}}$, where $\mathbf{D}_0^2$ is the diagonal matrix of top $K$ eigenvalues of $\mathbf{A}$.

2. Normalize $\widehat{\mathbf{D}} = (\mathbf{D}_0^2 + \alpha_1d_{0,\max}^2\mathbf{I})^{-1/2}\mathbf{D}_0$, where $d_{0,\max}^2$ is the largest eigenvalue.

3. Let $\widehat{\mathbf{z}}^{(0)}$ be the clustering assignments of $K$-means on $\widehat{\mathbf{D}}\mathbf{V}^{\mathsf{T}}$.

4. Let $\widehat{\lambda}_k^{(0)} = n^{-1}\sum_i \mathbb{I}(\widehat{z}_i^{(0)} = k)$ and $\widetilde{\boldsymbol{\mu}}_k^{(0)} = \{\sum_i \mathbb{I}(\widehat{z}_i^{(0)} = k)\}^{-1}\sum_i \mathbb{I}(\widehat{z}_i^{(0)} = k)\mathbf{X}_i$ for each $k$. Given $\widehat{\mathbf{z}}^{(0)}$ and $\widetilde{\boldsymbol{\mu}}_k^{(0)}$, let $(\widetilde{\sigma}_{\ell}^{(0)})^2 = n^{-1}\sum_{k=1}^{K}\sum_{i:\widehat{z}_i^{(0)}=k}(X_{i\ell}-\widetilde{\mu}_{k\ell}^{(0)})^2$ for each $\ell= 1,\ldots,p$.

5. Given $\widehat{\mathbf{z}}^{(0)}$, $\widehat{\boldsymbol{\lambda}}^{(0)}$, and $\{(\widetilde{\sigma}_{\ell}^{(0)})^2\}_{\ell=1}^{p}$, initialize $\widehat{g}^{(0)}$ by (2.2), $\widehat{\boldsymbol{\mu}}_k^{(0)} = \{\sum_i \mathbb{I}(\widehat{z}_i^{(0)} = k)\}^{-1}\sum_i \mathbb{I}(\widehat{z}_i^{(0)} = k)\widehat{g}^{(0)}(\mathbf{X}_i)$ for each $k$, and $\widehat{\boldsymbol{\Sigma}}^{(0)} = n^{-1}\sum_{k=1}^{K}\sum_{i:\widehat{z}_i^{(0)}=k}\{\widehat{g}^{(0)}(\mathbf{X}_i)-\widehat{\boldsymbol{\mu}}_k^{(0)}\}\{\widehat{g}^{(0)}(\mathbf{X}_i)-\widehat{\boldsymbol{\mu}}_k^{(0)}\}^{\mathsf{T}}$.

Output: Initial parameters $\widehat{\boldsymbol{\beta}}^{(0)} = (\widehat{\boldsymbol{\lambda}}^{(0)},\widehat{\boldsymbol{\mu}}^{(0)},\widehat{\boldsymbol{\Sigma}}^{(0)})$ and $\widehat{g}^{(0)}$.

how the initialization assumptions relate to the conditions for the subsequent estimation and clustering theory.

Simulations

Clustering accuracy

In this section, we compare CESME with a comprehensive set of methods, including CHIME, sparse $K$-means (SK-means), IF-PCA, spectral clustering (SC), sparse spectral clustering (SSC; Lu et al., (2016); Wang et al. [46]). Due to space limitations in the figures, we present the results of other competitors in Section E.1 in the supplementary file, including $K$-means, model-based clustering (Mclust; Fraley and Raftery, [11]), multivariate mixture model-based clustering (Nonpara; Benaglia et al., [4]), sparse convex clustering (SCV; Wang et al., [44]), Seurat (Hao et al., [18]), Sincera (Guo et al., [15]), significance analysis (SA; Grabski et al., [14]), and nonparametric mixture model (NPEM; Hall and Zhou, [16]).

For CESME, we use the penalty selection procedure described in § 2.3, with initial penalty $\rho^{(0)} = \sqrt{\log(np)/n}$, numerical decay parameter $\nu= 0.8$, and candidate grid $\mathcal{C} = \{2^j, j = -4,\ldots,4\}$. We also include CESME-oracle, which assumes the true transformations $g^*$ are known. As a baseline, we report the Bayes rate, which is the accuracy of the optimal Bayes rule under the true model.

To generate data $\mathbf{X}_i \in\mathbb{R}^p$ for $i = 1, 2, \ldots, n = 400$, we first draw latent variables $\mathbf{Z}_i$ independently from a Gaussian mixture $\sum_{k=1}^{K} \lambda_k \mathcal{N}(\boldsymbol{\mu}_k, \mathbf{\Sigma})$, with parameters defined below. For each $\ell$, the samples $Z_{1\ell}, \ldots, Z_{n\ell}$ are standardized to have zero mean and unit variance. The observed data are then obtained via component-wise inverse transformations: $X_{i\ell} = g_{\ell}^{-1}(Z_{i\ell})$ for each $i$, where the same transformation $g^{-1}$ is applied across all dimensions. We consider four transformations: (1) logarithm: $g^{-1}(x) = \exp(x)$; (2)

the rescaled inverse $t$-CDF: $g^{-1}(x) = F^{-1}((x - c)/d)$, where $F$ is the CDF of the $t_{3}$ distribution and $c$, $d$ rescale $x$ to $(0, 1)$; (3) polynomial: $g^{-1}(x) = x^{3}$; and (4) identity: $g^{-1}(x) = x$. We set $K = 4$ clusters with unequal proportions: $\lambda_{1} = 0.2$, $\lambda_{2} = 0.23$, $\lambda_{3} = 0.27$, and $\lambda_{4} = 0.3$. Two covariance structures are considered: the inverse AR (1) model, $\boldsymbol{\Sigma} = \boldsymbol{\Omega}^{-1}$ with $\Omega_{\ell\ell'} = 0.8^{|\ell-\ell'|}$; and the asymptotic quadratic decay (AQD) model with $\sigma_{\ell\ell'} = 0.7/(0.7 + |\ell-\ell'|)^{2}$. For each $\boldsymbol{\Sigma}$, the cluster means are defined as $\boldsymbol{\mu}_{1} = \boldsymbol{\Sigma}(\boldsymbol{\gamma}_{1} + \boldsymbol{\gamma}_{2} + \boldsymbol{\gamma}_{3})$, $\boldsymbol{\mu}_{2} = -\boldsymbol{\Sigma}\boldsymbol{\gamma}_{1}$, $\boldsymbol{\mu}_{3} = -\boldsymbol{\Sigma}\boldsymbol{\gamma}_{2}$, and $\boldsymbol{\mu}_{4} = -\boldsymbol{\Sigma}\boldsymbol{\gamma}_{3}$, where $\boldsymbol{\gamma}_{1}$, $\boldsymbol{\gamma}_{2}$, and $\boldsymbol{\gamma}_{3}$ are sparse signal vectors: $\boldsymbol{\gamma}_{1} = 4\delta(1,\ldots,1,0,\ldots,0)^{\mathsf{T}}$ with nonzero entries in positions 110; $\gamma_2 = 4\delta(0,\ldots,0,1,\ldots,1,0,\ldots,0)^\top$ in positions 11 to 20; and $\gamma_3$ similarly in positions 21 to 30. The parameter $\delta$ controls signal strength and thus governs cluster separability.

The theoretical conditions for initialization with Algorithm 2 are satisfied under both covariance settings, across all signal strengths and dimensions. For simulations, the true value of $K$ is provided to all competing methods.

Figure 3 and Figure 4 display clustering accuracies for $\delta= 2.5$ with $p$ ranging from 100 to 1000, and for varying $\delta\in[0.25, 2.5]$ with fixed $p = 800$, respectively. CESME consistently outperforms all competing methods under the three non-identity transformations.

Vector (PDF)

Although CHIME and SSC generally perform better than the other methods, they still fall notably short of CESME. The gap between CESME and CESME-oracle is negligible, suggesting that the transformation estimation (2.2) is satisfactory and robust, even with complex covariances and unbalanced clusters. When the transformation is identity, CESME performs nearly identically to CHIME, which in this case aligns with CESME-oracle. The numerical results support the theory in § 3. The comparison of CESME with additional competitors is presented in Figures S.3 and S.4 in Section E of the supplementary file. It is observed that $K$-means, Mclust, Nonpara, and NPEM perform worse than the competitors displayed in Figure 3 and Figure 4 even under identity transformations, since they do not exploit model sparsity. CESME continues to dominate these competitors.

On the selection of $K$ via the data-driven BIC-type criterion

To assess CESME’s ability to select $K$ using the data-driven BIC-type criterion $\mathrm{Cr}_{K}$, we used the same data generation models with an additional setting involving heterogeneous transformations across dimensions (denoted “combination”; see Section E of the supplementary file). For comparison, we included the gap statistics [41] which is a very popular technique for determining the number of clusters. Simulations were conducted with $n = 400$, $K = 4$, $p = 100, 200, \ldots, 1000$ and $\delta= 0.25, 0.5, \ldots, 2.5$.

For each setting, CESME was run with $K \in\{2, \ldots, 8\}$, and the value maximizing $\mathrm{Cr}_{K}$ was selected. Results based on 100 replicates are reported in Table 1. Overall, the proposed data-driven BIC-type criterion ($\mathrm{Cr}_{K}$) achieves strong accuracy in selecting the correct $K$, substantially outperforming the gap statistic.

	Covariance	logarithm	$t$-CDF	polynomial	identical	combination
Cr	inverse AR (1)	0.88	0.88	0.87	0.90	0.82
	AQD	0.85	0.88	0.84	0.92	0.79
Gap	inverse AR (1)	0.58	0.61	0.59	0.67	0.55
	AQD	0.50	0.55	0.52	0.59	0.51

We further compared our proposed data-driven BIC-type criterion with several additional methods for selecting $K$, including prediction strength (PS; [40]), gap statistic (GAP; [41]), rank detection (RD; [23]), and the eigenvalue-ratio test (ER). Figure 5 summarizes the results. Overall, the proposed data-driven BIC-type criterion outperforms the competing methods.

Vector (PDF)

Impact of initialization

Furthermore, we conduct numerical experiments to examine the performance of CESME under different initialization strategies with the same simulation settings in Section 4. For each setting, we run CESME initialized by Algorithm 2, $K$-means, sparse $K$-means (SK-means), spectral clustering (SC), and spectral clustering with the $(1+\epsilon)$-approximate $K$-means algorithm (SC-approx). Since the initialization methods except Algorithm 2 require the true number of clusters $K$ to provide valid initialization, we employ the true $K$ for this simulation. As shown in Figure 6, CESME is insensitive to initialization as long as the initialization is reasonable. The proposed initialization, namely Algorithm 2 is consistently slightly better than other initialization schemes.

Real Data Analysis

To demonstrate the performance of CESME, we apply it to six high-dimensional scRNA-seq benchmark data sets from Buettner et al. [5]; Deng et al. [9]; Schlitzer et al. [38]; Pollen et al. [37]; Ting et al. [42]; Treutlein et al. [43], which have also been used in prior studies benchmarking clustering methods in scRNA-seq analysis (Wang et al., [46]; Park and Zhao, [36]). Each dataset is preprocessed for normalization and includes ground truth cell annotations. Although clustering is performed in an unsupervised manner, the true labels are used post hoc to evaluate clustering accuracy.

For a comprehensive comparison, we include several methods that performed competitively in the simulation studies: CHIME, $K$-means, sparse $K$-means (SK-means), IF-PCA, spectral clustering (SC), and sparse spectral clustering (SSC). We exclude multivariate mixture model-based methods and sparse convex clustering due to their instability in high-dimensional settings. We also compare CESME with three widely used pipelines, as noted in the Introduction, for a more scientifically relevant evaluation: significance analysis (SA), Sincera, and Seurat. To determine the number of clusters $K$, we use the data-driven BIC-type criterion proposed in § 2.3, evaluated over $K = 2, \ldots, 12$, for both CESME and CHIME. For SSC, we apply the eigenvalue-ratio test, selecting $K = \arg\max_{k>1} \sigma_k(\mathbf{X}^{\mathsf{T}}\mathbf{X})/\sigma_{k+1}(\mathbf{X}^{\mathsf{T}}\mathbf{X})$, where $\sigma_k(\mathbf{X}^{\mathsf{T}}\mathbf{X})$ denotes the $k$th largest eigenvalue of $\mathbf{X}^{\mathsf{T}}\mathbf{X}$. For $K$-means and SK-means, we adopt the prediction strength test and set $K$ as the largest number of clusters for which the prediction strength exceeds 0.8 (Tibshirani and Walther, [40]). For IF-PCA, we follow its originally proposed rank detection method. For SA, Sincera, and Seurat, we use their default pipelines to select $K$.

As shown in Table 2, CESME consistently achieves the highest clustering accuracy on five of the six datasets and most accurately estimates $K$, including the challenging Pollen dataset with 11 true cell clusters, one of which contains only seven cell samples. Identifying $K$ in such cases is difficult, as small clusters tend to be absorbed into larger ones, leading to underestimation. On the Deng dataset, CESME performs comparably to the widely used benchmark, Seurat. The estimated sparsity levels by CESME are 325, 163, 62, 117, 53, and 230, across the datasets, suggesting that scRNA-seq data exhibit underlying sparse structures. It is worth noting that for datasets with a large number of true clusters, random clustering would result in very low accuracy, and methods that fail to accurately estimate $K$ also tend to perform poorly. Furthermore, substantial non-Gaussianity, as evidenced by gene-wise kurtosis and skewness in Figure S.12 in the supplementary file, makes accurate selection of $K$ even more difficult for other methods and highlights the advantage of CESME in these settings.

Dataset	$n, K$	$p$	CESME	Seurat	IF-PCA	SA	$K$-means	SK-means	CHIME	SC	Sincera	SSC
Schlitzer	251, 3	11835	83.9 (3)	53.8 (4)	59.0 (2)	63.3 (2)	64.5 (2)	66.1 (2)	47.1 (5)	57.8 (4)	30.7 (5)	41.4 (4)
Pollen	249, 11	14806	79.1 (9)	77.9 (7)	66.7 (6)	70.8 (9)	27.3 (2)	27.3 (2)	57.4 (2)	27.3 (2)	52.4 (6)	18.5 (2)
Deng	135, 7	12549	73.3 (7)	73.3 (4)	52.6 (3)	46.7 (2)	45.2 (2)	45.2 (2)	50.2 (5)	50.4 (9)	28.9 (3)	29.6 (9)
Ting	114, 5	14406	79.1 (5)	78.9 (3)	64.9 (3)	68.9 (5)	62.3 (2)	65.8 (2)	26.7 (2)	57.9 (2)	32.5 (4)	22.8 (2)
Treutlein	80, 5	959	68.5 (5)	66.2 (2)	60.0 (2)	60.3 (2)	66.2 (2)	65.0 (2)	61.6 (2)	45.0 (11)	36.2 (3)	30.0 (11)
Buettner	182, 3	8989	68.1 (3)	67.6 (4)	64.8 (2)	44.5 (2)	64.8 (2)	60.4 (2)	55.8 (2)	50.0 (4)	25.3 (6)	40.7 (4)

Conclusions and Discussions

In this paper, we introduce a novel semiparametric mixture model that forms the basis of CESME, a new high-dimensional clustering method, which enjoys both strong empirical performance, outperforming widely used pipelines for scRNA-seq analysis, and rigorous theoretical guarantees. Statistically, effective clustering naturally aligns with a mixture model, making model-based clustering a principled approach to the problem. The success of CESME reinforces this perspective. While CESME currently employs an EM-type alternating maximization algorithm, exploring more efficient optimization strategies remains an important direction for future work. The computational issue becomes even more critical when we extend the method to high-dimensional tensor data clustering. Recently, Mai et al. [30] developed DEEM based on a mixture of tensor normals for high-dimensional clustering. It is natural to consider the tensor generalization of CESME, which will be reported in a separate paper.

The online supplementary file provides the following materials: (A) proofs of Theorem 3.1 with necessary technical lemmas, (B) discussions of identifiability of the model, (C) proofs of technical lemmas to establish Theorem 3.1, (D) more technical lemmas and discussions on the initializations, (E) additional numerical results for the simulation study and real data analysis, (F) simulations of one-step EM variant of CESME and performance of CESME under mis-specified models.

Data Availability

The datasets analyzed in this study are publicly available from previously published studies and are described in § 5. R code for reproducing the numerical studies and data examples can be found at https://github.com/lyuouzhang/CESME.