Adaptively Robust and Sparse $K$-means Clustering
Authors: HAO LI
Shonosuke Sugasawa
Shota Katayama
Venue: Accepted by TMLR ()

Abstract: While $K$-means is known to be a standard clustering algorithm, its performance may be compromised due to the presence of outliers and high-dimensional noisy variables. This paper proposes adaptively robust and sparse $K$-means clustering (ARSK) to address these practical limitations of the standard $K$-means algorithm. For robustness, we introduce a redundant error component for each observation, and this additional parameter is penalized using a group sparse penalty. To accommodate the impact of high-dimensional noisy variables, the objective function is modified by incorporating weights and implementing a penalty to control the sparsity of the weight vector. The tuning parameters to control the robustness and sparsity are selected by $\rm Gap$ statistics.Through simulation experiments and real data analysis, we demonstrate the proposed method's superiority to existing algorithms in identifying clusters without outliers and informative variables simultaneously.

Published in Transactions on Machine Learning Research (11/2024)
Adaptively Robust and Sparse K-means Clustering
Hao Li
lee 19c1995hao 0keiioooo@@ko.pp Skeio.jjp
Graduate School of Economics Keão University
Shonosuke Sugasawa
sugosawofBcon..eio.cco
Faculty of Economics Keio University
Shota Katayama
katayama@ccon.eiio.com
Paculty of Economics
Keio University
Reviewed on OpenReview: https: httpe://openrevieel openreview.
net/ forum? id= f
Abstract
While K-means is known to be a standard clustering algorithm, its performance may be compromised due to the presence of outliers and high-dimensional noisy variables. This paper proposes adaptively robust and sparse K-means clustering (ARSK) to address these practical limitations of the standard K-means algorithm. For robustness, we introduce redundant error component for cach observation, and this additional parameter is penal- ized using group sparse penalty. To accommodate the impact of high-dinensional noisy variables, the objective function is modified by incorporating weights and implementing penalty to control the sparsity of the weight vector. The tuning parameters to control the robustness and sparsity are selected by Gap statistics. Through simulation experiments and real data analysis, we demonstrate the proposed hod's superiority existing algorithms in identifying clusters without outliers and informative variables simultaneously.
1 Introduction
Identifying clustering structures is recognized as an important task in pursuing potential heterogeneity in datasets. While thore are a variety of clustering mothods, K-means clustering (Forgy 1965 is the most standard clustering algorithm and is widely used in many scientific applications. However, real-world datasets present significant difficulties, such as the presence of outliers and high-dimensional noisy variables. Regarding the outlier issues, there some attempts to robustify the standard K-means such trimmed K-means (Cuesta-Alberte Albertos et al. 1997) and robust K-means (Klochkov et al. 2021). On the other hand, high-dimensional noisy characteristics are typically addressed through clustering approaches that focus on sparsity for variable selection. In particular, Wiitten & Tibshirani (2010) proposed a lasso-type penalty to the variable weights incorporated in the clustering objective function. However, the efficacy of this algorithm might be lost if the dataset contains a significant. number of outliers. While several methods address either of the two aspects (the existence of outlicrs and noisy variables), practically useful methods to address both aspects simultaneously are scarce. One exception is the method proposed by Kondo al. (2016). However, this approach requires the assumption of trimming level in advance. In contrast, Brodinová et al. (2019) introduces approach beginning with climinating outliers identified by the model. However, this methodology potentially leads to reduced sample size and may result in information loss. Moreover, the excessive variations in gradient magnitudes may destabilize when optimizing the variable weight process unless all outliers are eliminated in advance.
In addition clustering methods based on objective functions, probabilistic clustering using mixture models is also popular. It is typically done by fitting multivariate Gaussian mixtures. While several approaches

---

Published in Transactions on Machine Learning Research (11/2024)
have been proposed for stable estimation and clustering under existence of outliers (e.g. Coretto & Hennig
2016 Fujisawa & Eguchi 2006 Punzo & McNicholas 2016; Sugasawa & Kobayashi 2022 Yang et al. 2012
such approach suffers from a huge number of parameters (especially for modelling covariance matrix) to be estimated when the dimension is large. Therefore, these methods would not be a reasonable solution when our primary focus is on clustering.
To address the difficulty of a series of sparse trimming- and probabilistic-lluster approaches, we propose approach based on regularization techniques for sparsity well as robustness, called adaptively robust and sparse K-means (ARSK). We introduce parameter to absorb the undesirable effects caused by outliers, making the clustering algorithm robust. We employ group lasso penalty and a group smoothly clipped absolute deviation (SCAD) penalty (Wang et al. 2007) to estimate these parameters during clustering steps. Meanwhile, to reduce the data dimension, we also introduce a SCAD penalty (Fan & Li 2001) to weights for each variable to eliminate noise variables irrelevant to clustering. Based on our analysis of real-world datasets, in the ultra-high dimensional setting (i.e., the dataset. dimension exceeding 700), the result. procured through the implementation of the lasso penalty proposed by Witten & Tibshirani (2010) frequently poor performance. Therefore, have developed methodology based on the SCAD penalty. The strong sparsity property of the SCAD penalty effectively improves the selection of relevant variables in the ultra-high dimensional scenario. We then develop an efficient computation algorithm to minimize the proposed objective function. While the proposed method has two tuning parameters controlling robustness (the number of outliers to be detected) and sparsity (the level of the sparsity of the variable weight), we adopt modified version of the Gap statistics (Tibshirani et al. 2007) to determine the optimal values of them. Hence, the level of robustness and sparsity is adaptively determined in fully data-dependent manner.
component for each observation
additional
The paper is organized as follows. Section 2 reviews the K-means algorithm and its extended model intro- duces the proposed method ARSK and presents its details and pseudocode. We also present the selection approach of tuning parameters and its pscudocode. while Section 3 presents the tuning paramcter selection result. Both two algorithms are applied to various scenarios of artificial datasets and several real-world datasets. This paper ends in Section with concluding remarks. R code implementing the proposed method is available at the GitHub repository hettts////gttaab..cm//ee99999..
2 Robust and Sparse K-means Clustering
2.1 Conventional K-means algorithm
Suppose that we observe a p-dimensional vector, X;; = (X,).... for = I,..., and that we are interested in clustering X Given the mumber of clusters, K, the conventional K-means algorithm is to find the optimal clustering by minimizing the following objective function:
min
|
X2,111
(1)
C1 CL
Thk
i
where |X,,- - 1 is the squared value of Lisdttaance between two observed vectors. Here e is a set of indices representing a cluster of observations and nę the size of ce, such that UK_ = (1....,n) n) and
Σ n4 n, where c Scc = o for l * k., The objective function (1] can be easily optimized by iteratively updating cluster means and assignment (e.g. Lloyd 1982).
As demonstrated in Witten & Tibshirani (2010), the objective function can be reformulated in terms of the detween-cluster sum of squares, given by max.....c. Σj= (X θ), where
n (x ΣΣ (X xx,))²
Σ (Xy-X19)) Xx,)²,
(2)
k= nk
2

---

