Regularized k-POD Clustering

• Regularized k-POD clustering is a k-means variant that handles incomplete data by integrating imputation with feature-wise regularization.
• It utilizes a majorization–minimization framework to iteratively update cluster assignments and shrink noise features, reducing bias in high dimensions.
• Empirical studies, including simulations and single-cell RNA sequencing applications, demonstrate its superior accuracy over standard methods in handling missing data.

Regularized $k$-POD clustering is a sparse, regularized extension of $k$-means for incomplete data matrices. It starts from $k$-POD, which replaces the complete-data $k$-means objective by the squared Frobenius norm restricted to observed entries, and adds a feature-wise penalty on the cluster-center matrix to reduce bias in high-dimensional missing data with many irrelevant features. The method is motivated by two complementary observations: the original $k$-POD algorithm is computationally efficient and often effective for partially observed data, but its population target differs from full-data $k$-means under missingness, and in high dimensions the presence of noise features can substantially bias the estimated centers (Chi et al., 2014, Terada et al., 2024, Guan et al., 16 Jul 2025).

1. Observed-entry $k$-means and the origin of $k$-POD

Classical $k$-means assumes a fully observed matrix $X\in\mathbb{R}^{n\times p}$ and solves

$k$0

where $k$1 is the cluster-membership matrix with $k$2, and $k$3 contains the $k$4 cluster centers as rows. When entries are missing, this objective is not directly evaluable because distances and cluster means are undefined unless incomplete cases are removed or missing entries are imputed.

$k$5-POD addresses this by introducing the projection operator $k$6 that keeps observed entries and sets missing ones to zero. Its core objective is

$k$7

Equivalently, with missingness indicators $k$8, the empirical loss can be written as

$k$9

If there is no missingness, $k$0 is the full index set and $k$1-POD reduces exactly to ordinary $k$2-means (Chi et al., 2014).

The original $k$3-POD algorithm uses a majorization–minimization scheme. Missing entries are filled with current cluster-center predictions, producing a completed matrix, and ordinary $k$4-means is then run on that matrix. This outer MM loop and inner Lloyd-type clustering step give monotone descent of the observed-entry loss and keep the computational profile close to standard $k$5-means. A central practical attraction of this formulation is that it does not require a specified missingness mechanism, external information, or a separate imputation model (Chi et al., 2014).

2. Why unregularized $k$6-POD requires correction

The principal theoretical criticism of vanilla $k$7-POD is that, under Missing Completely At Random (MCAR), its population objective is not the full-data $k$8-means objective. If $k$9 denotes a generic missingness pattern and $k$0 is the expected $k$1-means loss computed only on coordinates with $k$2, then

$k$3

Thus the expected $k$4-POD loss is a weighted sum of $k$5-means losses over projected subspaces. The estimator is consistent for the minimizer of $k$6, but not in general for the minimizer of the full-dimensional $k$7-means loss $k$8. In the paper’s formulation, $k$9-POD is therefore “consistent for the wrong target” (Terada et al., 2024).

This decomposition yields a precise geometric interpretation. $k$0-POD implicitly clusters many coordinate projections at once, with weights determined by missingness-pattern probabilities. Variables with higher observation probabilities appear in more patterns and therefore exert greater influence on the objective. The resulting effective feature weighting can distort cluster boundaries relative to full-data $k$1-means. In a two-dimensional Gaussian mixture example with $k$2, the $k$3-POD boundary is described as substantially different from the full-data boundary (Terada et al., 2024).

The 2025 regularized formulation adds a second, high-dimensional motivation. In settings with many irrelevant features, $k$4-POD can produce biased cluster-center estimates even under MCAR because its weighted mixture of subspace objectives includes many subsets dominated by noise features. Empirically, when $k$5 and both features are relevant, $k$6-POD is nearly unbiased and matches full-data $k$7-means; when $k$8 with 2 relevant and 98 noise features, the estimated centers are severely biased in every feature (Guan et al., 16 Jul 2025).

The critique is not absolute. Theoretical inconsistency does not preclude practical utility. For high-dimensional data with low missing rate per variable, $k$9-POD can still be preferable to complete-case analysis. One simulation in the 2024 paper reports $k$0, $k$1, missing rate $k$2 in each variable, and $k$3, with MSE of centers relative to oracle $k$4-means equal to $k$5 for the oracle, $k$6 for complete-case $k$7-means, and $k$8 for $k$9-POD (Terada et al., 2024). Regularization is therefore introduced not because $k$0-POD is uniformly unusable, but because its bias becomes structurally important in high-dimensional, noisy-feature regimes.

3. Formal regularized $k$1-POD objective

Regularized $k$2-POD modifies the observed-entry loss by penalizing the cluster-center matrix columnwise: $k$3 with $k$4, $k$5, and $k$6. The first term is exactly the $k$7-POD loss; the second is a feature-wise regularizer designed to shrink or eliminate entire features across all clusters. The method reduces to $k$8-POD when $k$9, and to regularized $k$0-means when there is no missing data (Guan et al., 16 Jul 2025).