Katayama & Published in Transactions on Machine Learning Research (11/2024)
with X., = (X.j....., Xny) and θ being collection of C1.... CK. Witten & Tibshirani (2010) employed the modified version of the objective function (2), given by Σド= 4, Wj Qx (X.,; 6) with an appropriate constrained on wy; for conducting variable selection and clustering simultaneously.
2.2 Robust and sparse K-means algorithm
The objective function ? is sensitive to outliers since the approach relies on the L2 distance between each observation. As result, even small quantity of deviating observations can affect the K-means algorithm. In order to robustify the K-means approach, we modify the function (x as follows:
(
Σ
Σ
(Xx) - E.j) - -Xvs - E3
(3)
E 1(1 - -(Xus - Ex
i.PEc
where E is an n X p error matrix, which is a collection of Ei and Eg is an additional location parameter for each Xi such that Faj will be non-zero values for outlying observations Xij to prevent undesirable effects from outliers. While the number of elements of E is the same as the number of observations, it can be assumed that E is sparse in the sense that most elements of E are 0. That is, most elements are not outliers. To stably estimate E, we consider penalized estimation for Eg when optimizing the function (3) Such approaches are used in the robust fitting of regression models (e.g. She & Owen 2011 Fujisawa, 2017).
Based on the robust objective function 3 we propose an objective function for robust and sparse K- K-means algorithm as follows:
1 Σ E.. Σ (Pa(w) X2) ww]).
L(O, E. A; Σ Ox O E.g)
where P1 and P2 are group penalty functions and penalty functions for sparse (shrinkage) estimation of E,; and Wy,, respectively, and = (A1, set luning parameters. To find unique solution for the weight coefficients Wj, we add a quadratic term 1 ドー= w/ the objective function. Then, the robust and sparse K-means clustering can be defined as Ô such that
(θ, E, () = argmax E,WE L(O, E, w. A),
H = wi RR:
Σ 2=1
In what follows, we employ two specific forms of group penalty functions Pi [[[] E2, A1): one type is the convex group penalty function, which is the group lasso penalty as other type is the nonconvex group penalty function, which is the group SCAD penalty (Huang et al. 2012)
Yuan & Lin 2006), and the
as
A E,slle
if Er,: 1l2 I A1; if \ E S adı; if E. ll > ad1.
PSCAD (1))),1))=
(1) E,, 2adi E., + 12)/2(a 1)
=
(a + /2
We also consider both convex and nonconvex penalties for the penalty term Pa A2): the convex standard lasso penalty P2 (a A2) = da Se|wal and, the non-convex SCAD penalty expressed as Pe(Wy; A2) = PSCAD (w;),
where PSCAD (wy) define as:
d2 |wm1
if |wy] S A2;
PSCAD (wy) = - -(wy 2ad2 wyy + x3)/2(a N2
if X < [ S aAz;
(a+ + 1]23/2
if |wy] > a\2-
3

---

Published in Transactions on Machine Learning Research (11/2024)
When applying the L1 or SCAD penalties, seek the optimal weight vector of w (Wi,.... Wp), which is normalized to have unit Lo-norm, ie., VEi= wf 1. This approach helps to promote more balanced solution Zou & Hastie 2010
2.3 Optimization algorithm
The block-coordinated decent method can optimise the proposed objective function 4 for updating param- eters. Regarding the update of the clustering assignment θ we note that the maximization of L(O, E, w. A) with respect to @ given E and is equivalent to minimizing
p
L( w; A)
Σ - 1,3
Eis) (X,3 -E,s)}
The
The optimal solution to this equation is achieved by assigning each adjusted observation point to the sample mean of the adjusted observation with the smallest squared Euclidean distance (e.g. Lloyd 1982). given by
argming ΣΣΣ w Eig
ni Exi (x,y-Exj)
Enj)
(5)
k=1i=
which can be easily computed by using the standard K-means algorithm. Given clustering assignment and variable weight w, the error matrix E can be updated as
K
argming: Σ
(Xij E --1)) Σ EP((IE,.ll2)))))
(6)
k=1 iEc
where măj (Xix Ei)) and Big is the optimized value of Ei from the previous iteration.
As previously discussed, the variable E plays crucial role in robustness. an observation does not fall into any cluster, vector E;,: elements will cxhibit large values. The penalty Pi Ax) controls the robustnoss of the model. We also note that under Å₁ - 00, E,,- equal a p-dimensional zero vector. As a result, the objective function 6 will be equivalent traditional K K-means clustering.
In the minimization problem 6 as previously mentioned, employ the group lasso penalty and group Antoniadis & Fan 2001). The optimizers are obtained by applying the thresholding function corresponding to the group penalty to X,, pk.:: These two types of group penalty functions correspond to two distinct categorics of thresholding functions: the group lasso penalty function corresponds to the multivariate woft-threshold operator, where E,; can be calculated through
SCAD penally (Tibshirani 1996
A1
E,: = (Xi: (k::)) max 0, |X,,: /
(7)
while the group SCAD penalty function corresponds to the multivariate SCAD-threshold operator introduced by Huang et al. (2012). which E,, can be calculated as:
proposed by Witten 2013
S(z; X1). (a 2) -) Zi.
if || S 2\1. if 2\ < || S ad1, if | ad1-
E,:-
1)S(z: (a (1-)(a11)());; -aaa))
(8)
where zi = X - Pk. and S(2;; X) represents multivariate soft-threshold operator. this paper, we set a as equal to 3.7, as recommended in Fan & Li (0001)
In addition the previously mentioned multivariate soft-threshold operator and multivariate SCAD- threshold operator, The application of the aforementioned methods also be extended to the group

---

Published in Transactions on Machine Learning Research (11/2024)
minimax concave penalty (Huang & Zhang 2010
associated with the multivariate minimax concave penalty- threshold operator. However, due to its similarity to the SCAD penalty, we exclusively employed the SCAD penalty. Moreover, we suggest avoiding applying the group hard penalty [Antoniadis] 1996), as optimizing it is inherently difficult. Furthermore, during our simulation process, we discovered that is easy to converge local minima.
As discussed by Witten (2013) and She & Owen [0011) this formulation establishes connection between the framework of robust M-estimation and the proposed model (6) with clustering centre Pkj for each lo held fixed.
As described above, the matrix of errors acquired from the algorithm is weighted based on the current weights. Therefore, it is necessary to restore the weighted error matrix E to its original unweighted state by dividing the current weight before we go to the next step to calculate the new weight for each variable (if the weight of certain variable is zero, then E.g will divide one instead of zero).
Finally, given E and θ the minimization of L(O, E,w; A) with respect to is equivalent to
Pa(w);da) + 1 (w3) }
argmax WEH
my ( E.5) Σ
(9)
where H {Ws E RR+t ER+, denoted
/ w3 } We show in the Appendix that the final solution
(9) can be
Wi S(QR X O A2)); (SQ (X O, :))))
When applying the 2(ww):2)) = X2|w0/ penalty in 9 S can be represented by the soft-thresholding operator 88
i A2 x2
S(1;12)
if if S A2 if < -2
+ A2
When applying the Pa(wy;Ao) = (SCAD(ww)) penalty in (9,, S can be represented by the SCAD-thresholding operator as
sgn(x)) - >2)+ {(a 1)x- ad2 sign(z)) (a-
| S 2\2 2\₂ |x| a12
((r))))= =
if |0|
The parameter A2 has a crucial role in controlling the sparsity the variable weight vector. A higher value of Åa results increased sparsity of the variable vector u.. Specifically, a higher weight value wy indicates greater importance of the jth variable in the clustering process. Conversely, when wj is equal zero, it significs that the jth variable does not contribute to the clustering.
In order to gain a deeper comprehension of the algorithmic process we have proposed, we summarize the aforementioned procedures in a pseudocode in Algorithm 1. Moreover, the computation complexity of T iteration of Algorithm 1 can be represented (n)tTT)), assuming that the number of iterations in the 4th step (L).
2.4 Adaptation of the tuning parameters by the robust Gap statistics
There are two tuning parameters the proposed algorithm: A and Ag. The tuning parameter aı crucial for finding the error matrix and impacts the detection of outliers, and the tuning parameter d2 plays a crucial role in filtering out variables contributing to the clustering process. We here provide the robust Gap statistics to determine these parameters.
5

---

Published in Transactions on Machine Learning Research (11/2024)
Algorithm Iterative algorithm for ARSK
Initialize the weight vector w as W 1/ vD, set tolerance & and r
Initialize the error matrix by setting Fily of the 80% of data farthest from all data center points to Xiy: the others set as O.
repeat
Run the clustering algorithm on the weighted dataset X;j = (Xis Eijy) by
(Xis
K
X
Σ
x
max Ch
k=1
5: Calculate the new Eig by or S
which incorporates the cluster-specific means PR.:
Inxl Σj Eieg, Xij
6: Until the (6) converges, we can obtain the error matrix E and the cluster arrangement O(++1)
Keep those results for the next variable selection phase. Restore the error matrix by setting:
(r+1) 14
== -Ei ///t"/////"
/v
8: Arrange the cluster (++1)) 1) and compute the Q (,((,,,; 0((+1))) E) for different variables by:
K
0, ( E)
X
Σ
X
Σ
x() 2
71
1=1
k=
iE
(r11)
iECH
(rl1)
where Xij E E Compute new variable weight by:
S(Q) O E.)):
(r+1)
w
v (500 X ji O, E. 22))
10:
r=
11: until convergence of the criterion
E, 1 |/o ((((t)))) w0 (r) < is met.
12: Output: w, E, θ
The Gap statistics were originally proposed in Tibshirani et al. (2001) for selecting the number of clusters in the K-means algorithm. Note that the Gap statistics is constructed of squares, D = Ei= wn(n-1 Σx= E'=, (X,, dataset. However, since D is sensitive to outliers, it is not suitable for selecting the tuning parameters under existence of outliers. Hence, we propose a robust version of Gap statistics by adding the error part to each observation when calculating the D As we mentioned before, the error part can outweigh the influence of the outlier, so that we define Gapx,,,, as
function of the between-cluster sum )² - ΣK " Eives, (X.; X²,)²), for the original
B log( ((0) A
Gapa = -log(DR
ij O, E.g): DR Σ (X.;i
(10)
B
6=1
and the optimal value of X2 and Xı corresponding to the largest value of Gap Note that D A2, is a weighted robust. belween-cluster sum of squares defined in 3 and DR(()) is a version with bth permuted datasets. Here, B denotes the number of permuted datasets generated by randomly selecting from the original dataset. Inereasing the number of permuted datasets improves the accuracy of selecting the true parameters. The between-cluserr sum of squares, taking account of the variable weight, is adopted in Witten & Tibshirani (2010). our version used in 11) can be regarded its robust extension.

---

Published in Transactions on Machine Learning Research (11/2024)
Since it may be computationally intensive to search the optimal value of (1) by a grid search method, we conduct an alternating optimization algorithm for tuning parameter search. Specifically, we first set a suitable value Ao for Å and compute the optimal value X of X2 by maximizing GaP Then, we obtain the optimal value Kj of Xι by maximizing Gapx;.A1- Using this search algorithm instead of the grid search method can save a significant computation time. We also offer pseudocode for this procedure for easy understanding of this search algorithm given in Algorithm 2
Algorithm Selection of (A2, A1) by the robust Gap statistics
1: Set some value for Ao and maximize the following Gap statistics with respect to Ag:
Gapx; log(D ba
1
Σ
log(D
(D
R(b) Ag
B
b=
to obtain the optimal X1
2: Fix A X2 and optimize the following Gap statistics with respect to A
Gapx, A5
1
Σ log D
R(b) Ag
B
to obtain the optimal Xi.
3: Output: the optimal set of tuning parameters, XT)
3 Numerical Studies
3.1
Experimental setup
In this section, explore the ability of the proposed clustering method. We consider that each observation is generated independently from a multivariate normal distribution, given that the observation belongs to cluster k. Specifically, for an observation x, in cluster k, we have 2; ~ N(p Ep), where pkk, E R' denotes the mcan vector for cluster k and Σ, E R2xP represents the covariance matrix. Each element of Pa.: independently sampled from from either U(-6,-3) or U(3, 6). To simulate scenarios involving outliers, we introduce a contaminated error distribution with a multivariate normal mixture, represented as - (H Ep) + =N(us.. + by. Ep) where (l....,p) and π [0, 1] represents the proportion of outliers for each cluster. In the simulation study, we consider two types of Ep. In the first scenario, the assume there is a correlation
assumption is made that all variables are independent. In another scenario, structure among variables. To introduce randomness to the appearance of outliers, b, is generated with a random number from one of two uniform distributions: U3-13, -7) or U37,13). In this setup, we assume that some explanatory variables are significantly associated with the response, while the remaining variables are redundant for variable selection. To mimic this situation, we randomly set g elements of pk:: to non-zero values while the remaining elements are set to zero. Consequently, successful method should accurately identify the q positions where the weights Wy corresponding to the true significant variables are non-zero while ensuring that the weights wy of the remaining variables are set to zero.
In order to evaluate the approach's 's clustering accuracy capability, we employ the clustering error rate (CER) Rand 1971). The CER measures the extent to which the model's predicted partition & for a set of n observations is consistent with the true labels C.
The CER defined as
CER(C,
Elaau - 1o
(11)
where, Ien,e, and la,,, indicate whether the ith and i'th observations belong to the same cluster. Lower CER indicates better performance in clustering accuracy outlier detection.
7

---

Published in Transactions on Machine Learning Research (11/2024)
Furthermore, in order to evaluate each method's proficiency in variable sparsity, apply two criteria for vari- able selection. The variable true positive rate (TPR) indicates the approach's success in finding informative variables. The true negative rate (TNR) represents the successful identification of non-informative variables. If both criteria are closer to 1, it indicates that the model provides a more comprchensive explanation of the structure of the variables predicted in the study.
3.2 Simulation 1: selection of turning parameter
In this subsection explores the new Gap statistics' capability to directly search for the tuning parameters of robustness and variable vector sparsity using Algorithm 2. Initially, we considered clusters, each containing 50 observations, number of variables p 50, with q= 5, and all variables are independent, ie., Ep Ip. A proper of hyperparameters Åι and X2 should make the model accurately identify the structure of the clustering data in different contamination levels for π in {0, 0.1, 0.2)..3) The proposed tuning paramcter search strategy generates both the robustiness parameter A and the variable selection tuning parameter X2 using exponential decay. We consider the number of permuted datasets to be 25 (i.e., B = =22). We employed different thresholding function types, where the first thresholding function is used to control model sparsity, and the second thresholding function is used control model robustness. The tuning parameter search process is repeated 30 times for 30 different datasets, and the resulting table is presented in Table 3.2
threshold type
contamination level
number detected outliers
number of detected
informative variable
soft-soft
π = 0
0 (0)
4.300 (0.000)
0.1
15.56 (6.425)
4.533 (0.618)
π 0.2
28.73 (2.434)
4.200 (0.871)
π 0.3
33.10 (14.76)
13.60 (7.203)
soft-SCAD
= 0
0 (0)
4.400 (0.000)
π = 0.1
14.43 (0.495)
6.310 (1.854)
π 0.2
28.68 (1.319)
1.566 (0.495)
π 0.3
42.36 (2.287)
6.800 (2.916)
SCAD-soft
0
0 (0)
4.545 (0.687)
π = 0.1
14.30 (0.483)
5.400 (0.469)
π 0.2
27.11 (6.166)
5.941 (1.748)
π = 0.3
38.63 (16.63)
9.272 (4.540)
SCAD-SCAD
π = 0 π 0.1
0.090 (0.333)
4.272 (0.881)
13.33 (1.258)
5.833 (3.125)
π 0.2
24.10 (5.065)
4.500 (2.068)
0.3
44.00 (4.2110)
7.187 (1.223)
Table 1: Using the robust Gap statistic to select the optimal tuning parameters for the RSKC algorithm based on soft-soft-thresholding
BrAAD-soffttirrehoodddng and SEs--CAD-thhreee -ireecldDttn we report the average and its standard error (in parentheses) of the number of detected outliers and the number of detected informative variables.
Table 3.2 summarizes the results of the evaluation measures. Overall, a good tuning parameter λι should accurately identify the number of outliers under different contamination levels. As the contamination level æ increases from 0.3, we observe that the four thresholding methods demonstrate varying outlier detection capabilitics across different contamination levels. At = 0.1, the detected number of outlicrs ranges from 13.33 to 15.56, with
achieving the highest average of 15.56 (standard error: 6.425). When increases to 0.2, the detected number of outliers falls between 24.10 and 28.73, with soft-soft and soft-SCAD-thrrsholding showing similar performance. This demonstrates that the Gap statistics can help the ARKC method determine Y1 to some extent. However, when the contamination level reaches 0.3, the soft-soft-hhreehoodiing ARSK method detects 33.10 10 outliers on average, which is lower than the SCAD-SCAD- thresholding ARSK method (44.00 outliers) and deviates from the true number of outliers.
Regarding the selection of the tuning parameter X2, good tuning parameter should accurately identify the number of informative variables under different contamination levels. When the data is clean (i = 0),
8

---

Published in Transactions on Machine Learning Research (11/2024)
all four methods identify around 4 to 5 informative variables, which is close to 5. As the contamination level increases, the number of detected informative variables generally increases, with some fluciuations. At x== 0.3, the soft-soft-thresooiiin ARSK method detects 13.60 informative variables on average, which significantly deviates from the 5. In contrast, the SCAD-SCAD-thresooddnng ARSK method detects 7.187 informative variables, showing better performance in the presence of high contamination. Furthermore, errors in identifying informative variables can lead to significant biases in outlier detection.
Considering these results, can conclude that the proposed parameter selection algorithm based the Gap statistic works relatively well under low to moderate contamination levels. However, its performance may degrade when the contamination level is high, especially for the soft-thresholding-aaee methods. The SCAD-brresholding based methods generally show better robustness and accuracy in identifying outliers and informative variables under various contamination levels.
3.3 Simulation 2: comparison with other methods
In this subsoction, we present the findings of a comprehensive simulation study conducted to investigate the cluster data structure and properties by the ARSK algorithm based on soft-thresholdigg and SCAD- thresholding. The study extensively applied Gap statistics to estimate optimal parameter settings and com- pared with several benchmark approaches, including the original K-means (KC), PCA-K-meat (PCA-KC) a mixed approach that combines Principal Component Analysis (PCA) with the original K-means algo- rithm, as initially proposed by Ding & He (2004). It is noted that PCA is a widely recognized dimensionality reduction technique that transforms nigh-dimensional data into a liver-dimesiionll space. We also consid- ered Trimmed K-mcans (TKM), which exccutes standard K-mcans clustering and subsoquently climinates a predetermined percentage of points exhibiting the greatest Euclidean distance from their respective clus- ter centroids. Furthermore, we employ Robust and Sparse K-means (RSKC), whose fundamental premise involves adapting the SK-means algorithm of Witten & Tibshirani (2010) in conjunction with Trimmed K-means. Lastly, we apply the Weighted Robust K-means (WRCK) method, which shares similarities with RSKC but distinguishes itself through the incorporation of observation-speeifcc weights in its objec- tive function. These weights are reflective of cach observation's isolation level relative to its surrounding neighbourhood.
Considering the mixture error model mentioned at the beginning, we also generated data for 3 clusters, each containing 50 observations. We set the number of variables as p 50, with 9 5, and p = 500, with q = 50. Meanwhile, we consider π in {0, 0.1, 0.2) for specifying the distribution of outliers. As previously stated, we study two types of covariance matrices: one with Σn In and the another generated according to the method proposed by Hirose et al. (2117) ) as
1
Pt
Pt
Σp Q
Pt
Pt
Q
1
Pr
Pt
where the Q denotes a pxp random rotation matrix satisfying Q"' = =Q-1 and the Pr are randomly generated form an Uniform distribution 1)).
In order to evaluate the clustering accuracy and outlier detection capability of cach approach, we employ the clustering error rate (CER) as described above. To assess the outlier detection performance, the outliers identified by each model are assigned to the (K + 1)-th cluster group. Table 3.3 presents the results of 100 simulations for the scenario where all variables are independent, i.c., Σp Ip, while Table 3.3 summarizes the results of 100 simulations for the scenario where variables exhibit a certain correlation structure.
In Tables 3.3 and 3.3] the results show that KC, PCA-KC, and TKM approaches, which lack robustness for high dimensions or outliers, are the worst performers, particularly when the dimensionality p and proportion of outliers increase. For instance, in the independent case with p 500 and = 50, the CER of KC increases from 0.050 to 0.341 increases from to 0.2. Similarly, the CER of PCA-KC and TKM
9