A feature $k$1 is called a noise feature if the true optimal $k$2-means centers satisfy $k$3. The regularization acts directly on the column vector $k$4, so entire features are selected or shrunk jointly rather than cluster by cluster.

Penalty	Definition	Effect
$k$5	$k$6	Hard feature selection
$k$7	$k$8	Group-lasso shrinkage

The $k$9-type penalty $k$0 counts the number of nonzero columns of $k$1, while the group-lasso penalty $k$2 uses Euclidean norms in $k$3 and feature-specific weights $k$4. In the proposed implementation, these weights are adaptive: one first runs $k$5-POD to obtain $k$6, then sets

$k$7

Features with small $k$8-POD center norms receive stronger shrinkage, while features with larger norms are penalized less heavily (Guan et al., 16 Jul 2025).

This formulation can also be understood as a structured matrix-factorization problem with missing entries. The loss $k$9 resembles low-rank matrix completion, but the target is not full matrix recovery. The goal is clustering with sparse or shrunken centers.

4. Majorization–minimization and Lloyd-type optimization

The regularized algorithm preserves the MM structure of $X\in\mathbb{R}^{n\times p}$0-POD. Given current iterates $X\in\mathbb{R}^{n\times p}$1, it defines the majorizing function

$X\in\mathbb{R}^{n\times p}$2

With the imputed matrix

$X\in\mathbb{R}^{n\times p}$3

this becomes

$X\in\mathbb{R}^{n\times p}$4

Each outer iteration therefore reduces to solving a regularized $X\in\mathbb{R}^{n\times p}$5-means problem on a complete matrix whose missing entries have been filled by current cluster means. The resulting objective decreases monotonically, and the algorithm converges to a local minimum (Guan et al., 16 Jul 2025).

The inner optimization is Lloyd-style alternating minimization. Given current centers, the assignment update is the usual nearest-center rule on $X\in\mathbb{R}^{n\times p}$6: $X\in\mathbb{R}^{n\times p}$7 Given updated memberships $X\in\mathbb{R}^{n\times p}$8, the center update is columnwise.

For $X\in\mathbb{R}^{n\times p}$9, letting $k$00 and $k$01, the update is a hard-thresholded cluster mean: $k$02 For $k$03, a single inner MM step for the group-lasso penalty yields

$k$04

This is a ridge-regularized mean whose shrinkage strength depends on $k$05, $k$06, and the current column norm (Guan et al., 16 Jul 2025).

The algorithm inherits the computational economy of $k$07-POD. Per outer iteration, the imputation cost is $k$08, where $k$09 is the observed fraction, and the regularized $k$10-means step is dominated by $k$11 when $k$12 is large. Asymptotically, each outer iteration is about $k$13, where $k$14 is the number of inner assignment/center alternations. The reported runtime is similar to $k$15-POD and far cheaper than multiple imputation procedures (Guan et al., 16 Jul 2025, Chi et al., 2014).

5. Sparsity, shrinkage, and bias reduction

The central structural result for regularized $k$16-POD is a featurewise thresholding characterization of the estimated centers. Under MCAR, let $k$17 be the observed proportion of feature $k$18, let $k$19 be its marginal variance based on observed entries, let $k$20 be the within-cluster sum of squares in feature $k$21, and let $k$22 denote the minimum one-dimensional $k$23-means loss in that feature. Then the paper states the following (Guan et al., 16 Jul 2025).

For $k$24, if

$k$25

then $k$26. Otherwise $k$27 and the estimated centers equal the sample cluster means in that feature. For $k$28, if

$k$29

then $k$30; otherwise the estimate is nonzero and each center is a shrunken version of the sample cluster mean: $k$31

These threshold conditions make the regularizer’s role explicit. The quantity

$k$32

measures how much clustering in feature $k$33 reduces variance. Small values indicate that the feature contributes little to cluster separation, in which case $k$34 discards it and $k$35 either discards or strongly shrinks it. Large values preserve the feature, with either no shrinkage ($k$36) or controlled shrinkage ($k$37).

This directly addresses the bias mechanism of $k$38-POD. Since vanilla $k$39-POD optimizes a weighted sum of subspace $k$40-means losses, subsets composed mostly of noise features can pull the minimizer away from the true centers. Regularization restricts the admissible center matrix $k$41 by forcing many noise-feature columns toward zero. The paper’s interpretation is that many of the biased component losses are thereby “de-biased,” and the aggregate minimizer is pulled closer to the truth (Terada et al., 2024, Guan et al., 16 Jul 2025).

The method does not eliminate all limitations. The paper notes that bias reduction depends on effective tuning and weight selection, and that formal feature-selection guarantees are not fully established. In the absence of missingness, however, the stated sparsity and shrinkage results coincide with those for regularized $k$42-means (Guan et al., 16 Jul 2025).

6. Simulation evidence and single-cell RNA sequencing applications