---

Published in Transactions on Machine Learning Research (11/2024)
p 50
p 500
q =
= 50
π
π 0.1
π 0.2
f 0
5 0.1
I = 0.2
KC
0.073 (0.119)
0.191 (0.124)
0.285 (0.107)
0.050 (0.105)
0.228 (0.154)
0.341 (0.145)
PCA-KC
0.109 (0.113)
0.326 (0.167)
0.383 (0.132)
0.058 (0.113)
0.219(0.098)
0.300 123)
TKM (a = 0.1)
0.140 (0.027)
0.093 (0.036)
0.220 (0.087)
0.128 (0.005)
0.098 (0.044)
0.302 (0.146)
RSKC (a = 0.1)
0.119 (0.029)
0.008 (0.031)
0.169 (0.118)
0.105 (0.007)
(0)
0.123 (0.129)
RSKC (a 0.2)
0.208 (0.029)
0.127 (0.015)
0.005 (0.021)
0.187 (0.007)
0.123 (0.009)
0 (0)
WRCK
0,129 (0.037)
0.208 (0.104)
0.136 (0.071)
0.118 (0.024)
0.084 (0.028)
0.046 (0.018)
soft-Soft-ARSK
0,010 (0.043)
0.014 (0.043)
0.032 (0.064)
(0)
0 (0)
0 (0)
Poft-SCAD-Ae
0.010 (0.038)
0.013 (0.035)
0.026 (0.033)
0 (0)
0 (0)
0 (0)
SCAD-soft-ARE ARRK
0.021 (0.059)
0.016 (0.036)
0.031 (0.061)
0 (0)
0 (0)
0 (0)
SCAD-SCAD-ARSS
0.017 (0.057)
0.023 (0.052)
0.019 (0.045)
0 (0)
0 (0)
0 (0)
Table 2: When variables are independent, the average values of the CER and their standard errors (in parentheses) based on 100 Monte Carlo replications.
p 50
500
9 =
q 50
π ==00 0
π = 0.1
π = 0.2
. = 0
π = 0.1
. = 0.2
KC
(.078(0.12))
(..1))))
1.3)
0.037(0.096)
0.219 172)
0.306(0.136)
PCA-KC
0.284(0.169) 0.112(0.113)
0.372(0.126)
0.113(0.131)
0.278(0.016) 0.212(0.075)
TKM (a = 0.1)
(0.85))
0.039(0.080)
0.240(0.104)
0.065(0.002)
0.139 (0.053)
0.136 (0.286)
RSKC (ar 0.1)
0.082 (0.021)
0.005 (0.016)
0.235 (0.128)
0.082 (0.008)
0 (0)
0.071 (0.089)
RSKC (ar = 0.2)
0.194 (0.021)
0.123 (0.021)
0.008 (0.015)
0.162 (0.010)
0.122 (0.009)
(0)
WRCK
0.119 (0.021)
0.203 (0.116)
0.146 (0.071)
0.126 (0.019)
0.102 (0.020)
0.022 (0.081)
soft-soft-ARSK
0.009 (0.038)
0.014 (0.039)
0.037 (0.047)
(0)
(0)
0.000 (0.002)
Soft-SCAD-ARSK
0.005 (0.020)
0.015 (0.039)
0.032 (0.026)
0 (0)
0 (0)
0.001 (0.002)
SCAD-soft-
0.007 (0.021)
0.009 (0.014)
0.047 (0.069)
0 (0)
0 (0) 0.000 (0.002)
SCAD-SCAD-ARSK
0.005 (0.019)
0.001 (0.005)
0.029 (0.033)
0 (0)
( 0.002 (0.018)
Table 3: When variables correlated, the average values of the CER and their standard errors (in paren- theses) based on 100 Monte Carlo replications.
(a = 0.1) increase from 0.058 to 0.300 and 0.128 to 0.302, respectively, under the same scenario. Conversely, RSCK, WRCK, and ARSK perform excoptionally well in high-dimension and high-contamination scenarios. However, the effectiveness of the RSCK approach heavily depends on the choice of the parameter a. When a equals the true outlier proportion, RSCK performs impressively, with CER close to 0. For example, in the correlated case with p- 500 and q = 50, the CER of RSKC (a = 0.1) is when क = 0.1. However, the approach loses effectiveness if is incorrectly specified, leading to higher CER values. In the same scenario, the CER of RSKC (a 0.2) is 0. when = 0.1. The WRCK approach, based on the densily clustering method, generally outperforms RSCK when outliers are present. Still, it may erroneously classify more normal data points as outliers when no contamination exists. In contrast, the ARSK approach, with different combinations of thresholding functions (e.g., soft-soft, Soft-SCAD, SCAD-soft, and SCAD-SCAD). consistently maintains a low Clustering Error Rate (CER) across all simulated data scenarios.
In summary, the simulation results demonstrate the effectiveness and robustness of the ARSK method in clustering digh-dimensional datasets, particularly in the presence of outliers and variable correlations. Across diverse data scenarios, the ARSK variants consistently outperform traditional methods, including KC, PCA- KC, TKM, RSKC, and WRCK, highlighting the distinct advantages of the ARSK approach.
Both In Table 3.3 and Table 3.3 we observe that the CER of KC is worse than ARSK when there are no outliers in the datasct. Morcover, when π 0.1, the CER of the TKM method significantly higher than that of the RSCK method. We attribute this result to the presence of noisy variables, which contribute to the increased CER. Therefore, it is essential to identify these noisy variables to reduce their impact.
10

---

Published in Transactions on Machine Learning Research (11/2024)
Through 100 repetitions, the results of the TPR and TNR are presented in Tables 3.3 and 3.3 As the proportion of outliers increases from 0 to 0.2, both the TPR and TNR of most approaches demonstrate a decreasing trend across different scenarios, suggesting that the presence of outliers can negatively influence the variable selection performance. However, the TPR and TNR of cach type of ARSK approach variant remain relatively high and stable across various data scenarios, demonstrating ARSK's ' superior variable selection capabilities in the presence of outlier contamination in datascts. For example, in the independent case with p- 500 and q 50, the TPR of soft-soft-ARSK, Soft-SCAD-ARSK, SCAD-soft-ARSK, and SCAD-SCAD-ARSK remain above 0.797, 0.810, 0.798, and 0.810, respectively, even when π 0.2. Similarly, their TNR values stay above 0.977, 0.980, 0.978, and 0.981 in the same scenario. These results emphasize the superior ability of soft-ARSK and SCAD-ARSK methods to accurately identify informative and non- informative variables in the presence of outliers, making them valuable tools for handling contaminated clustering datasets.
In contrast, during the simulation process, we observed that traditional trimmed approaches such as WRCK and RSKC often struggled to maintain their variable selection performance under high outlier proportions. In many cases, these methods incorrectly assigned non-zero weights to variables that originally had zero weights, leading to an increased number of false positives. This limitation makes it challenging for WRCK and RSKC to effectively detect informative variables in the presence of outliers, emphasizing the need for more effective methods like A ARSK.
p 50
p 500
9 5
q 50
π
x 0.1
π 0.2
I 0.1
π 0,2
TPR
RSKC (a = 0.1)
0.924 (0.038)
0.890 (0.100)
0.710 (0.267)
0.845 (0.
0.908 (0.100)
0.780 (0.255)
RSKC (a 0.2)
0.654 (0.109)
0.728 (0.143)
0,850 (0.098)
0.690 (0.043)
0.993 (0.027)
0.917 (0.095)
WRCK
0.888 (0.099)
0.664 (0.172)
0,834 (0.106)
0.734 (0.044)
0.978 (0.048)
0.841 (0.101)
soft-soft-ARSK
0,922 (0.127)
0.842 (0.169)
0.810 (0.183)
0.961 (0.035)
0.853 (0.053)
0.797 (0.084)
-of-SSAD-ARSK
0,944 (0.103)
0.910 (0.162)
0.800 (0.171)
0.932 (0.020)
0.836 (0.084)
0.818 (0.018)
SCAD-soft-ARSS
0.962 (0.092)
0.906 (0.153)
0.766 (0.144)
0.964 (0.029)
0.969 (0.024)
0.798 (0.054)
SCAD-SCAD-ARSH
0.970 (0.076)
0.802 (0.150)
0.810 (0.160)
0.969 (0.026)
0.820 (0.054)
0.810 (0.055)
TNR
RSKC (a 0.1)
0.671 (0.427)
0.660 (0.475)
0.926 (0.102)
0.806 (0.377)
0.548 (0.490)
0.817 (0.235)
RSKC (a 0.2)
0.978 (0.142)
0,880 (0.326)
0.753 (0.430)
1 (0)
0.086 (0.259)
0.532 (0.482)
WRCK
0.660 (0.476)
0.911 (0.199)
0.837 (0.309)
1 (0)
0.261 (0.416)
0.536 (0.478)
soft-soft-
(0)
0.969 (0.048)
0.981 (0.042)
0.999 (0.001)
0.986 (0.012)
0.977 (0.009)
-of--SCAI-ARSK
(0)
0.954 (0.040)
0.996 (0.008)
0.972 (0.001)
0.988 (0.014)
0.980 (0.009)
SCAD-soft-ARS
(0)
0.906 (0.153)
0.976 (0.021)
0.999 (0.001)
0.999 (0.001)
0.978 (0.008)
SCAD-SCAD-ARSK
1 (0)
0.976 (0.021)
0.969 (0.017)
0.999 (0.000)
0.981 (0.010) 0.981 (0.010)
Table 4: When variables are independent, the average of 100 tests of TPR and TNP for variable selections and its standard error(in parentheses) are presented for the various data scenarios
3.4 Applications to real data
3.4.1
UCI data application
In this subsection, we apply ARSK with different thresholding configurations to real-life datasets and compare it with the benchmark methods previously mentioned the simulation study. We consider the dataset all from the UCT Machine Learning Repository (Dua & Graff (2019)). The datascts evaluated include various applications, such as glass identification, breast cancer diagnosis, acoustic signal classification, spam base detection, mortality analysis, and Parkinson's 's disease detection. The real-life dataset comprises data with wide range of dimensions, ranging from low to high. As in the previous analysis, we also applied Gap statistics to estimate the hyperparameter. Due to the different units of the continuous variables, we normalize all of continuous variables before conducting the study.
Table summarizes the result of the real-life dataset comparison. The reported values are CERs for cach approach. In order better illustrate the performance of that method, we have emphasized the first and second best CER values in bold when applied to a certain real-life dataset. While traditional methods such
11

---

Published in Transactions on Machine Learning Research (11/2024)
50
p 500
9
50
π7
R 0.1
πT = 0.2
T
= 0.1
π 0.2
TPR
RSKC (a = 0.1)
0.917 (0.013)
0.839 (0.125)
0.761 (0.236)
0.813 (0.125)
0.934 (0.092)
0.802 (0.294)
RSKC (a 0.2)
0.705 (0.132)
0.608 (0.260)
0.860 (0.015)
0.695 (0.232)
0.903 (0.057)
0.884 (0.072)
WRCK
0.813 (0.126)
0.604 (0.212)
0.734 (0.126)
0.613 (0.213)
0.932 (0.074)
0.812 (0.126)
soft-soft-ARSK
0.953 (0.086)
0.860 (0.175)
0.769 (0.042)
0.974 (0.031)
0.843 (0.071)
0.733 (0.095)
soft-SCAD-ARSK
0.940 (0.105)
0.855 (0.181)
0.806 (0.177)
0.968 (0.022)
0.849 (0-071)
0.773 (0.109)
SCAD-soft-RRK
0.928 (0.096)
0.830 (0.177)
0.772 (0.169)
0.973 (0.029)
0.877 (0.051)
0.813 (0.063)
SCAD-SCAD-ARSK
0.928 (0.104)
0.918 (0.120)
0.758 (0.186)
0.968 (0.036)
0.874 (0.053)
0.793 (0.058)
TNR
RSKC (a 0.1)
0.680 (0.387)
0.730 (0.325)
0.864 (0.226)
0.776 (0.218)
0.632 (0.363)
0.743 (0.206)
RSKC (a 0.2)
0.837 (0.232)
0.935 (0.102)
0.794 (0.400)
1 1
0.136 (0.454)
0.542 (0.436)
WRCK
0.697 (0.437)
0.954 (0.083)
0.803 (0.231)
1 (0)
0.476 (0.376)
0.644 (0.306)
soft-soft-RRKK
(0)
0.961 (0.085)
0.985 (0.026)
0.998 (0.006)
0.977 (0.011)
0.984 (0.015)
soft-SCA--RSKK
(0)
0.876 (0.080)
0.997 (0.008)
0.999 (0.001)
0.978 (0.012)
0.991 (0.008)
SCAD-softARSK SCAD-SCAD-ARSK
(0)
0.978 (0.035)
0.992 (0.017)
0.983 (0.010)
0.977 (0.011)
0.983 (0.010)
(0)
0.968 (0.020)
0.994 (0.013)
0.997 (0.005)
0.986 (0.010)
0.981 (0.009)
Table 5: When variables are correlated, the average of 100 tests of TPR and TNP for variable selections and standard error(in parentheses) are presented for the various data scenarios
RSKÖ
ARSK (proposed)
dataset
K
(n.
KC
PCA-KC
TKM
0.1
a = 0.2
WRCK
soft-soft
soft-SCAD
SCAD-soft
SCAD-SCAD
glass
(214,9)
0.334
0.299
0.324
0.324
0.304
0.319
0.259
0.297
0.159
0.142
Breast Cancer
(699.9)
0.080
0.080
0.107
0.071
0.081
0.104
0.065
0.057
0.079
0.084
Acoustic
(400, 50)
0.309
0.443
0.392
0.293
0.280
0.382
0.354
0.360
0.355
0.346
Spanbase
2
(1601, 57)
0.176
0.181
0.310
0.106
0.563
0.493
0.434
0.132
0.361
0.113
mortality
3
(198, 90)
0.227
0.185
0.228
0.153
0.149
0.106
0.131
0.136
0.143
0.149
DARWIN
(174,441)
0.391
0.400
0.447
0.402
0.434
0.465
0.385
0.398
0.454
0.401
Parkinson
(756, 754)
0.477
0.462
0.500
0.482
0.563
0.392
0.445
0.441
0.429
0.462
Table 6: Comparison of CERs of different algorithms for the real-world data. smallest CER values are shown in bold.
The smallest and sccond
as KC, PCA-KC, and TKM demonstrate effectiveness in certain contexts, they generally yield higher CERs across the datasets compared to the proposed ARSK algorithm and its variants. For instance, in the glass and Breast Cancer dataset, the CER values obtained by KC, PCA-KC, and TKM are considerably higher than those of the best-performing ARSK. In seven real-world datasets, our ARSK method ranked among the top two best methods in six of the datasets. Therefore, the proposed ARSK algorithm consistently achieves superior performance with lower CERs across most of the real datascts, particularly when using the multivariate soft-threshold operator to oblain the error matrix E. ARSK demonstrates competitive performance in this regard. Furthermore, Figure 3.4.1 summarizes the number of zero elements in the variable weight vector of cach method, where ARKC adopted a combination of thresholding corresponding data, the variable
to the minimum CER. Figure 3.4.1 shows that when analyzing weight vector predicted by ARSK exhibits strong sparsity.
12

---

Published in Transactions on Machine Learning Research (11/2024)
The number zaro element in variable weight vector
method
ARESI
RSSKC-0.1
RSSKC-22
WRRCX
glsss Conccokcourtc
Spambase monallyy
detaset
Figure The number of zero elements in the variable weight vector of each method, with results ranging from low-dimension to high-dimension dataset
This comparative analysis emphasizes the ARSK algorithm's success in significantly reducing classification crror rates across diverse real-world datasets, showcasing the versatility and adaptability of our proposed clustering method. Its different configurations (e.g., sofl-soft, soft-SCAD, SCAD-soft, and SCAD-SCAD) allow it to handle various data structures and sizes effectively. The consistent performance across different real-world datasets underscores the broad applicability of our approach in various fields.
3.4.2 Analyzing the group structure of recognition of handwritten digits
The dataset was downloaded from Kaggle, a data science competition platform. It consists of digitized images of the numbers O to 9 handwritten by 13 subjects. The images were divided into 64 non-overlapping blocks of x 4 bits. For each block, we observe the number black pixels. Hence, each image is represented by p = 64 variables recording the count of pixels in each block, taking integer values between and 16. Its names identify the true digit in the image, and the column names identify the position of each block in the original bitmap.
Dig stidecc
Figure 2: The images of handwritten digits data
Figure 3:
The heatmap of the variable weight vector obtained by SCAD-soft. ARSK
Our focus is to delect an entire variable weight and group structure, i.e,, ten groups. Given the superior performance of SCAD-soft. ARSK in subsection 3.4.1, we apply RSKC (a = 0.1), RSKC (a = 0.2), WRCK, and SCAD-soft ARSK this dataset. Although the precisc contribution of individual pixels the cluster- ing process remains indeterminate, visual inspection of the real images of handwritten digits data (0 and 9) Figure 2 shows that the pixels situated in both the initial and terminal columns exhibit uniform white
13

---

Published in Transactions on Machine Learning Research (11/2024)
RSKC-0.2
RSKC-0.1
SCAD-soft ARSK
WRCK
CER
0.086
0.077
0.079
0.156
Table 7: Evaluation of the clustering performance recognition of handwritten digils
RSKC-0.2 RSKC-0.1
SCAD-soft ARSK WRCK
CFR
13
Table 8: The number of zero elements within the weight vector derived from each method
colouration. Thus, these 16 pixels can be reasonably deemed redundant and inconsequential for the clus- tering. The optimal robustness and sparsity parameter for ARSK also selected using the Gap statistic. The evaluation of the resulting clustering solution is presented in Table 3.4.2) We also examine the final variable weights obtained by the sparse aaheans-based algorithms. Table 3.4.2 illustrates each method's mumber of zero elements within the weight vector derived. The analysis of the weight values reveals that WRCK is ineffective in achieving sparsity within the variable weight vector, as evidenced by the exclusive presence of nonzero elements throughout. Correspondingly, it also exhibits the highest. CER. As shown in Table 3.4.2 the performance of RSKC and ARSK is similar, with CER values around 0.8. However, the RSKC (a 0.2) produces weight vectors with only zero elements, indicating a modcrate sparsity level. In contrast, the proposed ARSK approach generales weight. vectors with 13 zero elements that are close to 16. We also presented heatmap, in Figure 3 visualization of the variable weight vector derived from the SCAD-soft ARSK algorithm Figure 2 wherein the intensity of the chromatic scale correlates positively with the magnitude of the weights. We can conclude that these findings correspond with those shown in Figure 2 These results demonstrate that ARSK exhibits both spansity and accuracy of feature selection capabilities.
4 Concluding remarks
In this paper, we propose an Adaptively Robust and Sparse K-means Clustering algorithm designed to handle multiple tasks within a unified framework for comprehensive analysis of contaminated high-dimensional datasets. Our approach demonstrates greater fiexibility in handling outliers within real-world datasets. We also proposed a modified Gap statistic with an alternating optimization algorithm to search for luning parameters. This approach significantly reduces the computational time required for parameter tuning compared to traditional methods. The experimental results demonstrate the significant capabilities of the proposed approach in revealing group structures in high-dimensional data. Our model not only performs well in detecting outliers but also identifies significant. variables directly, leading to more thorough understanding of complex clustering datasets.
In our current study, we assume that the number of clusters is known and consistent with the traditional K-means method. Additionally, Chen & Witten (2023) proposed a test for detecting differences in means between each cluster estimated from K-means clustering. This technique can applied to the modified dataset Xij - Eig; where noise variables (i.e., Wj = 0) are removed help determine which clusters should be retained.
Furthermore, have consistently utilized single penalty for outliers in the objective function This implements uniform approach to detecting outliers across various groups. However, given the complexity and diversity of the data, more sophisticated strategy may be required. This might involve using distinct group penalty functions for different groups. Incorporating this approach into unsupervised learning presents challenges in selecting group penalty functions. These extensions are left to our future work.
Acknowledgments
This work is supported by the Japan Society of the Promotion of Science (JSPS KAKENHI) grant numbers, 20H000S0 and 21 H00699.
14

---

Published in Transactions on Machine Learning Research (11/2024)
References
Scandinavian Journal of Statistics, Anestis Antoniadis. Smoothing noisy data with tapered coiflets series.
pp. 313-330, 1996.
Anestis Antoniadis and Jianqing Fan. Regularization of wavelet approximations. Journal of the American Slalislical Associalion, 96(455):933 967, 2001.
Sárka Brodinová, Peter Filzmoser, Thomas Ortner, Christian Breiteneder, and Maia Rohm. Robust and sparse k-means clustering for high-dimensional data. Advances in Data Analysis and Classification, 13: 005-932, 2019.
Yiqun T Chen and Daniela M Witten. Selective inference for k-means clustering. Journal of Machine Learning Research, 24(152):1 41, 2023.
Pictro Coretto and Christian Hennig. Robust improper maximum likelihood: tuning, computation, and comparison with other methods for robust gaussian clustering. Journal of the American Statistical Associalion, 111(516):1644-1666 2016.
Juan Antonio Auesta-Albertos Alfonso Gordaliza, and Carlos Matrán. Trimmed k-means: an attempt robustify quantizers. The Annols of Statistics, 25(2):553355) 576, 1997.
Chris Ding and Xiaofeng He. K-means clustering via principal component analysis. Tn Proceedings of the twenty-fors international conference on Machine learning, pp. 29, 2004.
Dheeru Dua and C Graff. Uci machine learning repository irvine. University of California, 2019.
Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle properties.
Journal of the Americon statistical Association, 96(45)) 1360, 2001.
Edward W Forgy. Cluster analysis of multivariate data: efficiency versus interpretability of classifications. biometrics, 21:768 769, 1965.
Hironori Fujisawa and Shinto Eguchi. Robust estimation in the normal mixture model. Journal of Statistical Planning and Inference, 136(11) 3989 4011, 2006.
Kei Hirose, Hironori Fujisawa, and Jun Sese. Robust sparse gaussian graphical modeling. Journal of Multi- variate Analysis, 161:172-190, 2017.
Jian Huang, Patrick Breheny, and Shuangge Ma. A selective review of group selection in high-dimensional models. Slatistical science: a review journal of the Institude of Mathematical Statistics, 27(4), 2012.
Junzhou Huang and Tong Zhang. The benefit of group sparsity. 2010.
Shota Katayama and Hironori Fujisawa. Sparse and robust linear regression: An optimization algorithm and its statistical properties. Statistica Sinica, pp. 1243 1264, 2017.
Yegor Klochkov, Alexey Kroshuin, and Nikita Zhivotovskiy. Robust k-means clustering for distributions with two moments. The Annals of Statisfics, 29(4))2206-222 2021.
Yumi Kondo, Matias Allibian-Barreaa and Ruben Zamar. RSKC: An R. package for robust and sparse k-means clustering algorithm. Journal of Statistical Software, 72:1-26, 2016.
Stuart Lloyd. Least squares quantization in PCM. IEEE transactions on information theory, 28(2):129-13) 1982.
Antonio Punzo and Paul D McNicholas. Parsimonious mixtures of multivariate contaminated normal distri- butions. Biometrical Journal, 16)6:11073) 2016.
William M Rand. Objective criteria for the evaluation of clustering methods. Journal of the American Statistical association, 66(336):846 850, 1971.
15

---

Published in Transactions on Machine Learning Research (11/2024)
Yiyuan She and Art B Owen. Outlier detection using nonconvex penalized regression. Journal of the
American Stalistical Associalion, 106(494):626 639, 2011.
Shonosuke Sugasawa and Genya Kobayashi. Robust fitting of mixture models using weighted complete estimaling equations. Compulational Statistics Dala Analysis, 174:107526, 2022.
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series Statistionl Methodology, 58(1):267 288, 1996.
Robert Tibshirani, Guenther Walther, and Trevor Hastic. Estimating the number of clusters in a data set via the gap statistic. Journal of the Royal Statistical Society: Series B (Statislical Methodology), 63(2): 411-423, 2001.
Lifeng Wang, Guang Chen, and Hongzhe Li. Group scad regression analysis for microarray time course gene expression data. Bioinformatics, 23(12):1486-144) 1494,, 2007.
Daniela M Witten. Penalized unsupervised learning with outliers. Statistics and its Interface, 6(2):211, 2013.
Daniela M Witten and Robert Tibshirani. A framework for feature selection in clustering. Journal of the American Statistical Association, (09())0) 726, 2010.
Miin-Shen Yang, Chien-Yo Yo Lai, and Chih-Ying Lin. A robust EM clustering algorithm for Gaussian mixture models. Pattern Recognition, 45(11):3950-399) 3961, 2012.
Ming Yuan and Yi Lin. Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society Series B: Statistical Methodology, 68(1):49-67, 2006.
Hui Zou and Trevor Hastic. Addendum: 'regularization and variable selection via the clastic net" [j. stat. soc. ser. b stat. methodol. 67 (2005), no. 2, 301-320; mr2137327]. journal of the royal statistical society, 67(5):768 768, 2010.
A Appendix: Derivation of the solution to (9)
Consider maximizing the following regularized objective function for estimating the weights wy.
p L(w,, = Ew Σ
w
Pa(ly:de) + {
j=1
5=1
where Q, (R O, E.)). Taking the derivative with respect to wy and setting it to zero, we obtain:
Q w - r = O.
When the L1 penalty Pa(Wy; X2) applied, I represents the subgradient of the L1 penalty with respect to wy- Solving the above equation yields
Qi A2. if 0. 1
Ssoft (x; (x;)2) A2)
Wj
Q + d2
if if < -2
When the SCAD penalty Pa A2) is applied, I'y represents the subgradient of the SCAD penalty with respect to wy. Solving the above equation yields
SBCAD (x; d2) w0 = {
sgn( (gn(Q))(11)l - A2)
if [22 < 2A₂
(a (a {(a- 1)Q; adzsign (Q,))
if 2 Σ [0, <
a\х
Qj
if Taa2 < Q1
Projecting vřy onto the parameter space
v /E== W/ 1, the final solution can be rewritten as
wi S(Q5i
(SQ3,
16


Reviews:
[{"review_id": "9q0n5SCQQb", "reviewer": "Reviewer_BVFi", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n\n- It solves a well-defined problem on adaptively robust and sparse clustering.\n- A paper is generally well-written.\n- The proposed method is analyzed appropriately.\n\n### Weaknesses\n\n- Presentation can be improved more.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "- This work proposes a $K$-means clustering algorithm.\n- It proposes adaptively robust and sparse $K$-means clustering.\n- The authors introduce an error component for each observation to enhance robustness.\n- Several experimental results are demonstrated to show the effectiveness of the proposed method.\n- In addition, source code has been provided in the submission.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- What is SCAD? It is not defined in the submission.\n- I think the notation of $x_{i, :}$ is too complex. Why don't you simply use $x_{i}$?\n- $(x\\_{i, :} - x\\_{i', :})^2$ should be the squared value of $L_2$-distance.\n- Could you provide the computational complexity of the iterative algorithm, which is described in Algorithm 1?\n- I think that Sections 2.4 and 3.2 are useful for understanding how to tune $\\lambda_1$ and $\\lambda_2$, but could you provide a practical way to tune them?\n- Could you add some visual example of your algorithm that can show the strengths of the proposed algorithm? It will be helpful to efficiently understand your method.\n- What is $E_{\\cdot j}$ in the last paragraph of Page 4?", "broader_impact_concerns": "I don't have particular concerns on the broader impact of this work, but the authors might be able to discuss them.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "NKb38aYekX", "reviewer": "Reviewer_n9Js", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "*Strength*: I think the main idea is interesting and relatively easy-to-follow. If my interpretation is correct, that's a neat solution to find outliers. \n\n*Weakness*: Text quality could be improved. Specially, the algorithm's presentation and motivation could be more elaborated. Instead of motivating the problem, the introduction reads like a \"related work\" section, which sets your formulation as an incremental combination of the previous idea. The introduction would be stronger if the authors referred to practical problems that now can be solved by the proposed formulation and were not possible before.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a k-means formulation for clustering data points while detecting outliers and unimportant features. The proposed algorithm is called adaptively robust and sparse K-means (ARSK). The basis for their work is the algorithm proposed by Witten & Tibshirani (2010) that rephrases the k-means objective function in terms of the \"between-cluster sum of squares.\" In that paper, the goal was to perform variable selection with k-means. In this paper, the authors go one step further and propose a formulation that is also robust to outliers. The main is to estimate a new parameter `E_{ij}` for each observation `x_{ij}`, where `i` is the index of the observation and `j` is the coordinate. `E_{ij}` will contribute to the objective function when `x_{ij}` is an outlier. The authors propose enforcing sparsity to the matrix `E` so that the optimization can't \"mark\" all points as outliers. The algorithm is validated on synthetic experiments and over UCI datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "My main recommendation is a text revision. There are couple of typos that a minor revision would fix it, for example, capital letters in the middle of the text and phrases that could be shorter and clearer.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "KdBcTEP5J6", "reviewer": "Reviewer_S4PT", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. The proposed algorithm combines robustness and sparsity, enhancing its resistance to outliers and noise, thereby improving the stability and reliability of clustering results. The algorithm is supported by theoretical.\n2. It effectively reduces computational costs and enhances practical utility, with the source code being readily available.\n\nWeaknesses:\n1. The use of datasets from the UCI Machine Learning Repository may not fully represent the diversity and complexity of real-world data. \n2. The recent comparison algorithms are published in 2016 and 2019, which might not provide a comprehensive view of the current state-of-the-art in clustering algorithms. \n3. The comparison algorithms used in the paper lack references to their original papers and necessary explanations. This makes it difficult for readers to understand these algorithms.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents a novel adaptively robust and sparse K-means clustering (ARSK) algorithm, which overcomes the limitations of the conventional K-means algorithm in handling outliers and high-dimensional noisy data. ARSK can improve the robustness of clustering outcomes by incorporating a redundant error component and a group sparse penalty. Furthermore, ARSK reduces the computational time required for parameter tuning compared to traditional methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please refer to the Weaknesses.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers agree the method proposed in this work is of interest and deserve publication at TMLR. However, they were also unanymous in pointing out that the paper would benefit from a throughout revision. During the discussion period, the reviewers made some suggestions that the authors took into account. 

However, in their recommendation reviewer *n9Js* pointed out that the text is still confusing, giving the following examples:

> *"Moreover, if the approach learns to handle unusual data by adjusting outliers instead of removing them, it particularly helps to preserve the continuity of estimating parameters"*

> *"Meanwhile, to reduce the dimension of data, we also introduce a SCAD penalty to weights for each variable to eliminate noise variables irrelevant to clustering. Based on our analysis of real-world datasets, in the ultra-high dimensional setting (i.e., p > 700), the result procured through the implementation of the lasso penalty proposed by Witten & Tibshirani (2010) frequently poor performance. Therefore, we have developed a methodology based on the SCAD penalty. The strong sparsity property of the SCAD penalty effectively improves the sparsity of estimating the variable weight vector in the ultra-high dimensional scenario."*

- p is only defined later.
- donated -> denoted

> *"A proper of hyperparameters $\lambda_{1}$ and $\lambda_{2}$ should enable the model to accurately identify the number of outliers and the number of informative variables in different contamination levels for"*

Following the reviewers recommendation, I am recommending the manuscript for publication to TMLR. However, I do urge the authors to make a throughout revision of the text, including a spelling check to improve readability.
Decision: Accept (poster)