The simulation study in the 2025 paper uses a $k$43-component Gaussian mixture with diagonal covariance and equal mixing weights $k$44. The main settings are $k$45, $k$46, with either $k$47 or $k$48. The first $k$49 features are relevant and the remaining $k$50 are noise. Missingness mechanisms include MCAR, MAR, MNAR1, and MNAR2. Evaluation uses MSE of centers, CER, and predictive CER, with competitors including complete-case $k$51-means, multiple imputation with mice followed by $k$52-means, $k$53-POD, and regularized $k$54-POD with both $k$55 and $k$56 (Guan et al., 16 Jul 2025).

The main empirical findings are strongly asymmetric between moderate and high dimension. In almost all settings, regularized $k$57-POD dramatically reduces MSE compared with $k$58-POD, complete-case clustering, and mice. For $k$59, the $k$60 version often yields the smallest MSE. For $k$61, the group-lasso version is described as more stable and typically best: under MCAR with $k$62–$k$63 missingness, its MSE is around $k$64–$k$65, while $k$66-POD has MSE $k$67–$k$68. CER and predictive CER are also systematically lower in most cases; under $k$69, MCAR $k$70, predictive CER drops from about $k$71 for $k$72-POD to about $k$73 for the group-lasso regularized version. One reported exception is that, in some MCAR $k$74 settings with very large missingness and strong relationships among relevant features, mice can outperform regularized $k$75-POD in MSE because its assumptions are satisfied and it can impute accurately (Guan et al., 16 Jul 2025).

Two real-data applications are single-cell RNA-seq datasets with heavy dropout. In the Usoskin dataset, with $k$76 cells, $k$77 genes, 4 known sensory neuron subtypes, and about $k$78 zeros, the reported CERs are $k$79 for zero-imputation plus $k$80-means, $k$81 for scImpute plus $k$82-means, $k$83 for $k$84-POD, $k$85 for regularized $k$86-POD with $k$87, and $k$88 for regularized $k$89-POD with group lasso. In the Treutlein dataset, with $k$90, $k$91, $k$92 states, and about $k$93 missingness, the CERs are $k$94, $k$95, $k$96, $k$97, and $k$98, respectively. In both datasets, the group-lasso variant achieves the lowest error (Guan et al., 16 Jul 2025).

The application-level interpretation in the paper is that feature-wise regularization can select informative genes and mitigate the effect of thousands of noisy or dropout-affected coordinates. By contrast, vanilla $k$99-POD is reported to misclassify many cells in these settings, especially when dropout and high dimensionality interact with irrelevant genes.

7. Tuning, initialization, related methods, and limitations

The regularization strength $k$00 is selected from candidate values typically on a log scale, such as $k$01. Two criteria are described. The preferred one is instability: split the data into training subsets and a validation subset, fit the clustering separately on the training subsets for a fixed $k$02, assign validation points under both fitted clusterings, and measure disagreement probability. The alternative is

$k$03

where $k$04 is the number of active features. Instability is reported to be more robust, especially for the group-lasso penalty and under MNAR settings, whereas BIC is less reliable when its missingness assumptions are violated (Guan et al., 16 Jul 2025).

Initialization materially affects performance. Two strategies are used: comp, which runs $k$05-means++ on complete cases and then assigns all observations using observed features only, and impt, which pre-imputes each feature by its observed mean, samples $k$06 rows from the pre-imputed matrix as initial centers, and then assigns points to those centers. The paper reports that impt is generally preferable, especially under high missingness or large $k$07. Multiple random starts are recommended, and for very high dimensions a sparse initialization based on $k$08-POD centers with thresholded columns is suggested (Guan et al., 16 Jul 2025).

Within the broader clustering literature, regularized $k$09-POD occupies a distinct position. Sparse $k$10-means is designed for complete data and typically regularizes feature weights rather than cluster centers. EM for Gaussian mixtures with missing data can integrate over missing entries but imposes a parametric model and can be computationally heavy in high dimension. Matrix completion followed by clustering targets whole-matrix reconstruction rather than center estimation. Regularized $k$11-POD instead integrates imputation and clustering in a single objective built around sparse centers (Guan et al., 16 Jul 2025).

A broader regularization paradigm for center-based clustering introduces a garbage cluster and solves

$k$12

so that points may be discarded at fixed cost $k$13. This formulation is not a missing-data method, but it provides an adjacent view of regularization as an explicit trade-off between center fit and non-assignment (Kushagra et al., 2017). This suggests a broader design space for $k$14-POD-type methods in which regularization can target not only noise features, as in the 2025 method, but also unreliable or outlying observations.

The main limitations reported for regularized $k$15-POD are those of $k$16-means-type procedures more generally: sensitivity to highly overlapping or non-spherical clusters and to cluster imbalance. Performance can deteriorate when the missing proportion is very large or the missingness mechanism is very complex, especially under MNAR, because the method relies on cluster-wise mean imputation that is described as unbiased primarily under MCAR/MAR-like mechanisms. The $k$17 version is more sensitive to initialization and may require more random starts, whereas the group-lasso version is reported as more stable in high dimensions (Guan et al., 16 Jul 2025).