# Tandem clustering with invariant coordinate selection

Andreas Alfons<EMAIL_ADDRESS>Aurore Archimbaud  Klaus Nordhausen  Anne
Ruiz-Gazen  Erasmus School of Economics, Erasmus University Rotterdam,
Netherlands  Department of Mathematics and Statistics, University of
Jyväskylä, Finland  Toulouse School of Economics, Université Toulouse 1
Capitole, France

###### Abstract

For multivariate data with noise variables, tandem clustering is a well-known
technique that aims to improve cluster identification by first reducing the
dimension. However, the usual approach using principal component analysis
(PCA) has been criticized for focusing only on inertia so that the first
components do not necessarily retain the structure of interest for clustering.
To overcome this drawback, a new tandem clustering approach based on invariant
coordinate selection (ICS) is proposed. By jointly diagonalizing two scatter
matrices, ICS is designed to find structure in the data while returning affine
invariant components. Some theoretical results have already been derived and
guarantee that under some elliptical mixture models, the group structure can
be highlighted on a subset of the first and/or last components. Nevertheless,
ICS has received little attention in a clustering context. Two challenges are
the choice of the pair of scatter matrices and the selection of the components
to retain. For clustering purposes, it is demonstrated that the best scatter
pairs consist of one scatter matrix that captures the within-cluster structure
and another that captures the global structure. For the former, local shape or
pairwise scatters are of great interest, as is the minimum covariance
determinant (MCD) estimator based on a carefully selected subset size that is
smaller than usual. The performance of ICS as a dimension reduction method is
evaluated in terms of preserving the cluster structure present in data. In an
extensive simulation study and in empirical applications with benchmark data
sets, different combinations of scatter matrices as well as component
selection criteria are compared in situations with and without outliers.
Overall, the new approach of tandem clustering with ICS shows promising
results and clearly outperforms the approach with PCA.

###### keywords:

$k$means, Principal Component Analysis, Linear Discriminant Analysis, Scatter
Matrices, Minimum Covariance Determinant, Robustness

††journal: Econometrics and Statistics

##  1 Introduction

Clustering or cluster analysis is the task of finding groups of observations
in a data set. Clustering problems arise in many application fields such as
marketing with customer segmentation or biostatistics with gene expression
data analysis. In the present paper we mainly focus on hard clustering where
the observations are divided into subsets called clusters. The clusters form a
partition of the data set and every observation belongs to one and only one
cluster. As stated for example in Hennig et al. [2015, p. 2], the definition
of a cluster is not unique nor precise. Thus, many clustering algorithms have
been developed and studied, and their usefulness and relevance depend on the
context, the type of data, and the objectives of the data analyst. Several
books and surveys exist on the subject, e.g., Kaufman & Rousseeuw [1990], Xu &
Wunsch [2005], Aggarwal & Reddy [2014], and Hennig et al. [2015].

The generally accepted objective for clustering is that similar observations
should be gathered in the same cluster while dissimilar ones should belong to
different clusters. Among the many existing approaches, common clustering
methods are either distance-based, model-based, or density-based. In the
present paper, we focus on centroid-based clustering techniques that belong to
the family of distance-based methods, as well as clustering via Gaussian
mixture models that belong to the family of model-based methods. The data
under consideration are observations characterized by numerical variables or
features.

When the cluster structure lies in a low-dimensional latent subspace, or when
some variables are noisy and not relevant for finding the cluster structure, a
preprocessing step such as feature selection or dimension reduction can be
helpful. Tandem analysis, as called by Arabie & Hubert [1994], is a well-known
method that consists of carrying out principal component analysis (PCA)
followed by clustering on the selected principal components. However, several
authors including Chang [1983], Soete & Carroll [1994], Vichi & Kiers [2001],
and Radojičić et al. [2021] warn against this approach, as it may mask the
cluster structure of the data. The explanation is that PCA focuses on the
conservation of inertia, which is linked to the variability of the data and
does not necessarily extract the relevant information for clustering. Many
papers propose alternatives to the tandem analysis (see, e.g., Ding & Li,
2007, Vichi & Saporta, 2009, Timmerman et al., 2013, van de Velden et al.,
2019a; and references therein).

In the present paper, we propose a new tandem approach where PCA is replaced
by invariant coordinate selection (ICS) introduced by Tyler et al. [2009]. ICS
is a dimension reduction method in the same spirit as PCA with computation and
selection of linear combinations of the original features. While PCA can be
derived from the eigendecomposition of a unique scatter matrix, ICS is based
on the simultaneous diagonalization of two scatter matrices [see, e.g.,
Nordhausen & Ruiz-Gazen, 2022]. Contrary to PCA, ICS is able to extract
components that are relevant for clustering. More precisely, Tyler et al.
[2009] show that for a location mixture of elliptical distributions and for
any pair of scatter matrices, ICS is able to find Fisher’s linear discriminant
subspace without knowing the class labels. Moreover, ICS is affine invariant
while PCA is only orthogonally invariant.

The use of ICS as a preprocessing step for clustering is not new. ICS and its
earlier versions are considered as a method for dimension reduction prior to
group identification in Art et al. [1982], Caussinus & Ruiz [1990],Caussinus &
Ruiz-Gazen [1993, 1995, 2007], Gnanadesikan et al. [1993], Tyler et al.
[2009], Hennig [2009], Peña et al. [2010], Fekri & Ruiz-Gazen [2015],
Alashwali & Kent [2016], and Fischer et al. [2017, 2020]. However, there is no
in-depth study of ICS in combination with clustering and the present paper
aims to fill this gap. Following the study by Archimbaud et al. [2018] on ICS
for outlier detection, we first highlight the relevance of ICS in the presence
of noise. We then compare several pairs of scatter matrices and discuss the
selection of invariant coordinates. It should be noted that ICS, just like the
standard approach for PCA, is not well-suited to high-dimensional data. In
this article, we assume that the number of observations is large compared to
the number of variables and that the number of variables is moderate.

In the clustering framework, ICS can be compared to linear discriminant
analysis (LDA) but in an unsupervised situation where the groups are unknown.
When the groups are known, and assuming that the clusters have similar shape,
it is possible to compute a within-group and a between-group scatter matrix,
and LDA [see Mardia et al., 1982] is nothing else than a joint diagonalization
of these particular scatter matrices (which is equivalent to the joint
diagonalization of the within-group and total covariance matrices). Keeping
the analogy with LDA, when the groups are unknown, we can expect ICS to give
interesting results when preprocessing data prior to clustering if we jointly
diagonalize a scatter matrix close to the within-group covariance matrix and
another scatter close to the total covariance matrix. Concerning the within-
group matrix, in this paper we further consider local scatter matrices
[Hennig, 2009], weighted pairwise scatter matrices [Caussinus & Ruiz, 1990,
Tyler et al., 2009], and minimum covariance determinant [Rousseeuw, 1985]
scatter matrices based on a small percentage of observations. Interestingly,
Kettenring [2006] claims that cluster analysis and LDA are “at the opposite
ends of a spectrum”. While no information is available about the number or the
types of clusters in clustering, the number of groups and their content is
known in LDA. If we follow this idea, ICS is in the middle of the spectrum
with some assumptions on the type of clusters (ideally with Gaussian
distributions with different means and the same covariance matrices), but
without information on the number of groups and on how the observations are
distributed among the groups. As for LDA, some theoretical properties of ICS
can be derived for mixtures of Gaussian distributions with different means and
equal covariance matrices [see Tyler et al., 2009, Archimbaud et al., 2018],
and we consider this particular model in a privileged manner throughout the
paper. However, theoretical properties of ICS can also be derived for other
types of mixtures [see Archimbaud, 2018], and in practice, ICS can be useful
beyond this particular framework, as we illustrate in artificial and empirical
examples.

The paper is organized as follows. In Section 2, we consider a mixture of
Gaussian distributions with different means and equal covariance matrix. We
give a property stating that, if the number of noise variables increases while
the number of groups is kept fixed, the noise variables will completely mask
the clusters when using the total or the within-group Mahalanobis distances,
while ICS allows for solving the problem. Section 3 gives a detailed
presentation of ICS together with the definition of several scatter matrices
that are relevant in a clustering context, as well as a description of
component selection criteria. Section 4 provides results from an extensive
simulation study for Gaussian mixtures with two, three and five groups. We
compare a large number of dimension reduction strategies using ICS or PCA
followed by some centroid-based and model-based clustering algorithms. Section
5 illustrates the use of tandem clustering with ICS in three empirical
examples, and Section 6 concludes the paper and discusses further
perspectives.

##  2 Highlighting the relevance of ICS in the presence of noise

Consider the following probability model introduced in Caussinus & Ruiz-Gazen
[1993]. Let $\mathbf{X}$ be a $d$-variate real random vector distributed
according to the following mixture distribution:

$\int{\cal N}_{d}(\boldsymbol{x},\boldsymbol{\Sigma}_{W})\
dP(\boldsymbol{x}),$ (1)

where $P$ is concentrated on a unknown $l$-dimensional linear manifold with
$l<d$, and ${\cal N}(\boldsymbol{x},\boldsymbol{\Sigma}_{W})$ denotes the
$d$-dimensional normal distribution with mean $\boldsymbol{x}$ and covariance
$\boldsymbol{\Sigma}_{W}$. The mixing probability $P$ corresponds to the
unknown structure of interest, which is said to be of dimension $l$, while the
Gaussian distribution with unknown covariance matrix $\boldsymbol{\Sigma}_{W}$
corresponds to noise. One example of such a model is the case where $P$ is a
mixture of Dirac probability measures at points $\boldsymbol{\mu}_{1}$,…,
$\boldsymbol{\mu}_{q}$. This example is a finite mixture of Gaussian
distributions with different means and equal covariance matrices, and $l$ is
equal to the dimension of the subspace spanned by the mean vectors
$\boldsymbol{\mu}_{1},\dots,\boldsymbol{\mu}_{q}$. In this paper, we focus on
the case where $l$ is small compared to $d$. We derive a result for the case
of a Gaussian mixture distribution with different means and equal covariance
matrices when considering Mahalanobis distances. The result is analogous to
Proposition 1 in Archimbaud et al. [2018] but adapted to the context of
clustering instead of outlier detection. The number of groups is fixed while
the dimension of the data set increases to infinity, which means that the
intrinsic dimension $l$ is small compared to the original dimension $d$ in
model (1). In such a context, ICS is particularly suitable for data reduction
before clustering since it may allow to get rid of the Gaussian noise.

The question of using the Mahalanobis distance instead of the Euclidean
distance when using the $k$means algorithm has been studied in the literature
for a long time. Gnanadesikan et al. [1993] give an example with three groups
in three dimensions where the Euclidean distance is not appropriate, contrary
to the Mahalanobis distance based on the within-group dispersion, which takes
into account the elliptical nature of the groups (i.e., correlations). Another
advantage of using the Mahalanobis distance is that it is invariant under
affine transformations of the data, which means that there is no need to
standardize the variables prior to $k$means [see Cerioli, 2005, for more
details]. The use of the Mahalanobis distance is however far from being the
most common choice when using $k$means. The main reason is that it is not easy
to adapt the $k$means algorithm to take into account the within-group
covariance matrices (see Cerioli, 2005, Lapidot, 2018; and references
therein). In the tandem approach that we propose, Euclidean distances computed
on the ICS components are Mahalanobis distances in the sense of one of the
scatter matrices on the original variables. Choosing a scatter matrix close to
the within-group scatter matrix allows to take into account elliptical
clusters without modifying the $k$means algorithm.

Let $\mathbf{X}$ be a $d$-variate real random vector distributed according to
a mixture of $q$ Gaussian distributions with $q<d$, different location
parameters $\boldsymbol{\mu}_{1},\dots,\boldsymbol{\mu}_{q}$, and the same
positive definite covariance matrix $\boldsymbol{\Sigma}_{W}$:

${\mathbf{X}}\sim\sum_{h=1}^{q}\epsilon_{h}\,{\cal
N}(\boldsymbol{\mu}_{h},\boldsymbol{\Sigma}_{W}),$ (2)

where $\epsilon_{1},\dots,\epsilon_{q}$ are mixture weights with
$\epsilon_{1}+\cdots+\epsilon_{q}=1$. Let $\boldsymbol{\Sigma}$ denote the
total covariance matrix of $\mathbf{X}$. Let us consider the following squared
pairwise Mahalanobis distances between two random vectors $\mathbf{X}$ and
$\mathbf{Y}$ with the Gaussian mixture distribution from (2):

$\displaystyle r^{2}(\mathbf{X},\mathbf{Y})$
$\displaystyle=(\mathbf{X}-\mathbf{Y})^{\top}{\boldsymbol{\Sigma}}^{-1}(\mathbf{X}-\mathbf{Y}),$
$\displaystyle r_{\boldsymbol{W}}^{2}(\mathbf{X},\mathbf{Y})$
$\displaystyle=(\mathbf{X}-\mathbf{Y})^{\top}{\boldsymbol{\Sigma}_{W}}^{-1}(\mathbf{X}-\mathbf{Y}).$

These distances are affine invariant in the sense that
$r^{2}(\mathbf{A}\mathbf{X}+\mathbf{b},\mathbf{A}\mathbf{Y}+\mathbf{b})=r^{2}(\mathbf{X},\mathbf{Y})$,
for any full rank $d\times d$ matrix $\mathbf{A}$ and any $d$-dimensional
vector $\mathbf{b}$.

In order to identify groups, we should observe that the distance between two
independent random vectors $\mathbf{X}_{h}$ and $\mathbf{Y}_{h^{\prime}}$ from
two different groups $h\neq h^{\prime}$ is larger than the distance between
two independent random vectors $\mathbf{X}_{h}$ and $\mathbf{Y}_{h}$ from the
same group $h$. However, we can prove that asymptotically the variance of the
difference between these distances increases with the dimension $d$ if the
number of groups $q$ is kept fixed, which makes cluster identification
increasingly difficult.

###### Proposition 1

Let $\mathbf{X}_{h},\mathbf{Y}_{h}\sim{\cal
N}(\boldsymbol{\mu}_{h},\boldsymbol{\Sigma}_{W})$ and
$\mathbf{Y}_{h^{\prime}}\sim{\cal
N}(\boldsymbol{\mu}_{h^{\prime}},\boldsymbol{\Sigma}_{W})$ be independent
random vectors with $h,h^{\prime}\in\\{1,\ldots,q\\}$ and $h\neq h^{\prime}$.
Assume that $q>1$ is fixed and that

$\sum_{h=1}^{q}\epsilon_{h}(\boldsymbol{\mu}_{h}-\boldsymbol{\mu})^{\prime}\boldsymbol{\Sigma}_{W}^{-1}(\boldsymbol{\mu}_{h}-\boldsymbol{\mu})=o(d).$

Then under model (2),

$\frac{c}{\sqrt{d}}\Big{[}\left(r^{2}(\mathbf{X}_{h},\mathbf{Y}_{h^{\prime}})-r^{2}(\mathbf{X}_{h},{\mathbf{Y}}_{h})\right)-\mathbb{E}\left(r^{2}(\mathbf{X}_{h},\mathbf{Y}_{h^{\prime}})-r^{2}(\mathbf{X}_{h},{\mathbf{Y}}_{h})\right)\Big{]}$

and

$\frac{c_{\boldsymbol{W}}}{\sqrt{d}}\Big{[}\left(r^{2}_{\boldsymbol{W}}(\mathbf{X}_{h},\mathbf{Y}_{h^{\prime}})-r^{2}_{\boldsymbol{W}}(\mathbf{X}_{h},{\mathbf{Y}}_{h})\right)-\mathbb{E}\left(r^{2}_{\boldsymbol{W}}(\mathbf{X}_{h},\mathbf{Y}_{h^{\prime}})-r^{2}_{\boldsymbol{W}}(\mathbf{X}_{h},{\mathbf{Y}}_{h})\right)\Big{]}$

converge in distribution to standard Gaussian distributions when $d$ goes to
infinity with some constants $c$ and $c_{\boldsymbol{W}}>0$. Moreover, the
expectations of the differences of the Mahalanobis distances remain bounded
when $d$ goes to infinity.

The proof is similar to the proof of Proposition 1 in Archimbaud et al. [2018]
and is not detailed further. The key element of the proof is that the
dimension $l\leq q-1$ of the subspace of interest is fixed (where $q$ denotes
the number of groups) while the total number of dimensions $d$ is growing to
infinity, which means that the number of noise dimensions $d-l\geq d-q+1$ is
also growing to infinity. In the noise dimensions, the expectation of the
difference of distances does not differ but its variance inflates and masks
the difference in expectation in the subspace of interest. If the relevant
subspace of dimension $l$ is known, it is easy to avoid the problem of the
$d-l$ noise dimensions by projecting the data set on this subspace and
calculating a distance based only on the relevant $l$ dimensions. This is
exactly what ICS is all about, providing the data-analyst with the ability to
select a subspace revealing the clusters in an unsupervised way, and to
project the data on this subspace.

For illustration of the above discussion, Figure 1 demonstrates that dimension
reduction can already be useful for $k$means clustering in two dimensions. It
shows an example of two groups, with 850 observations generated from the
distribution $\mathcal{N}_{d=2}(0,\mathbf{I}_{d})$ and 150 observations from
$\mathcal{N}_{d=2}(\mathbf{(}10,0)^{\top},\mathbf{I}_{d})$. The difference
between the groups lies in the first dimension while the second dimension
contains only noise.

Figure 1: Illustration of $k$means clustering on simulated data. The columns
correspond to the standardized data (left), the PCA components (center), and
the ICS components (left). In the top row, the respective scatterplots are
shown. In the middle row, pairwise Euclidean distances are visualized. In the
bottom row, scatterplots of standardized data are shown, with colors
indicating the clusters obtained by $k$means on the initial data and on the
reduced data from PCA and ICS.

We compare the results of the usual $k$means algorithm when applied to (i) the
standardized data (variables Std $X_{1}$ and Std $X_{2}$), (ii) the first
component of PCA ($\text{PC}_{1}$) using the sample correlation matrix, and
(iii) the first component of ICS ($\text{IC}_{1}$) using the
$\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$ pair of scatter matrices
(see Section 3.1). Figure 1 consists of 9 plots. The first row shows
scatterplots of the standardized data (left), the PCA components (center), and
the ICS components (right). While ICS is able to recover the dimension of
interest on the first component, this is not the case for PCA. The middle row
of Figure 1 visualizes the pairwise Euclidean distances between the
observations. For the standardized data (left), it is not possible to
distinguish the distances between observations from the same group and
observations from different groups, which is similar to the result from
Proposition 1. For PCA and ICS, the pairwise distances are computed on
$\text{PC}_{1}$ (center) and $\text{IC}_{1}$ (right), respectively. Only for
$\text{IC}_{1}$, we observe a gap between the two types of distances. The
bottom row of Figure 1 contains scatterplots of the standardized data with
observations colored according to the different $k$means clusterings. The two
groups are not recovered when $k$means is applied to the original standardized
data (left) or $\text{PC}_{1}$ (center), while they are recovered correctly
for $\text{IC}_{1}$ (right).

##  3 Invariant coordinate selection

In Section 3.1, the principle of ICS and various scatter matrices are
detailed. Section 3.2 discusses ICS in a clustering context, while Section 3.3
presents criteria for selecting which invariant coordinates to keep.

###  3.1 Scatter matrices and the principle of ICS

ICS was first introduced as generalized PCA in Caussinus & Ruiz [1990],
Caussinus & Ruiz-Gazen [1993, 1995, 2007], Caussinus et al. [2003], and its
current name was introduced in Tyler et al. [2009], where also many of its
properties were derived. The main ingredients of ICS are two scatter matrices.
A scatter matrix in this context is any $d\times d$ matrix-valued functional
$\boldsymbol{V}(F_{\boldsymbol{X}})=\boldsymbol{V}(\boldsymbol{X})$ which is
affine equivariant in the sense that

$\boldsymbol{V}(F_{\boldsymbol{A}\boldsymbol{X}+\boldsymbol{b}})=\boldsymbol{A}\boldsymbol{V}(F_{\boldsymbol{X}})\boldsymbol{A}^{\top}$

for all $d$-dimensional vectors $\boldsymbol{b}$ and all full rank $d\times d$
matrices $\boldsymbol{A}$. Here $\boldsymbol{X}$ denotes a $d$-variate random
vector with cumulative distribution function (CDF) $F_{\boldsymbol{X}}$. For a
sample
$\boldsymbol{X}_{n}=(\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})^{\top}$,
the scatter functional is applied to the empirical CDF
$F_{\boldsymbol{X}_{n}}$ and called a scatter matrix.

The literature is full of suggestions of scatter matrices, with the most
popular being the sample covariance matrix

$\mathrm{COV}(\boldsymbol{X}_{n})=\frac{1}{n-1}\sum_{i=1}^{n}(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}}_{n})(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}}_{n})^{\top},$

where $\bar{\boldsymbol{x}}_{n}=\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{x}_{i}$
is the sample mean. Other popular proposals come from the robustness
literature and include for example the M-estimator of scatter based on the
likelihood of the Cauchy distribution [Kent & Tyler, 1991], denoted by
$\mathrm{MLC}(\boldsymbol{X}_{n})$. It can be computed using the simultaneous
estimation equations

$\displaystyle\boldsymbol{m}_{C}(\boldsymbol{X}_{n})$
$\displaystyle=\sum_{i=1}^{n}\frac{d+1}{r^{2}_{C}(\boldsymbol{x}_{i})+1}\boldsymbol{x}_{i}\bigg{/}\sum_{i=1}^{n}\frac{d+1}{r^{2}_{C}(\boldsymbol{x}_{i})+1},$
$\displaystyle\mathrm{MLC}(\boldsymbol{X}_{n})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{d+1}{r^{2}_{C}(\boldsymbol{x}_{i})+1}(\boldsymbol{x}_{i}-\boldsymbol{m}_{C}(\boldsymbol{X}_{n}))(\boldsymbol{X}_{i}-\boldsymbol{m}_{C}(\boldsymbol{X}_{n}))^{\top},$

where
$r^{2}_{C}(\boldsymbol{x}_{i})=(\boldsymbol{x}_{i}-\boldsymbol{m}_{C}(\boldsymbol{X}_{n}))^{\top}\mathrm{MLC}(\boldsymbol{X}_{n})^{-1}(\boldsymbol{x}_{i}-\boldsymbol{m}_{C}(\boldsymbol{X}_{n}))$
and $\boldsymbol{m}_{C}(\boldsymbol{X}_{n})$ is the corresponding location
estimate.

One of the most widely-used robust scatter matrices is the minimum covariance
determinant estimator ($\mathrm{MCD}$) [Rousseeuw, 1985]. For a tuning
parameter $\alpha$, the $\mathrm{MCD}$ selects out of the $n$ observations
those $n_{\alpha}=\lceil\alpha n\rceil$ observations
$\boldsymbol{x}_{i_{1}},\ldots,\boldsymbol{x}_{i_{n_{\alpha}}}$ for which the
sample covariance matrix has the smallest determinant, i.e.,

$\mathrm{MCD}_{\alpha}(\boldsymbol{X}_{n})=c_{\alpha}\frac{1}{n_{\alpha}}\sum_{j=1}^{n_{\alpha}}(\boldsymbol{x}_{i_{j}}-\bar{\boldsymbol{x}}_{\alpha,n})(\boldsymbol{x}_{i_{j}}-\bar{\boldsymbol{x}}_{\alpha,n})^{\top},$

where $\bar{\boldsymbol{x}}_{\alpha,n}$ is the sample mean of the selected set
of observations and $c_{\alpha}$ is a consistency factor. In the robustness
literature, we typically have $\alpha\in[0.5,1]$. The $\mathrm{MCD}_{\alpha}$
is sometimes also called the raw $\mathrm{MCD}$. For a value of $\alpha$ close
to 0.5, it has a high breakdown point but a low efficiency, and therefore it
is often combined with a reweighting step to increase efficiency. To avoid
confusion, we call the reweighted version $\mathrm{RMCD}_{\alpha}$. For more
details about $\mathrm{MCD}_{\alpha}$ and $\mathrm{RMCD}_{\alpha}$, see for
example Hubert & Debruyne [2010], Cator & Lopuhaä [2012], and Hubert et al.
[2018].

There are, however, scatter matrices for which robustness is not the focus,
such as the matrix of fourth moments

$\mathop{\mathrm{COV}_{4}}\nolimits(\boldsymbol{X}_{n})=\frac{1}{n}\frac{1}{d+2}\sum_{i=1}^{n}r^{2}(\boldsymbol{x}_{i})(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}}_{n})(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}}_{n})^{\top},$

where
$r^{2}(\boldsymbol{x}_{i})=(\boldsymbol{x}_{i}-\bar{\mathbf{x}}_{n})^{\top}\mathrm{COV}(\boldsymbol{X}_{n})^{-1}(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}}_{n})$
is the squared Mahalanobis distance.

Some other scatter matrices were suggested specifically in the context of ICS.
For a positive, decreasing weight function $w$ and a tuning constant
$\beta>0$, a one-step M-estimator of scatter can be computed as

$\mathrm{SCOV}_{\beta}(\boldsymbol{X}_{n})=\frac{\sum_{i=1}^{n}w(\beta\,r^{2}(\boldsymbol{x}_{i}))(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}}_{n})(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}}_{n})^{\top}}{\sum_{i=1}^{n}w(\beta\,r^{2}(\boldsymbol{x}_{i}))},$

where $r^{2}(\boldsymbol{x}_{i})$ is again the squared Mahalanobis distance.
Closely related to $\mathrm{SCOV}_{\beta}(\boldsymbol{X}_{n})$ is the
estimator

$\mathrm{TCOV}_{\beta}(\boldsymbol{X}_{n})=\frac{\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}w(\beta\,r^{2}(\boldsymbol{x}_{i},\boldsymbol{x}_{j}))(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})^{\top}}{\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}w(\beta
r^{2}(\boldsymbol{x}_{i},\boldsymbol{x}_{j}))},$

where
$r^{2}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})^{\top}\mathrm{COV}(\boldsymbol{X}_{n})^{-1}(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})$.
Note that $\mathrm{TCOV}_{\beta}(\boldsymbol{X}_{n})$ is based on pairwise
differences and therefore does not require a location estimate. Both
estimators are for example studied in great detail in Caussinus & Ruiz-Gazen
[1993, 1995]. Ruiz-Gazen [1996] suggests in this context to consider also

$\mathrm{UCOV}_{\beta}(\boldsymbol{X}_{n})=(\mathrm{SCOV}_{\beta}(\boldsymbol{X}_{n})^{-1}-\beta\,\mathrm{COV}(\boldsymbol{X}_{n})^{-1})^{-1}.$

In the following, we use the weight function $w(x)=\exp(-x/2)$ for these
scatters, and we set $\beta=2$ for $\mathrm{TCOV}$ and $\beta=0.2$ for
$\mathrm{UCOV}$.

The last scatter matrix we introduce is the local shape matrix
$\mathrm{LCOV}_{\boldsymbol{V}_{0},\beta}$ that was suggested by Hennig
[2009]. For $\mathrm{LCOV}_{\boldsymbol{V}_{0},\beta}$, first a scatter
$\boldsymbol{V}_{0}$ is computed and the distance matrix of the corresponding
pairwise Mahalanobis distances are calculated. Then for each observation, the
sample covariance matrix is computed based on the $n_{\beta}=\lceil\beta
n\rceil$ nearest neighbors. All $n$ covariance matrices are subsequently
averaged, for which they first need to be appropriately scaled to be
comparable. We standardize them in such a way that they all have determinant
one. Note that this is not a scatter in our original sense but only a shape
matrix [Nordhausen & Tyler, 2015], which suffices for ICS. We use
$\mathrm{LCOV}_{\mathrm{COV},0.1}$, but one could also use
$\mathrm{LCOV}_{\mathrm{MCD}_{0.8},0.1}$ as in Hennig [2009].

In a clustering framework, $\mathrm{TCOV}$ and $\mathrm{LCOV}$ can be seen as
estimators of the within-group covariance matrix (cf. Figure 12). An
alternative estimator proposed by Art et al. [1982] and Gnanadesikan et al.
[1993] is based on the smallest within-group pairwise differences of
observations. However, it will not be considered further, as the algorithm is
iterative and involves a graphical aid to determine the number of pairwise
differences to consider.

A general result is that all scatter matrices are proportional to each other
(if they exist) if the data comes from an exchangeable sign-symmetric
distribution, which includes all elliptical distributions [Nordhausen & Tyler,
2015]. This means that for such models, when whitening the data [e.g.,
Nordhausen et al., 2008] or performing principal component analysis (PCA), it
does not matter which scatter matrix and location estimate are used, as the
results differ only in scale. Similarly, Tyler [2010] showed that a comparison
of two scatter matrices is only meaningful if the sample size $n$ is larger
than the dimension $d$, hence this is a requirement for ICS.

Given a sufficient sample size, in a clustering framework assuming, for
example, a Gaussian mixture model, different scatter matrices will measure
different population quantities and the idea of ICS is to exploit these
differences. An ICS matrix $\boldsymbol{W}(\boldsymbol{X}_{n})$ is a matrix
that simultaneously diagonalizes two scatter matrices
$\boldsymbol{V}_{1}(\boldsymbol{X}_{n})$ and
$\boldsymbol{V}_{2}(\boldsymbol{X}_{n})$ in the following way:

$\boldsymbol{W}(\boldsymbol{X}_{n})\boldsymbol{V}_{1}(\boldsymbol{X}_{n})\boldsymbol{W}(\boldsymbol{X}_{n})^{\top}=\boldsymbol{I}_{d}\hskip
8.00003pt\text{and}\hskip
8.00003pt\boldsymbol{W}(\boldsymbol{X}_{n})\mathbf{\boldsymbol{V}}_{2}(\boldsymbol{X}_{n})\boldsymbol{W}(\boldsymbol{X}_{n})^{\top}=\boldsymbol{D}(\boldsymbol{X}_{n}),$

where $\boldsymbol{D}(\boldsymbol{X}_{n})$ is a diagonal matrix with diagonal
elements $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{d}$. ICS can be
interpreted as a transformation where $\boldsymbol{X}_{n}$ is whitened with
respect to $\mathbf{V}_{1}$ and then a subsequent PCA is performed with
respect to $\boldsymbol{V}_{2}$ to see what structure is still left after the
whitening with $\boldsymbol{V}_{1}$. ICS is usually computed via a generalized
eigenvalue decomposition [see Tyler et al., 2009].

The so-called invariant coordinates are then obtained as

$\boldsymbol{Z}_{n}=(\boldsymbol{X}_{n}-\mathbf{1}_{n}\boldsymbol{T}(\boldsymbol{X}_{n})^{\top})\boldsymbol{W}({\boldsymbol{X}_{n}})^{\top},$

where $\mathbf{1}_{n}$ denotes an $n$-variate vector full of ones and
$\boldsymbol{T}(\boldsymbol{X}_{n})$ denotes a location estimate, usually the
one that goes along with $\boldsymbol{V}_{1}$. In case $\boldsymbol{V}_{1}$
has no natural accompanying location, any location estimate can be used for
the purpose of centering, or centering can be omitted when computing the ICs.
The name of ICS is motivated by a property of the (centered) invariant
coordinates, which are often also called invariant components or simply ICs.
These ICs are invariant in the sense that for any linear transformation of
$\boldsymbol{X}_{n}$, i.e.,
$\boldsymbol{X}_{n}^{*}=\mathbf{\boldsymbol{X}}_{n}\mathbf{A}^{\top}+\mathbf{1}_{n}\mathbf{b}^{\top}$,
it holds that

$(\boldsymbol{X}_{n}^{*}-\mathbf{1}_{n}\boldsymbol{T}(\boldsymbol{X}_{n}^{*})^{\top})\boldsymbol{W}(\boldsymbol{X}_{n}^{*})^{\top}=(\boldsymbol{X}_{n}-\mathbf{1}_{n}\mathbf{T}(\boldsymbol{X}_{n})^{\top})\boldsymbol{W}({\boldsymbol{X}_{n}})^{\top}\boldsymbol{J}$

if the diagonal elements of
$\boldsymbol{D}(\boldsymbol{X}_{n})=\boldsymbol{D}(\boldsymbol{X}_{n}^{*})$
are all distinct, where $\boldsymbol{A}$ and $\boldsymbol{b}$ are as above and
$\boldsymbol{J}$ is a $d\times d$ sign-change matrix, i.e., a diagonal matrix
with $\pm 1$ on its diagonal. Thus, if ICS is performed on the same data but
measured in different coordinates systems, the resulting ICs will differ at
most by their signs. Note that in the following, we will drop the dependency
of the matrices on $\boldsymbol{X}_{n}$ in the notation when the context is
clear.

Not only do the values $\lambda_{1},\ldots,\lambda_{d}$ fix the order of the
components, they also have a concrete interpretation as a measure of kurtosis
for the corresponding components in the sense of $\boldsymbol{V}_{1}$ and
$\boldsymbol{V}_{2}$, since they correspond to the ratio of two scale
measures, i.e.,

$\kappa(\boldsymbol{w}^{\top}\boldsymbol{X})=\frac{\boldsymbol{w}^{\top}\boldsymbol{V}_{2}(\boldsymbol{X})\boldsymbol{w}}{\boldsymbol{w}^{\top}\boldsymbol{V}_{1}(\boldsymbol{X})\boldsymbol{w}},$

where $\boldsymbol{w}$ is a $d$-dimensional vector. If
$\boldsymbol{V}_{1}=\mathrm{COV}$ and
$\boldsymbol{V}_{2}=\mathop{\mathrm{COV}_{4}}\nolimits$, this measure relates
to the classical kurtosis measure of Pearson [Nordhausen et al., 2011]. Tyler
et al. [2009] show then that the maximal and minimal kurtosis in the sense of
$\boldsymbol{V}_{1}$ and $\boldsymbol{V}_{2}$ that can be obtained for
$\boldsymbol{X}_{n}\boldsymbol{w}^{\top}$ is captured by $\lambda_{1}$ and
$\lambda_{d}$, respectively.

###  3.2 ICS as preprocessing step for clustering

While ICS has many applications, e.g., in a transformation-retransformation
framework, for outlier detection, or as a method for independent component
analysis [Tyler et al., 2009, Nordhausen et al., 2011, Archimbaud et al.,
2018, Nordhausen et al., 2008], we study ICS in the framework of clustering.

The most common framework of clustering are location mixture models where, if
the class labels would be known, a natural approach would be to maximize the
ratio of the between-group and within-group covariance matrices to obtain
Fisher’s linear discriminant subspace. However, even without knowing the class
labels, the aim of ICS is that the two scatter matrices, which are in general
not proportional with each other in this model, capture this ratio. Indeed,
Tyler et al. [2009] proved that in an elliptical mixture model where the $q$
populations have proportional scatter matrices and different means, ICS in
general will have at most $q-1$ values in $\lambda_{1},\ldots,\lambda_{d}$
that are different from the remaining non-informative values, which are all
equal. Denote the value of the non-informative values by $\lambda_{0}$. For
ICS as a preprocessing step for clustering, the aim is therefore to use only
the invariant coordinates which span Fisher’s linear discriminant subspace,
and to drop the other coordinates. Thus, we want the $l\times d$ matrix
$\boldsymbol{W}_{DS}$ that contains the rows of $\boldsymbol{W}$ corresponding
to interesting values of $\lambda_{1}\ldots,\lambda_{d}$, i.e, those values
with $\lambda_{i}\neq\lambda_{0}$, $i\in\\{1,\ldots,d\\}$. The problem is that
the values of $\lambda_{0}$ and $l\leq q-1$ depend on the distributions of the
groups and their mixture weights, as well as on which scatters
$\boldsymbol{V}_{1}$ and $\boldsymbol{V}_{2}$ are used. Also, unlike PCA where
usually the first few PCs are of interest, the interesting ICs can be the
first or the last ones, or combinations thereof. Assuming that
$\boldsymbol{V}_{1}$ is more robust than $\boldsymbol{V}_{2}$, then
heuristically speaking, based on the generalized kurtosis interpretation, it
is to be expected that groups with quite different group sizes will have large
$\boldsymbol{V}_{1}$-$\boldsymbol{V}_{2}$ kurtosis values while groups with
similar group sizes will have small $\boldsymbol{V}_{1}$-$\boldsymbol{V}_{2}$
kurtosis values. This means large and small clusters will be found in the
different extremes [e.g., Peña et al., 2010].

It is meanwhile well established that there is no general best scatter
combination for ICS, but that some combinations are better suited for certain
applications than others. For example, Peña et al. [2010] consider the
combination $\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$, which is known
as FOBI [Cardoso, 1989, Nordhausen & Virta, 2019], to reveal clusters.
However, this combination is highly sensitive to extreme observations and
therefore, for example, Archimbaud et al. [2018] recommend this combination to
detect outliers. This combination is also of interest due to being completely
moment-based and therefore easy to study. For example, one can show that in a
Gaussian mixture framework, the value of the non-interesting (Gaussian)
components is $\lambda_{0}=1$. For a mixture of two Gaussians with the same
covariance, it can be shown that if the mixture weight for the smaller cluster
is $(3+\sqrt{3})^{-1}\approx 21$%, then the kurtosis of the mixture is the
same as for a pure Gaussian component, independent of the difference in
locations, and therefore also the generalized eigenvalue of such a component
equals 1 [see, e.g., Tyler et al., 2009].

The general motivation is that when using ICS as a preprocessing step for
clustering, one scatter matrix should measure the within-group scatter and the
other the between-group scatter (or equivalently, the total scatter).
Therefore, it seems natural to use one scatter matrix that aims to estimate
the local structure and another that estimates the global structure. For
instance, Hennig [2009] suggests $\mathrm{LCOV}$ as the former. But also
robust scatter matrices could be an interesting option to capture the
dispersion of a single cluster. Specifically, for $\mathrm{MCD}_{\alpha}$, it
might be worth to consider also values of $\alpha$ that are smaller than 0.5
for this purpose.

Note that while different scatter combinations have been considered for ICS in
a clustering framework, a thorough comparison is still missing from the
literature. To conclude this part, we emphasize that ICS is performed in an
unsupervised manner without using class label information. If that information
is available, one scatter matrix involved in ICS can use this information so
that it becomes a supervised scatter matrix, which leads to supervised
invariant coordinate selection (SICS). Many supervised dimension reduction
methods, such as linear discriminant analysis (LDA), sliced inverse regression
(SIR), and sliced average variance estimation (SAVE) can be seen as special
cases of SICS [Liski et al., 2014].

###  3.3 Selection of the invariant coordinates to retain

Another topic often ignored in the literature is the selection of the ICs that
should be retained for clustering. If the number of groups $q$ is known, we
expect from the theoretical properties of ICS, that the number of invariant
coordinates to retain is less than or equal to $q-1$. In what follows, we
compare several criteria for invariant coordinates selection. Two of the
criteria that we consider further are based on the assumption that $q$ is
known, while one criterion makes no such assumption. In a specific data
analysis, the most natural approach might be to look at the scatterplot matrix
of the invariant coordinates and at the generalized eigenvalues to decide
which $l<d$ invariant coordinates carry relevant information for clustering.
However, one may be interested in a more formal approach, therefore we discuss
possible rules in the following.

Assume that the data come from a Gaussian mixture model having equal
covariance matrices with $q<d$ components. In that case, non-interesting
components will have a Gaussian distribution and the informative ones will be
non-Gaussian. This setting also fits into the framework of non-Gaussian
component analysis (NGCA), where the goal is to find a non-Gaussian subspace.
Nordhausen et al. [2017, 2022] show that $\mathrm{COV}-\mathrm{COV}4$ can be
seen as an NGCA method that recovers this subspace, and they derive asymptotic
and bootstrap tests for the number of non-Gaussian components, which can be
used in a successive testing strategy to obtain an estimate for $l$. The main
idea here is that all Gaussian components have a generalized eigenvalue of 1,
hence the variance among these generalized eigenvalues should be small.
Similarly, resampling-based estimates of $l$ for this scatter combination are
considered in Luo & Li [2016, 2021]. These tests with a joint estimation
strategy were extended to a resampling framework for all scatter combinations
in an NGCA setting in Radojičić & Nordhausen [2019]. The difficulty there is
that the value of the non-interesting (Gaussian) generalized eigenvalues is
not known, only that they should all have the same generalized eigenvalue. The
idea in Radojičić & Nordhausen [2019] is to go for a specific $l_{0}$ through
all combinations of $l_{0}$ neighboring generalized eigenvalues to choose the
ones with minimal variance, and to use bootstrapping to decide if assuming
equality is reasonable. However, bootstrapping is quite computationally
expensive in this context. If the number of clusters $q$ can be assumed known,
one rule could be to find the $d-q+1$ generalized eigenvalues with the
smallest variance and keep the components belonging to the non-included
generalized eigenvalues. We refer to this rule as the _var criterion_. While
this approach is outlined above in an NGCA framework, it can also be applied
in a more general framework of elliptical mixtures, where again the
generalized eigenvalues of the purely elliptical part will be equal, and the
generalized eigenvalues of components that contain cluster information will
differ.

A perhaps simpler criterion than the _var criterion_ in an elliptical mixture
model framework is to assume $q\leq d/2$, which implies that the majority of
the generalized eigenvalues should be equal. Then a rule could be to choose
the $q-1$ components whose generalized eigenvalues deviate the most from the
median of all generalized eigenvalues. We refer to this approach as the _med
criterion_. In a Gaussian mixture model, a very simple approach is to apply
marginal tests for normality and to choose all non-normal components, which we
call accordingly the _normal criterion_. This approach was successfully used
in an ICS-outlier detection framework in Archimbaud et al. [2018]. Note that
since those authors work in an outlier detection framework, they look for non-
Gaussian components only among the first components. In a clustering context,
it is natural to take the first and last components into consideration.

To summarize the section about ICS we would like to stress that the key
feature for ICS as a prepocessing step for clustering is that the components
are related to generalized kurtosis measures. Extreme generalized kurtosis
values can be seen as indicators of possible interesting structures for
clustering. This motivation does not strictly depend on the assumption of an
underlying elliptical location-mixture model but can also be argued for other
models with group structures. Nevertheless, the theory and many of the
component selection methods discussed above rely on a Gaussian mixture model,
which is one of the most popular clustering models. Furthermore, ICS is the
most meaningful if the clustering information is contained in lower
dimensional subspace, and the resulting dimension reduction is independent of
the clustering method that is subsequently applied. In the following
simulation study, we apply some popular clustering methods after component
selection using ICS, but in general any other clustering method could be used
instead. This makes ICS a general tool for dimension reduction in a clustering
context. Therefore, this is quite different from the many subspace estimation
methods developed and tailored for specific clustering methods, see for
instance Soete & Carroll [1994], Vichi & Kiers [2001], Bouveyron et al.
[2019].

##  4 Simulation study

In order to investigate the performance of tandem clustering with ICS and PCA,
we perform simulations. The simulation design is described in Section 4.1, and
the main results are presented in Section 4.2. Additional findings are
discussed in Section 4.3, while Section 4.4 provides a brief summary and
discusses some limitations of the simulation study.

###  4.1 Simulation design

Focusing on a data generating process in which the clustering structure lies
in a low-dimensional subspace, we compare the performance of ICS and PCA with
respect to reducing the dimensionality of the data while keeping the relevant
structure for clustering. We study the impact of the scatter pair and the
component selection criterion in different cluster settings and in the
presence of outliers. For assessing the performance of the methods, we compute
$\eta^{2}=1-\Lambda$, where $\Lambda$ denotes Wilks’ lambda. That is, we
compute

$\eta^{2}=1-\frac{\det(\boldsymbol{E})}{\det(\boldsymbol{T})},$

where $\boldsymbol{E}$ is the within-group sum of squares and cross-products
matrix and $\boldsymbol{T}$ is the total sum of squares and cross-products
matrix. The $\eta^{2}$ is a measure of discriminatory power commonly used in
discriminant analysis (see, e.g., McLachlan, 1992; as well as Todorov, 2007,
for a robust version of Wilks’ lambda). It is in the interval $[0,1]$, with a
higher value indicating better discriminatory power.

After dimension reduction, we apply clustering methods to the obtained
components in order to confirm the results on the discriminatory power. Our
aim is not to compare clustering methods, hence we limit our study to five
well-known methods: partitioning around medoids (PAM) [Kaufman & Rousseeuw,
1990], $k$means [Hartigan & Wong, 1979], t$k$means (a trimmed version of
$k$means that identifies a specified proportion of observations as outliers)
[Cuesta-Albertos et al., 1997], as well as model-based clustering with
Gaussian mixtures [Fraley & Raftery, 2002] without allowing for noise (mclust)
and with allowing for noise (rmclust). These clustering methods require
specifying the number of clusters $k$ as an input, which we set equal to the
true number of clusters $q$ to keep the simulations simple. As evaluation
measure, we compute the adjusted Rand index (ARI) [Hubert & Arabie, 1985],
which has a maximum value of 1 in the case of perfect clustering and a value
of 0 if the clustering is similar to a random assignment. For comparison, we
also apply the clustering methods to the simulated data without dimension
reduction, in which case we standardize the data using the mean and standard
deviation before $k$means, and we employ robust standardization using the
median and median absolute deviation (MAD) before PAM and t$k$means.

We generate $n=1000$ observations on $d=10$ variables from the model in (2)
for different numbers of clusters $q\in\\{2,3,5\\}$ with
$\boldsymbol{\Sigma}_{W}=\boldsymbol{I}_{d}$,
$\boldsymbol{\mu}_{1}=\mathbf{0}$, and
$\boldsymbol{\mu}_{h+1}=\delta\boldsymbol{e}_{h}$, $h=1,\dots,q-1$, where
$\delta=10$ and $\boldsymbol{e}_{h}$ is a $d$-dimensional vector with one in
the $h$-th coordinate and zero elsewhere. In this setting, the clustering
structure lies in a low-dimensional subspace of dimension $l=q-1$. We thereby
consider 22 different combinations of the mixture weights
$\epsilon_{1},\dots,\epsilon_{q}$ (see A for details).

In addition to a baseline setting without outliers, we consider settings in
which 2% and 5% of the observations are replaced by outliers. We draw the
outliers from a uniform distribution on a large hyperrectangle that has a
central cutout in the shape of a smaller hyperrectangle around the original
data. The length of each side of the smaller hyperrectangle is given by the
range of the original data in the respective variable, while the corresponding
side of the larger hyperrectangle is twice as long. Note that we treat
outliers as an additional group in the computation of the $\eta^{2}$ and the
ARI.

Regarding scatter matrices, we apply ICS with the following scatter pairs:
$\mathrm{LCOV}-\mathrm{COV}$, $\mathrm{TCOV}-\mathrm{COV}$,
$\mathrm{TCOV}-\mathrm{UCOV}$, $\mathrm{MCD}_{\alpha}-\mathrm{COV}$ and
$\mathrm{RMCD}_{\alpha}-\mathrm{COV}$ with
$\alpha\in\\{0.1,0.2,0.25,0.5,0.75\\}$,
$\mathrm{MCD}_{0.25}-\mathrm{MCD}_{0.95}$ and
$\mathrm{RMCD}_{0.25}-\mathrm{RMCD}_{0.95}$, $\mathrm{MLC}-\mathrm{COV}$, and
$\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$ (see Section 3.1 for details
on those scatter matrices). Note that $\mathrm{MCD}_{\alpha}$ and
$\mathrm{RMCD}_{\alpha}$ are computed using the FAST-MCD algorithm of
[Rousseeuw & Van Driessen, 1999]. Furthermore, we apply PCA using the
correlation matrices derived from $\mathrm{COV}$ and $\mathrm{RMCD}_{0.75}$
[the latter being advocated by, e.g., Croux & Haesbroeck, 1999]. Other robust
methods for PCA such as ROBPCA [Hubert et al., 2005] or MacroPCA [Hubert et
al., 2019] are beyond the scope of our analysis.

Concerning the selection of components, we use four criteria for ICS (see
Section 3.3 for details): the med criterion and the var criterion (which
select $k-1$ components), the normal criterion using the D’Agostino test of
skewness [D’Agostino, 1970] at a 5% significance level (which is not
restricted to selecting a certain number of components), and an oracle
criterion selecting the $k-1$ components (among the first and last components)
that maximize the discriminatory power for the true clusters (as measured by
the $\eta^{2}$). For PCA, we use two criteria: the $80\%$ criterion (which
keeps as many components as necessary to explain at least $80\%$ of the
variability in the data), and the $k-1$ criterion (which keeps the first $k-1$
components).

We simulate 100 data sets for each of the 66 different cluster and outlier
settings. Overall, we compare 72 different dimension reduction strategies: 68
for ICS (17 scatter pairs and 4 criteria) and 4 for PCA (2 scatters and 2
criteria).

###  4.2 Main results

Due to the large number of dimension reduction strategies, we only show a
selection of the most relevant results here. Section 4.2.1 presents the
results for the best performing methods in terms of discriminatory power of
the selected components, while Section 4.2.2 discusses clustering performance
of selected clustering methods after dimension reduction.

A detailed discussion on the selected dimension reduction methods and
component selection criteria is given in B.1. For ICS, the best performing
scatter pairs are $\mathrm{LCOV}-\mathrm{COV}$, $\mathrm{TCOV}-\mathrm{COV}$,
as well as the robust scatter pair $\mathrm{TCOV}-\mathrm{UCOV}$. Furthermore,
we compare the med criterion with the oracle criterion for component
selection. The former is a representative selection since there is not much
difference in the selection criteria for those best performing scatter pairs,
while the latter is included for reference purposes as it uses information on
the true clusters and can therefore not be computed in practice. For PCA, the
best performing method is $\mathrm{RMCD}_{0.75}$ combined with the $80\%$
criterion for selecting the components.

####  4.2.1 Discriminatory power of selected components

We first focus on the overall results across the investigated cluster
settings. Figure 2 shows boxplots of the discriminatory power of the selected
methods for different outlier settings.

Figure 2: Discriminatory power as measured by the $\eta^{2}$ of the best
performing dimension reduction methods for ICS and PCA, respectively, using
the corresponding best component selection criteria. Boxplots display the
results across the different cluster settings with $100$ simulation runs each.
Results for different outlier settings are shown in separate panels.

With and without outliers, ICS with $\mathrm{LCOV}-\mathrm{COV}$ and
$\mathrm{TCOV}-\mathrm{COV}$ clearly outperforms PCA with
$\mathrm{RMCD}_{0.75}$ in terms of discriminatory power. This is not
surprising, as PCA is not designed to pick up the relevant structure for
clustering. Although ICS with the robust scatter pair
$\mathrm{TCOV}-\mathrm{UCOV}$ also clearly outperforms PCA with
$\mathrm{RMCD}_{0.75}$ when there are no outliers or 5% outliers, a conclusion
cannot be drawn in the setting with 2% outliers. While the median $\eta^{2}$
is much higher for $\mathrm{TCOV}-\mathrm{UCOV}$ than for
$\mathrm{RMCD}_{0.75}$, the former also exhibits many points with a very small
$\eta^{2}$ in Figure 2. This behavior is investigated in more detail below in
the context of different cluster settings.

As the amount of outliers increases, we observe that the discriminatory power
decreases for all methods, together with a slight increase in variability.
Nevertheless, the discriminatory power of ICS with
$\mathrm{LCOV}-\mathrm{COV}$ and $\mathrm{TCOV}-\mathrm{COV}$ remains good in
all settings. For PCA, the decrease in discriminatory power is much more
pronounced despite using the robust scatter matrix $\mathrm{RMCD}_{0.75}$.
This indicates that PCA as a preprocessing technique for clustering is much
more affected by outliers than ICS with a suitable scatter pair.

Next, we take a look at specific cluster settings to gain further insights.
Since we investigate a large number of cluster settings, we present boxplots
for a representative selection of cluster settings in Figure 3. The main
finding is that $\mathrm{TCOV}-\mathrm{COV}$ may be preferred over
$\mathrm{LCOV}-\mathrm{COV}$, as the former has lower variability for two
perfectly balanced clusters (50–50). However, both scatter pairs yield very
similar and excellent discriminatory power in the investigated cluster
settings. Furthermore, note that the med criterion essentially yields the same
discriminatory power as the oracle criterion for those two scatter pairs.

Figure 3: Discriminatory power as measured by the $\eta^{2}$ of the best
performing dimension reduction methods for ICS and PCA, respectively, using
the corresponding best component selection criteria. Boxplots display the
results for selected cluster settings over $100$ simulation runs. Results for
different outlier settings are shown in separate rows, and results for
different methods are shown in separate columns.

Figure 3 also provides more insights into the behavior of the robust scatter
pair $\mathrm{TCOV}-\mathrm{UCOV}$. Specifically, we observe that this scatter
pair leads to similar discriminatory power as $\mathrm{LCOV}-\mathrm{COV}$ and
$\mathrm{TCOV}-\mathrm{COV}$, except in the setting with $q=2$ highly
unbalanced clusters (95–5) and with 2% of outliers. In this specific
situation, the discriminatory power considerably drops towards 0 and the
variability increases. This is due to an artifact of our simulation design,
where the outlyingness structure is likely to be constrained to a low-
dimensional subspace if the number of outliers is small. Both $\mathrm{TCOV}$
and $\mathrm{UCOV}$ treat a small cluster as outliers, meaning that this
scatter pair tends to ignore (some of) the relevant structure in such
settings. There is of course a philosophical discussion as to where the line
is between a group of outliers and a small cluster of observations, but such a
discussion is beyond the scope of this paper.

For PCA with $\mathrm{RMCD}_{0.75}$, we find much larger variability compared
to the other cluster settings when we have two clusters and both clusters are
larger than the trimming proportion (in Figure 3, we show 50–50 and 70–30).
Moreover, in all investigated cluster settings, the discriminatory power of
PCA with $\mathrm{RMCD}_{0.75}$ is lower or at best comparable to that of ICS
with $\mathrm{TCOV}-\mathrm{COV}$ and $\mathrm{LCOV}-\mathrm{COV}$.

####  4.2.2 Clustering performance

In order to validate our findings with respect to discriminatory power, we
apply five popular and relatively simple clustering methods. To keep
computation time low, we apply the clustering methods only for a subset of
cluster settings (see A for details).

Figure 4: Clustering performance as measured by the ARI after the best
performing dimension reduction methods for ICS and PCA, respectively, using
the corresponding best component selection criteria. Boxplots display the
results across the different cluster settings with $100$ simulation runs each.
Results for different clustering methods are shown in separate columns, and
results for different outlier settings are shown in separate rows.

Figure 4 displays overall results of the ARI of the selected methods. It
contains boxplots across the selected cluster settings for the different
clustering methods and outlier settings. Compared to the ARI on the original
data without dimension reduction, ICS with $\mathrm{LCOV}-\mathrm{COV}$ and
$\mathrm{TCOV}-\mathrm{COV}$ improves the performance of PAM, $k$means,
t$k$means, and mclust. For PAM and $k$means, there is not much change in the
median ARI as it is already near perfect for the observed data, but
variability is reduced and there are fewer instances where clustering fails.
For t$k$means, we observe a similar effect, taking into account that this
method trims away some observations so that the resulting outlier component is
reflected accordingly in a somewhat lower ARI. For mclust, ICS yields a
substantial improvement in ARI when outliers are present. For rmclust, on the
other hand, the ARI is almost perfect for the observed data, as expected in
this simulation design with Gaussian clusters and outliers. Nevertheless,
there are some instances where rmclust fails. Note that after ICS with
$\mathrm{LCOV}-\mathrm{COV}$ or $\mathrm{TCOV}-\mathrm{COV}$ and a simple
clustering method such as $k$means lead to an ARI as high as for rmclust for
any situation and never completely breaks down.

As the amount of outliers increases, the median ARI of PAM and $k$means
decreases while the median ARI of t$k$means increases accordingly. This
behavior is expected and can be observed on the original data and after ICS.
Also note that in the setting with 5% outliers, t$k$means on the observed data
performs almost perfectly in most instances (but fails in some instances).
Keep in mind that for t$k$means, we used the default trimming proportion of
5%. In practice, the trimming proportion needs to be carefully tuned, and the
interplay of the trimming proportion and the number of clusters is quite
complex [García-Escudero et al., 2011]. Hence, this result could be due to the
trimming proportion being equal to the true amount of outliers.

Otherwise, PAM, $k$means, and t$k$means after ICS are not much affected by the
outliers. This reflects that the discriminatory power of the selected
components remains high, even though the aforementioned clustering methods
fail in some instances for $\mathrm{TCOV}-\mathrm{UCOV}$. This phenomenon is
more pronounced for mclust and rmclust, where all scatter pairs yield a
considerable number of unsuccessful clustering instances.

PCA with $\mathrm{RMCD}_{0.75}$ deteriorates the overall performance of all
clustering methods compared to clustering the observed data, either by a quite
substantially lower median or a higher variability in ARI. This reflects that
PCA struggles to capture the relevant structure for clustering.

From the overall results across cluster settings, we can already see that the
results for clustering performance reflect the results for discriminatory
power very well. A detailed discussion of the clustering performance for the
selected cluster settings is therefore moved to B.2.

###  4.3 Additional findings

As we have seen from the previous section that clustering performance nicely
follows the results for the discriminatory power of the selected components,
this section is focused on discriminatory power to keep the paper concise.
Section 4.3.1 presents insights into the performance of
$\mathrm{MCD}_{\alpha}-\mathrm{COV}$ scatter pairs for various values of the
parameter $\alpha$, while Section 4.3.2 discusses results for the popular
$\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$ scatter pair.

####  4.3.1 $\mathrm{MCD}_{\alpha}-\mathrm{COV}$ scatter pairs

The MCD is typically applied in the context of robust statistics, where a
subset size based on the parameter $\alpha\geq 0.5$ is used to capture the
covariance structure of the majority of the data. However, in the context of
ICS as a preprocessing step for clustering, $\alpha\leq 0.5$ is better suited
to capture the local within-cluster structure. Hence, we focus on the scatter
pairs $\mathrm{MCD}_{0.1}-\mathrm{COV}$, $\mathrm{MCD}_{0.25}-\mathrm{COV}$,
and $\mathrm{MCD}_{0.5}-\mathrm{COV}$. Details on the selection of these
scatter pairs can be found in B.1.

Figure 5: Discriminatory power as measured by the $\eta^{2}$ of ICS with
$\mathrm{MCD}_{\alpha}-\mathrm{COV}$ scatter pairs, using the best component
selection criteria. Boxplots display the results for selected cluster settings
over $100$ simulation runs. Results for different outlier settings are shown
in separate rows, and results for different methods are shown in separate
columns.

Figure 5 shows boxplots of the discriminatory power for the different outlier
settings and a representative selection of cluster settings. The main finding
is that the performance of $\mathrm{MCD}_{\alpha}-\mathrm{COV}$ scatter pairs
in different cluster settings depends on the parameter $\alpha$. We first
focus on the setting without outliers. For $q=2$ clusters, for instance, a
lower $\alpha$ performs better for balanced clusters but yields larger
variability for unbalanced clusters (90–10, 95–5), while a larger $\alpha$ is
preferable for unbalanced clusters but yields lower discriminatory power or
larger variability for relatively balanced clusters (50–50, 55–45, 60–40).
However, as the number of clusters increases, this becomes less of an issue.

Regarding the selection of components, the med criterion yields a performance
similar to that of the oracle criterion, except that the aforementioned effect
is more pronounced. Specifically, the med criterion fails or exhibits large
variability for some cases (50-50 and 55-45 for $\alpha=0.25$ and
$\alpha=0.5$, and additionally 60–40 for $\alpha=0.5$). With outliers, the
discriminatory power somewhat drops but remains good overall. Otherwise,
results remain similar to the setting without outliers.

Overall, the results reveal that ICS with $\mathrm{MCD}_{\alpha}-\mathrm{COV}$
has potential in a clustering context, but further research is needed. First,
since results in different cluster settings depend on the parameter $\alpha$,
it could be interesting to consider $\alpha$ as a tuning parameter to be
selected, e.g., together with the number of clusters. Second, the results for
the oracle criterion indicate that the scatter pair
$\mathrm{MCD}_{\alpha}-\mathrm{COV}$ succeeds in revealing the relevant
structure for clustering, but component selection criteria that can be
computed in practice fail in certain cluster settings.

####  4.3.2 $\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$ scatter pair

As one of the most widely used scatter pairs in ICS, we take a closer look at
the performance of $\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$ in a
clustering context [see also Peña et al., 2010]. Boxplots of the
discriminatory power for different cluster and outlier settings are displayed
in Figure 6.

Figure 6: Discriminatory power as measured by the $\eta^{2}$ of ICS with the
$\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$ scatter pair, using the best
component selection criteria. Boxplots display the results for all
investigated cluster settings over $100$ simulation runs. Results for
different outlier settings are shown in separate rows.

We first discuss the setting without outliers. For $q=2$ clusters, the results
are as expected from the theory explained in Section 3.2. In that case, it is
known that the relevant structure is contained in the first component for very
unbalanced clusters and in the last component for (somewhat) balanced
clusters, while neither the first nor the last component can separate the
clusters if the smaller cluster contains about 21% of observations. Indeed, we
find that discriminatory power is very high except for the cluster settings
75–25 and 80–20, where the method fails. For $q=3$ clusters, the method still
cannot recover the cluster with 20% of observations (20–50–30, 10–70–20) but
the discriminatory power is quite good since the two remaining clusters are
highlighted. For $q=5$ clusters, discriminatory power is very high throughout
even though there are various cluster settings with clusters of size 20%. This
is a promising result that needs to be verified from a theoretical point of
view in further research. Furthermore, the med criterion is close to the
oracle criterion in terms of performance, except in settings where the method
struggles either way.

With 2% of outliers, there is a considerable drop in performance for the
oracle criterion. But more importantly, the med criterion fails almost
completely. In particular, the discriminatory power is close to 0 for $q=2$
clusters, while the discriminatory power is very low and variability is large
for $q=3$ clusters. In these settings, the med criterion fails to pick the
components that carry the relevant structure. Surprisingly, the results of the
med criterion improve again for 5% of outliers, being more closely aligned
again with the results of the oracle criterion. In fact, with 2% of outliers
and our sample size, the outlyingness structure is artificially lying in a
low-dimensional space. $\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$ is
designed to discover such a structure and so favors its highlighting on its
first component. However, the clusters are well visible in the last component,
but this one is not kept with the med criteria based on only $k-1$ components.
When the percentage of outliers increases, the outlying observations span the
entire space and so, the low-dimensional clustering structure is again found
as the most important one, explaining the good results with 5%.

###  4.4 Discussion of the simulation study

The main findings from our simulations are that (i) ICS with
$\mathrm{TCOV}-\mathrm{COV}$ and $\mathrm{LCOV}-\mathrm{COV}$ performs best
for tandem clustering, (ii) the scatter pair
$\mathrm{MCD}_{\alpha}-\mathrm{COV}$ with a value of $\alpha\leq 0.5$ could be
an interesting alternative when $\alpha$ is treated as a tuning parameter but
further research is needed, and (iii) PCA is not a suitable preprocessing step
for clustering purposes and yields large variability in clustering
performance.

However, there are several limitations of our simulations. First, we only
consider one value each for the number of observations $n$ and the number of
variables $d$. We compare a large number of scatter pairs for ICS with
different component selection criteria, and our main aim is to study those for
different numbers of clusters with various cluster sizes. As such, we had to
make some concessions to other parameters to keep the paper concise. Second,
we only consider a simple mean shift model for $q$ clusters that are defined
in a subspace of dimension $q-1$. An important consideration for this choice
is that it allows us to use simple clustering methods such as $k$means to
validate our findings for the discriminatory power of the selected components.
More complex models with, e.g., different covariance structures in different
clusters or non-normal distributions are beyond the scope of this paper and
are left for further research.

As a short glimpse into the performance of tandem ICS in a non-normal setting,
we provide an illustrative example in B.3 with a simulated data set coming
from the so-called barrow wheel distribution. This distribution goes beyond
Gaussian mixtures, and it is mentioned in the discussion of Tyler et al.
[2009], where the authors noticed that the theoretical properties of ICS are
also valid for this distribution [see also Archimbaud, 2018, for more
theoretical details].

##  5 Empirical applications

To further study the performance of tandem clustering with ICS, we consider
three publicly available benchmark datasets. First, the _crabs_ data set
[Campbell & Mahon, 1974] contains morphological measurements on crabs of both
sexes and from two color-different species. We take log-transformations of all
variables, which is common for this data set [e.g. Archimbaud et al., 2023].
Second, the well-known _iris_ data set [Anderson, 1936, Fisher, 1936] gathers
measurements of flowers from three different species. Third, the _Philips_
data set contains measurements of a certain part for television sets,
collected in 1997 by Philips Mecoma (The Netherlands) and studied by Rousseeuw
& Van Driessen [1999] and Hubert et al. [2008] for anomaly detection purposes.
In this specific case, the engineers agreed on two clear abnormal phenomena
regarding 15% and 11% of the data, which we treat as small clusters. An
overview of the dimensions and cluster sizes of the three data sets is given
in Table 1.

Name | $n$ | $d$ | $q$ | Cluster sizes (%) |  Source (R package)
---|---|---|---|---|---
crabs | 200 | 5 | 4 | 25–25–25–25 |  MASS [Venables & Ripley, 2002]
iris | 150 | 4 | 3 | 33–33–33 |  datasets [R Core Team, 2021]
Philips | 677 | 9 | 3 | 74–15–11 |  cellWise [Raymaekers & Rousseeuw, 2022]
Table 1: Overview of the three benchmark data sets, where $n$ denotes the
number of observations, $d$ the number of variables, and $q$ the number of
true clusters.

We focus on the clustering performance on these data sets using the same five
clustering methods and the best performing dimension reduction methods from
Section 4, thereby selecting the components using the med, var, and normal
criteria for ICS, as well as the $80\%$ criterion for PCA. The clustering
performance is again evaluated by the Adjusted Rand Index (ARI) [Hubert &
Arabie, 1985]. The results are shown in Figure 7, where reference lines
indicate the ARI on the initial data without dimension reduction. Even though
the data sets are relatively low-dimensional, ICS improves the identification
of clusters in most cases, as reflected by an increase in the ARI. As in the
simulation results, the med and var criteria lead to similar performance. By
contrast, PCA yields roughly the same or a lower ARI in all cases except for
mclust and rmclust for the Philips data. Furthermore, the Philips data set is
an interesting case where mclust and rmclust yield a lower ARI on the observed
data than simple methods such as PAM, $k$means, and t$k$means. Note that the
ARI slightly deteriorates in some instances also for ICS, hence the selection
of the scatter pair and the component selection criterion is crucial. For a
better understanding, we take a deeper look at the components obtained via ICS
in order to relate them to the clustering performance. We thereby focus on the
performance of $k$means and mclust, as the results are qualitatively similar
for $k$means, t$k$means and PAM, as well as for mclust and rmclust.

Figure 7: Adjusted Rand Index (ARI) in the empirical applications for
different clustering methods after selected dimension reduction methods for
ICS and PCA, respectively, using different component selection criteria.
Results for different clustering methods are shown in separate rows, and
results for different data sets are shown in separate columns. Black reference
lines indicate the ARI obtained on initial data without dimension reduction.
Figure 8: First two and last two components from ICS with
$\mathrm{TCOV}-\mathrm{COV}$ on the crabs data set: scatterplot matrix with
density estimates on the diagonal. The true clusters are indicated by
different colors.

For the crabs data set, the ARI on the initial data is is 0.04 for $k$means
and 0.37 for mclust. The improvement is substantial for ICS with the scatter
pairs $\mathrm{TCOV}-\mathrm{COV}$ and $\mathrm{TCOV}-\mathrm{UCOV}$, each
yielding an ARI between 0.79 and 0.88 depending on the criterion used and the
clustering method. The $\mathrm{MCD}_{\alpha}-\mathrm{COV}$ scatter pairs, on
the other hand, achieve almost the same performance when $\alpha$ is small but
the ARI decreases when $\alpha$ increases. Note that
$\mathrm{LCOV}-\mathrm{COV}$ works well for $k$means but not for mclust.
Another interesting fact is that the normal criterion fails to select any
components for any of the scatter pairs. To investigate this further, Figure 8
displays the first two and last two components from ICS with the scatter pair
$\mathrm{TCOV}-\mathrm{COV}$, which clearly shows that the first two
components reveal the relevant structure for clustering. The plot also
visualizes the nested nature of the $q=4$ clusters in the crabs data: there
are two species (blue and orange) and two sexes (male and female). ICS
highlights the two balanced groups regarding the species on one component
(here, the first) and the two balanced clusters regarding the sex on the
another component (here, the second). However, for two balanced groups, the
D’Agostino test of skewness [D’Agostino, 1970] does not reject the null
hypothesis of normality, explaining why the normal criterion results in the
selection of no components. The med criterion, on the other hand, is based on
the selection of $k-1=3$ components (where $k=4$ denotes the number of
clusters specified in the clustering algorithm). It selects the first two
components, which carry the relevant structure, but also the last component.
It is possible that clustering can be further improved by selecting only the
first two components.

For the iris data set, the ARI on the initial data is already excellent for
mclust with a value of 0.90. However, it is noteworthy that for $k$means, the
ARI of 0.62 on the initial data with be improved with the scatter pairs
$\mathrm{LCOV}-\mathrm{COV}$, $\mathrm{TCOV}-\mathrm{COV}$,
$\mathrm{TCOV}-\mathrm{UCOV}$, and $\mathrm{MCD}_{0.25}-\mathrm{COV}$ to
achieve competitive ARI values between 0.87 and 0.92. For the scatter pairs
involving $\mathrm{TCOV}$, the normal criterion selects only the first
component, whereas the med criterion selects the first and the last
components,and the var criterion keeps the first two components, with the
former resulting in a much higher ARI. Figure 9 displays the four components
from ICS with the scatter pair $\mathrm{TCOV}-\mathrm{COV}$, confirming that
the relevant structure is indeed contained only in the first component. For
$\mathrm{LCOV}-\mathrm{COV}$, the med criterion also selects the first and
last component, here resulting in an excellent ARI, while the normal criterion
does not select any component. It should be noted that the iris data set was
used as an example in Tyler et al. [2009], too. Those authors also find that
the relevant structure is captured by the first component, but they conjecture
that the results are similar “for almost any pair of scatter matrices that we
may choose”. Our results are an indication that the choice of scatter pair
matters more for this data set than initially thought.

Figure 9: All four components from ICS with $\mathrm{TCOV}-\mathrm{COV}$ on
the iris data set: scatterplot matrix with density estimates on the diagonal.
The true clusters are indicated by different colors.

For the Philips data set, the ARI on the initial data is rather low with a
value of 0.26 for $k$means, but even lower for mclust (0.06). The most
striking result in this example is that the ARI is by far the highest after
ICS with $\mathrm{MCD}_{0.5}-\mathrm{COV}$ with a value of 0.90 for $k$means,
while the other scatter pairs yield only a moderate improvement in the ARI
(which is also the case for mclust). Furthermore, the med, var, and normal
criteria yield similar results in terms of ARI for all considered scatter
pairs. This is surprising as the med criterion always selects the first two
components, whereas the normal criterion selects many more (the first five in
the case of $\mathrm{MCD}_{0.5}-\mathrm{COV}$). Figure 10 indeed confirms that
the first two components are necessary to reveal the clusters for the
$\mathrm{MCD}_{0.5}-\mathrm{COV}$ scatter pair. For this specific data set, it
is important to recall that the clusters have been labeled a-posteriori by
engineers at Philips based on Mahalanobis distances computed via the
$\mathrm{MCD}$. This might explain why other scatter pairs do not succeed in
the same way as $\mathrm{MCD}_{0.5}-\mathrm{COV}$ to reveal the structure.

Figure 10: First two and last two components from ICS with
$\mathrm{MCD}_{0.50}-\mathrm{COV}$ on the Philips data set: scatterplot matrix
with density estimates on the diagonal. The true clusters are indicated by
different colors.

Note that further results for the three considered data sets can be found in
C, where we include two methods that integrate dimension reduction into the
$k$means algorithm. On these data sets, tandem clustering with ICS performs
better than those two approaches.

##  6 Conclusions and discussion

Overall, ICS outperforms PCA for tandem clustering. More specifically, the
scatter pairs $\mathrm{TCOV}-\mathrm{COV}$ and $\mathrm{LCOV}-\mathrm{COV}$
yield the best performance, but also $\mathrm{MCD}_{\alpha}-\mathrm{COV}$ with
a carefully selected parameter $\alpha\leq 0.5$ is promising. However, there
is no scatter pair that is uniformly best in every situation. Regarding the
component selection criteria, note that we restrict the med criterion to
selecting $k-1$ components, where $k$ denotes the number of clusters to be
specified in the clustering algorithm. The main reason for this restriction is
to have a simple, automatic component selection rule for the simulations. In
practice, an even smaller subspace may carry the relevant structure, hence the
med criterion is not well adapted in this context and may select too many
components. On the other hand, the normal criterion can estimate the dimension
of the subspace of interest, but the D’Agostino test [D’Agostino, 1970]
struggles when there are two balanced groups on one component, which can
result in selecting none of the components, or in presence of outliers. Other
normality tests may be more suitable. Further research on suitable component
selection criteria is necessary, but in practice, researchers and
practitioners may also select the components based on a scatterplot matrix of
the components or a screeplot of the generalized eigenvalues. A more advanced
procedure could be to consider the choice of the invariant coordinates as a
tuning parameter of the subsequent clustering method. In this context, the
best subspace could be chosen based on an internal validation index such as
the silhouette coefficient [Kaufman & Rousseeuw, 1990] or those described in
Hennig [2019].

In addition to improving the component selection criteria, the present work
has many perspectives, including the treatment of mixed data [Caussinus &
Ruiz-Gazen, 2006, van de Velden et al., 2019b] and functional data [Archimbaud
et al., 2022], but also the improvement of the ICS calculation algorithm by
allowing to take into account (nearly) ill-conditioned data [Archimbaud et
al., 2023], or by refining the research of the components using a projection
pursuit algorithm as suggested in Dümbgen et al. [2023]. Generalizing ICS to
high-dimensional data through some regularization approach is another
interesting perspective in the context of clustering [see Bouveyron & Brunet-
Saumard, 2014, for a review for model-based clustering]. In the context of
rating-scale data from surveys in the social and behavioral sciences, tandem
clustering with ICS may further be useful for detecting different response
styles among respondents [e.g. Schoonees et al., 2015].

## Computational details

All computations are performed with R version 4.2.2 [R Core Team, 2021].
Replication files for the simulations and the empirical applications will be
made available online upon acceptance of this manuscript, and a link will be
included here.

## Acknowledgements

This work was partly supported by a grant of the _Dutch Research Council_
(NWO, research program Vidi, project number VI.Vidi.195.141), by the Austrian
Science Fund P31881-N32, by the HiTEc COST Action (CA21163), and by the French
_Agence Nationale de la Recherche_ under grant ANR-17-EURE-0010
(Investissements d’Avenir program).

## References

  * Aggarwal & Reddy [2014] Aggarwal, C. C., & Reddy, C. K. (Eds.) (2014).  Data Clustering: Algorithms and Applications.  CRC Press.  URL: http://www.charuaggarwal.net/clusterbook.pdf.
  * Alashwali & Kent [2016] Alashwali, F., & Kent, J. T. (2016).  The use of a common location measure in the invariant coordinate selection and projection pursuit.  Journal of Multivariate Analysis, 152, 145–161. doi:10.1016/j.jmva.2016.08.007.
  * Anderson [1936] Anderson, E. (1936).  The species problem in iris.  Annals of the Missouri Botanical Garden, 23, 457–509. doi:10.2307/2394164.
  * Arabie & Hubert [1994] Arabie, P., & Hubert, L. (1994).  Cluster analysis in marketing research.  In R. Bagozzi (Ed.), Advanced methods in marketing research (pp. 160–189).  Oxford: Blackwell.
  * Archimbaud [2018] Archimbaud, A. (2018).  Statistical methods for outlier detection for high-dimensional data.  Ph.D. thesis Thèse de doctorat, Université Toulouse 1 Capitole.
  * Archimbaud et al. [2022] Archimbaud, A., Boulfani, F., Gendre, X., Nordhausen, K., Ruiz-Gazen, A., & Virta, J. (2022).  ICS for multivariate functional anomaly detection with applications to predictive maintenance and quality control.  Econometrics and Statistics, in press. doi:10.1016/j.ecosta.2022.03.003.
  * Archimbaud et al. [2023] Archimbaud, A., Drmač, Z., Nordhausen, K., Radojičić, U., & Ruiz-Gazen, A. (2023).  Numerical considerations and a new implementation for ICS.  SIAM Journal on Mathematics of Data Science (SIMODS), 5, 97–121. doi:10.1137/22M1498759.
  * Archimbaud et al. [2018] Archimbaud, A., Nordhausen, K., & Ruiz-Gazen, A. (2018).  ICS for multivariate outlier detection with application to quality control.  Computational Statistics & Data Analysis, 128, 184–199. doi:10.1016/j.csda.2018.06.011.
  * Art et al. [1982] Art, D., Gnanadesikan, R., & Kettenring, J. (1982).  Data-based metrics for cluster analysis.  Utilitas Mathematica A, 21, 75–99.
  * Bouveyron & Brunet-Saumard [2014] Bouveyron, C., & Brunet-Saumard, C. (2014).  Model-based clustering of high-dimensional data: A review.  Computational Statistics & Data Analysis, 71, 52–78.
  * Bouveyron et al. [2019] Bouveyron, C., Celeux, G., Murphy, T. B., & Raftery, A. E. (2019).  Model-Based Clustering and Classification for Data Science: With Applications in R.  Cambridge: Cambridge University Press.  doi:10.1017/9781108644181.
  * Campbell & Mahon [1974] Campbell, N., & Mahon, R. (1974).  A multivariate study of variation in two species of rock crab of the genus leptograpsus.  Australian Journal of Zoology, 22, 417–425. doi:10.1071/zo9740417.
  * Cardoso [1989] Cardoso, J.-F. (1989).  Source separation using higher order moments.  In Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (pp. 2109–2112).  IEEE.  doi:10.1109/ICASSP.1989.266878.
  * Cator & Lopuhaä [2012] Cator, E. A., & Lopuhaä, H. P. (2012).  Central limit theorem and influence function for the mcd estimators at general multivariate distributions.  Bernoulli, 18, 520 – 551. doi:10.3150/11-BEJ353.
  * Caussinus et al. [2003] Caussinus, H., Fekri, M., Hakam, S., & Ruiz-Gazen, A. (2003).  A monitoring display of multivariate outliers.  Computational Statistics & Data Analysis, 44, 237–252. doi:10.1016/S0167-9473(03)00059-8.
  * Caussinus & Ruiz [1990] Caussinus, H., & Ruiz, A. (1990).  Interesting projections of multidimensional data by means of generalized principal component analyses.  In K. Momirović, & V. Mildner (Eds.), Compstat (pp. 121–126).  Heidelberg: Physica-Verlag HD.  doi:10.1007/978-3-642-50096-1_19.
  * Caussinus & Ruiz-Gazen [1993] Caussinus, H., & Ruiz-Gazen, A. (1993).  Projection pursuit and generalized principal component analysis.  In S. Morgenthaler, E. Ronchetti, & W. A. Stahel (Eds.), New directions in statistical data analysis and robustness. Monte Verita, Proceedings of the Centro Stefano Franciscini Ascona Series.  Springer.
  * Caussinus & Ruiz-Gazen [1995] Caussinus, H., & Ruiz-Gazen, A. (1995).  Metrics for finding typical structures by means of principal component analysis.  In Data Sci. and its Applications (pp. 177–192).  Japan: Harcourt Brace.
  * Caussinus & Ruiz-Gazen [2006] Caussinus, H., & Ruiz-Gazen, A. (2006).  Projection pursuit approach for categorical data.  Multiple Correspondence Analysis and Related Methods. Chapman and Hall/CRC, London, (pp. 405–418).
  * Caussinus & Ruiz-Gazen [2007] Caussinus, H., & Ruiz-Gazen, A. (2007).  Classification and generalized principal component analysis.  In Selected contributions in data analysis and classification (pp. 539–548).  Springer.  doi:10.1007/978-3-540-73560-1_50.
  * Cerioli [2005] Cerioli, A. (2005).  K-means cluster analysis and mahalanobis metrics: a problematic match or an overlooked opportunity.  Statistica Applicata, 17, 61–73.
  * Chang [1983] Chang, W.-C. (1983).  On using principal components before separating a mixture of two multivariate normal distributions.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 32, 267–275. doi:10.2307/2347949.
  * R Core Team [2021] R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing Vienna, Austria.  URL: https://www.R-project.org/.
  * Croux & Haesbroeck [1999] Croux, C., & Haesbroeck, G. (1999).  Influence function and efficiency of the minimum covariance determinant scatter matrix estimator.  Journal of Multivariate Analysis, 71, 161–190. doi:10.1006/jmva.1999.1839.
  * Cuesta-Albertos et al. [1997] Cuesta-Albertos, J. A., Gordaliza, A., & Matrán, C. (1997).  Trimmed $k$-means: An attempt to robustify quantizers.  The Annals of Statistics, 25, 553–576. doi:10.1214/aos/1031833664.
  * D’Agostino [1970] D’Agostino, R. B. (1970).  Transformation to normality of the null distribution of g1.  Biometrika, (pp. 679–681). doi:10.2307/2334794.
  * Ding & Li [2007] Ding, C., & Li, T. (2007).  Adaptive dimension reduction using discriminant analysis and k-means clustering.  In Proceedings of the 24th international conference on Machine learning (pp. 521–528).  doi:10.1145/1273496.1273562.
  * Dümbgen et al. [2023] Dümbgen, L., Gysel, K., & Perler, F. (2023).  Refining invariant coordinate selection via local projection pursuit.  In M. Yi, & K. Nordhausen (Eds.), Robust and Multivariate Statistical Methods: Festschrift in Honor of David E. Tyler (pp. 121–136).  Cham: Springer.  doi:10.1007/978-3-031-22687-8_6.
  * Fekri & Ruiz-Gazen [2015] Fekri, M., & Ruiz-Gazen, A. (2015).  A B-robust non-iterative scatter matrix estimator: Asymptotics and application to cluster detection using invariant coordinate selection.  In K. Nordhausen, & S. Taskinen (Eds.), Modern Nonparametric, Robust and Multivariate Methods: Festschrift in Honour of Hannu Oja (pp. 395–423).  Cham: Springer International Publishing.  doi:10.1007/978-3-319-22404-6_22.
  * Fischer et al. [2017] Fischer, D., Honkatukia, M., Tuiskula-Haavisto, M., Nordhausen, K., Cavero, D., Preisinger, R., & Vilkki, J. (2017).  Subgroup detection in genotype data using invariant coordinate selection.  BMC Bioinformatics, 18, 173–181. doi:10.1186/s12859-017-1589-9.
  * Fischer et al. [2020] Fischer, D., Nordhausen, K., & Oja, H. (2020).  On linear dimension reduction based on diagonalization of scatter matrices for bioinformatics downstream analyses.  Heliyon, 6, e05732. doi:10.1016/j.heliyon.2020.e05732.
  * Fisher [1936] Fisher, R. A. (1936).  The use of multiple measurements in taxonomic problems.  Annals of Eugenics, 7, 179–188. doi:10.1111/j.1469-1809.1936.tb02137.x.
  * Fraley & Raftery [2002] Fraley, C., & Raftery, A. E. (2002).  Model-based clustering, discriminant analysis, and density estimation.  Journal of the American statistical Association, 97, 611–631.
  * García-Escudero et al. [2011] García-Escudero, L. A., Gordaliza, A., Matrán, C., & Mayo-Scar, A. (2011).  Exploring the number of groups in robust model-based clustering.  Statistics and Computing, 21, 585–599. doi:10.1007/s11222-010-9194-z.
  * Gnanadesikan et al. [1993] Gnanadesikan, R., Harvey, J., & Kettenring, J. (1993).  Mahalanobis metrics for cluster analysis.  Sankhyā: The Indian Journal of Statistics, Series A, (pp. 494–505). URL: https://www.jstor.org/stable/25050956.
  * Hampel et al. [1986] Hampel, F., Ronchetti, E., Rousseeuw, P., & Stahel, W. (1986).  Robust statistics john wiley & sons.  New York, .
  * Hartigan & Wong [1979] Hartigan, J. A., & Wong, M. A. (1979).  Algorithm as 136: A k-means clustering algorithm.  Journal of the royal statistical society. series c (applied statistics), 28, 100–108. doi:10.2307/2346830.
  * Hennig [2009] Hennig, C. (2009).  Discussion of “Invariant Co-ordinate Selection”, by D. E. Tyler, F. Critchley, L. Dümbgen, and H. Oja.  Journal of the Royal Statistical Society B, 71, 579–583. doi:10.1111/j.1467-9868.2009.00706.x.
  * Hennig [2019] Hennig, C. (2019).  Cluster validation by measurement of clustering characteristics relevant to the user.  In Data Analysis and Applications 1 chapter 1. (pp. 1–24).  John Wiley & Sons, Ltd.  doi:10.1002/9781119597568.ch1.
  * Hennig et al. [2015] Hennig, C., Meila, M., Murtagh, F., & Rocci, R. (2015).  Handbook of cluster analysis.  CRC Press.  doi:10.1201/b19706.
  * Hubert & Arabie [1985] Hubert, L., & Arabie, P. (1985).  Comparing partitions.  Journal of classification, 2, 193–218. doi:10.1007/BF01908075.
  * Hubert & Debruyne [2010] Hubert, M., & Debruyne, M. (2010).  Minimum covariance determinant.  WIREs Computational Statistics, 2, 36–43. doi:https://doi.org/10.1002/wics.61.
  * Hubert et al. [2018] Hubert, M., Debruyne, M., & Rousseeuw, P. J. (2018).  Minimum covariance determinant and extensions.  WIREs Computational Statistics, 10, e1421. doi:https://doi.org/10.1002/wics.1421.
  * Hubert et al. [2019] Hubert, M., Rousseeuw, P. J., & den Bossche, W. V. (2019).  Macropca: An all-in-one pca method allowing for missing values as well as cellwise and rowwise outliers.  Technometrics, 61, 459–473. doi:10.1080/00401706.2018.1562989.
  * Hubert et al. [2005] Hubert, M., Rousseeuw, P. J., & Branden, K. V. (2005).  Robpca: A new approach to robust principal component analysis.  Technometrics, 47, 64–79. doi:10.1198/004017004000000563.
  * Hubert et al. [2008] Hubert, M., Rousseeuw, P. J., & Van Aelst, S. (2008).  High-breakdown robust multivariate methods.  Statistical science, 23, 92–119. doi:10.1214/088342307000000087.
  * Kaufman & Rousseeuw [1990] Kaufman, L., & Rousseeuw, P. J. (1990).  Finding groups in data: an introduction to cluster analysis.  New York: John Wiley & Sons Inc.  doi:10.1002/9780470316801.
  * Kent & Tyler [1991] Kent, J. T., & Tyler, D. E. (1991).  Redescending $M$-Estimates of Multivariate Location and Scatter.  The Annals of Statistics, 19, 2102 – 2119. doi:10.1214/aos/1176348388.
  * Kettenring [2006] Kettenring, J. R. (2006).  The practice of cluster analysis.  Journal of classification, 23, 3–30. doi:10.1007/s00357-006-0002-6.
  * Lapidot [2018] Lapidot, I. (2018).  Convergence problems of mahalanobis distance-based k-means clustering.  In 2018 IEEE International Conference on the Science of Electrical Engineering in Israel (ICSEE) (pp. 1–5).  IEEE.  doi:10.1109/ICSEE.2018.8646138.
  * Liski et al. [2014] Liski, E., Nordhausen, K., & Oja, H. (2014).  Supervised invariant coordinate selection.  Statistics: A Journal of Theoretical and Applied Statistics, 4, 711–731. doi:10.1080/02331888.2013.800067.
  * Luo & Li [2016] Luo, W., & Li, B. (2016).  Combining eigenvalues and variation of eigenvectors for order determination.  Biometrika, 103, 875–887. doi:10.1093/biomet/asw051.
  * Luo & Li [2021] Luo, W., & Li, B. (2021).  On order determination by predictor augmentation.  Biometrika, 108, 557–574. doi:10.1093/biomet/asaa077.
  * Mardia et al. [1982] Mardia, K., Bibby, J., & Kent, J. (1982).  Multivariate analysis.  Probability and mathematical statistics.  Acad. Press.  URL: https://books.google.fr/books?id=1nLonQEACAAJ.
  * McLachlan [1992] McLachlan, G. (1992).  Discriminant Analysis and Statistical Pattern Recognition.  Wiley Series in Probability and Statistics.  Wiley.
  * Nordhausen et al. [2011] Nordhausen, K., Oja, H., & Ollila, E. (2011).  Multivariate models and the first four moments.  In Nonparametric Statistics and Mixture Models (pp. 267–287).  Hackensack: World Scientific.  doi:10.1142/9789814340564\\_0016.
  * Nordhausen et al. [2008] Nordhausen, K., Oja, H., & Tyler, D. E. (2008).  Tools for exploring multivariate data: The package ICS.  Journal of Statistical Software, 28, 1–31. doi:10.18637/jss.v028.i06.
  * Nordhausen et al. [2022] Nordhausen, K., Oja, H., & Tyler, D. E. (2022).  Asymptotic and bootstrap tests for subspace dimension.  Journal of Multivariate Analysis, 188, 104830. doi:10.1016/j.jmva.2021.104830.
  * Nordhausen et al. [2017] Nordhausen, K., Oja, H., Tyler, D. E., & Virta, J. (2017).  Asymptotic and bootstrap tests for the dimension of the non-gaussian subspace.  IEEE Signal Processing Letters, 24, 887–891. doi:10.1109/LSP.2017.2696880.
  * Nordhausen & Ruiz-Gazen [2022] Nordhausen, K., & Ruiz-Gazen, A. (2022).  On the usage of joint diagonalization in multivariate statistics.  Journal of Multivariate Analysis, 188, 104844. doi:10.1016/j.jmva.2021.104844.
  * Nordhausen & Tyler [2015] Nordhausen, K., & Tyler, D. E. (2015).  A cautionary note on robust covariance plug-in methods.  Biometrika, 102, 573–588. doi:10.1093/biomet/asv022.
  * Nordhausen & Virta [2019] Nordhausen, K., & Virta, J. (2019).  An overview of properties and extensions of FOBI.  Knowledge-Based Systems, 173, 113–116. doi:10.1016/j.knosys.2019.02.026.
  * Peña et al. [2010] Peña, D., Prieto, F. J., & Viladomat, J. (2010).  Eigenvectors of a kurtosis matrix as interesting directions to reveal cluster structure.  Journal of Multivariate Analysis, 101, 1995–2007. doi:10.1016/j.jmva.2010.04.014.
  * Radojičić & Nordhausen [2019] Radojičić, U., & Nordhausen, K. (2019).  Non-gaussian component analysis: Testing the dimension of the signal subspace.  In Workshop on Analytical Methods in Statistics (pp. 101–123).  Springer.  doi:10.1007/978-3-030-48814-7_6.
  * Radojičić et al. [2021] Radojičić, U., Nordhausen, K., & Virta, J. (2021).  Large-sample properties of unsupervised estimation of the linear discriminant using projection pursuit.  Electronic Journal of Statistics, 15, 6677–6739. doi:10.1214/21-EJS1956.
  * Raymaekers & Rousseeuw [2022] Raymaekers, J., & Rousseeuw, P. (2022).  cellWise: Analyzing Data with Cellwise Outliers.  URL: https://CRAN.R-project.org/package=cellWise r package version 2.4.0.
  * Rousseeuw [1985] Rousseeuw, P. (1985).  Multivariate estimation with high breakdown point.  In Mathematical statistics and applications, Vol. B (Bad Tatzmannsdorf, 1983) (pp. 283–297).  Reidel, Dordrecht.
  * Rousseeuw & Van Driessen [1999] Rousseeuw, P. J., & Van Driessen, K. (1999).  A fast algorithm for the minimum covariance determinant estimator.  Technometrics, 41, 212–223. doi:10.1080/00401706.1999.10485670.
  * Ruiz-Gazen [1996] Ruiz-Gazen, A. (1996).  A very simple robust estimator of a dispersion matrix.  Computational Statistics & Data Analysis, 21, 149–162. doi:https://doi.org/10.1016/0167-9473(95)00009-7.
  * Schoonees et al. [2015] Schoonees, P. C., van de Velden, M., & Groenen, P. J. F. (2015).  Constrained dual scaling for detecting response styles in categorical datas.  Psychometrika, 80, 968–994. doi:10.1007/s11336-015-9458-9.
  * Soete & Carroll [1994] Soete, G. D., & Carroll, J. D. (1994).  K-means clustering in a low-dimensional euclidean space.  In New approaches in classification and data analysis (pp. 212–219).  Springer.  doi:10.1007/978-3-642-51175-2_24.
  * Stahel & Mächler [2009] Stahel, W., & Mächler, M. (2009).  Comment on “invariant co-ordinate selection”.  Journal of the Royal Statistical Society B, 71.
  * Timmerman et al. [2013] Timmerman, M. E., Ceulemans, E., De Roover, K., & Van Leeuwen, K. (2013).  Subspace k-means clustering.  Behavior research methods, 45, 1011–1023. doi:10.3758/s13428-013-0329-y.
  * Todorov [2007] Todorov, V. (2007).  Robust selection of variables in linear discriminant analysis.  Statistical Methods and Applications, 15, 395–407.
  * Tyler [2010] Tyler, D. E. (2010).  A note on multivariate location and scatter statistics for sparse data sets.  Statistics & Probability Letters, 80, 1409–1413. doi:10.1016/j.spl.2010.05.006.
  * Tyler et al. [2009] Tyler, D. E., Critchley, F., Dümbgen, L., & Oja, H. (2009).  Invariant co-ordinate selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71, 549–592. doi:10.1111/j.1467-9868.2009.00706.x.
  * van de Velden et al. [2019a] van de Velden, M., D’Enza, A. I., & Yamamoto, M. (2019a).  Special feature: dimension reduction and cluster analysis.  Behaviormetrika, 46, 239–241. doi:10.1007/s41237-019-00092-6.
  * van de Velden et al. [2019b] van de Velden, M., Iodice D’Enza, A., & Markos, A. (2019b).  Distance-based clustering of mixed data.  Wiley Interdisciplinary Reviews: Computational Statistics, 11, e1456.
  * Venables & Ripley [2002] Venables, W. N., & Ripley, B. D. (2002).  Modern Applied Statistics with S.  (4th ed.).  New York: Springer.  URL: https://www.stats.ox.ac.uk/pub/MASS4/ iSBN 0-387-95457-0.
  * Vichi & Kiers [2001] Vichi, M., & Kiers, H. A. (2001).  Factorial k-means analysis for two-way data.  Computational Statistics & Data Analysis, 37, 49–64. doi:10.1016/S0167-9473(00)00064-5.
  * Vichi & Saporta [2009] Vichi, M., & Saporta, G. (2009).  Clustering and disjoint principal component analysis.  Computational Statistics & Data Analysis, 53, 3194–3208. doi:10.1016/j.csda.2008.05.028.
  * Xu & Wunsch [2005] Xu, R., & Wunsch, D. (2005).  Survey of clustering algorithms.  IEEE Transactions on neural networks, 16, 645–678.

##  Appendix A Cluster settings in the simulation study

We consider the following 22 combinations of the mixture weights
$\epsilon_{1},\dots,\epsilon_{q}$ in the simulation study from Section 4. We
refer to the different cluster settings using the notation
$100\epsilon_{1}$–…–$100\epsilon_{q}$:

  * 1.

$q=2$ clusters: 50–50, 55–45, 60–40, 65–35, 70–30, 75–35, 80–20, 85–15, 90–10,
95–5,

  * 2.

$q=3$ clusters: 33–33–33, 30–40–30, 20–50–30, 10–50–40, 10–60–30, 10–70–20,
10–80–10,

  * 3.

$q=5$ clusters: 20–20–20–20–20, 10–20–20–20–30, 10–10–20–20–40,
10–10–10–30–40, 10–10–20–30–30.

The dimension reduction strategies, which are the main focus of this paper,
are evaluated in all 22 settings in terms of the discriminatory power of the
selected components.

However, to keep computation time low, we apply the clustering methods on the
obtained components only for the following subset of the cluster settings:

  * 1.

$q=2$ clusters: 50–50, 70–30, 80–20, 90–10, 95–5,

  * 2.

$q=3$ clusters: 33–33–33, 20–50–30, 10–80–10,

  * 3.

$q=5$ clusters: 20–20–20–20–20, 10–10–20–20–40.

Hence, results for the clustering performance (see Section 4.2.2 and B.2) are
based on those 10 cluster settings.

##  Appendix B Additional simulation results

###  B.1 Overall comparison of methods in setting without outliers

To establish which methods should be discussed in more detail, we take a look
at their performance in the baseline scenario without outliers. For all
investigated dimension reduction methods and component selection criteria,
Figure 11 provides boxplots of the discriminatory power of the selected
components across the different cluster settings and simulation runs. The
discriminatory power is thereby measured by the $\eta^{2}$.

Figure 11: Discriminatory power as measured by the $\eta^{2}$ of the
investigated dimension reduction methods with different component selection
criteria. Boxplots display the results across the different cluster settings
with $100$ simulation runs each.

We clearly observe that both the highest median and the lowest variability in
discriminatory power is obtained with ICS using the scatter pairs
$\mathrm{LCOV}-\mathrm{COV}$ and $\mathrm{TCOV}-\mathrm{COV}$, followed by the
robust scatter pair $\mathrm{TCOV}-\mathrm{UCOV}$. For those scatter pairs,
there is not much difference in the performance of different component
selection criteria, which are all close to the oracle criterion. For
illustration purposes, the first row of Figure 12 visualizes the shapes of
those scatters as ellipses based on a sample from a bivariate Gaussian mixture
distribution with three balanced components. $\mathrm{LCOV}$ and
$\mathrm{TCOV}$ nicely capture the local structure while $\mathrm{COV}$ and
$\mathrm{UCOV}$ recover the global structure. Suitable combinations therefore
lead to highlighting the structure of interest.

Figure 12: Shapes of different scatter matrices for a sample from a bivariate
Gaussian mixture distribution with three balanced components. In the first
row, the local scatters (together with $\mathrm{COV}$) are shown in the left
panel and the global scatters in the right panel. In the second row,
$\mathrm{MCD}_{\alpha}$ is shown in the left panel and
$\mathrm{RMCD}_{\alpha}$ in the right panel for different values of $\alpha$
(together with $\mathrm{COV}$).

Going back to the results in Figure 11, we find for PCA that the $80\%$
criterion yields components with a higher median and lower variability in
discriminatory power than the $k-1$ criterion. Using the $80\%$ criterion,
$\mathrm{RMCD}_{0.75}$ clearly outperforms the sample covariance matrix
$\mathrm{COV}$, as the median discriminatory power is higher while variability
is lower. However, ICS with the aforementioned three scatter pairs clearly
outperforms PCA with $\mathrm{RMCD}_{0.75}$.

Overall, ICS with a scatter pair involving the $\mathrm{MCD}$ performs
surprisingly well, but the results depend strongly on the parameter $\alpha$.
When paired with $\mathrm{COV}$, somewhat smaller values of $\alpha$ tend to
perform better in general. This is expected since a smaller subset size
increases the chances that the shape of the MCD estimate corresponds to the
local covariance structure within one of the clusters versus the global
covariance structure captured by $\mathrm{COV}$ [cf. the discussion in Hennig,
2009]. Another important question for $\mathrm{MCD}$ scatters is whether to
use the raw estimator $\mathrm{MCD}_{\alpha}$ or the reweighted estimator
$\mathrm{RMCD}_{\alpha}$. Using $\mathrm{MCD}_{\alpha}$ gives (almost) full
control over the subset size, which can be desirable in a clustering context
to capture the local covariance structure. While scatter pairs with
$\mathrm{RMCD}_{\alpha}$ in general yield a higher median $\eta^{2}$ than
scatter pairs with $\mathrm{MCD}_{\alpha}$, small differences in the median
are not that relevant as long as the $\eta^{2}$ is large enough. On the other
hand, scatter pairs with $\mathrm{RMCD}_{\alpha}$ fail to pick up the relevant
structure in more instances compared to their counterparts with
$\mathrm{MCD}_{\alpha}$, i.e., there are more points with small $\eta^{2}$ in
Figure 11, which is most visible for the oracle criterion. For a rather large
subset size (e.g., $\alpha=0.75$), using $\mathrm{RMCD}_{\alpha}$ can be
considered detrimental to performance according to our results. Nevertheless,
for a small subset size (e.g., $\alpha=0.1$), using $\mathrm{RMCD}_{\alpha}$
rather than $\mathrm{MCD}_{\alpha}$ can be beneficial to prevent that the
estimate is based on too few observations. This could be of particular
relevance when $n$ is small, although we emphasize that we did not study such
a setting. The second row of Figure 12, which visualizes the shapes of
$\mathrm{MCD}_{\alpha}$ and $\mathrm{RMCD}_{\alpha}$ in our illustrative
example of a bivariate Gaussian mixture, provides intuition for the above
discussion by exposing the dramatic impact of the parameter $\alpha$ and the
reweighting step on finding the local structure.

In Figure 11, the scatter pairs $\mathrm{MLC}-\mathrm{COV}$ and
$\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$ clearly perform worse than
the other scatter pairs for ICS in terms of discriminatory power, as they have
a lower median $\eta^{2}$ and fail in many instances. The first row of Figure
12 further illustrates that those scatter pairs are not the best choices for
finding the local structure. Nevertheless, since
$\mathrm{COV}-\mathop{\mathrm{COV}_{4}}\nolimits$ is likely the most widely
used scatter pair for ICS, we find it relevant to discuss its behavior in the
main text.

Regarding the component selection criteria for ICS, keep in mind that in this
paper, we only consider clustering methods that take the number of clusters
$k$ as an input. The med and var criteria, which are based on selecting $k-1$
components, are therefore conceptually appealing and, across the different
scatter pairs, lead to the best overall performance in discriminatory power.

###  B.2 Clustering performance in selected cluster settings

Figures 13, 14, and 15 display the clustering results for the setting without
outliers, with 2% outliers, and with 5% outliers, respectively. The figures
show boxplots of the ARI of the five clustering methods in the selected
cluster settings, applied after the best performing dimension reduction
methods for ICS and PCA.

Figure 13: Clustering performance as measured by the ARI after the best
performing dimension reduction methods for ICS and PCA, respectively, using
the corresponding best component selection criteria. Boxplots display the
results for the different cluster settings over $100$ simulation runs in the
setting without outliers. Results for different clustering methods are shown
in separate rows, and results for different dimension reduction methods are
shown in separate columns. Figure 14: Clustering performance as measured by
the ARI after the best performing dimension reduction methods for ICS and PCA,
respectively, using the corresponding best component selection criteria.
Boxplots display the results for the different cluster settings over $100$
simulation runs in the setting with 2% outliers. Results for different
clustering methods are shown in separate rows, and results for different
dimension reduction methods are shown in separate columns. Figure 15:
Clustering performance as measured by the ARI after the best performing
dimension reduction methods for ICS and PCA, respectively, using the
corresponding best component selection criteria. Boxplots display the results
for the different cluster settings over $100$ simulation runs in the setting
with 5% outliers. Results for different clustering methods are shown in
separate rows, and results for different dimension reduction methods are shown
in separate columns.

In Figure 13 for the setting without outliers, we observe a clear improvement
in clustering performance with ICS for PAM and $k$means. Specifically,
$\mathrm{TCOV}-\mathrm{COV}$ and $\mathrm{TCOV}-\mathrm{UCOV}$ yield a nearly
perfect ARI in all cluster settings. $\mathrm{LCOV}-\mathrm{COV}$ performs
similarly in most cluster settings, but fails on occasion for perfectly
balanced clusters (50-50). For t$k$means, we find a substantial improvement
only for certain cluster settings (50-50, 20-50-30). On the other hand, mclust
and rmclust, perform almost perfectly on the observed data. While ICS keeps
this near-perfect performance in most settings, it can increase failures in
some settings (in particular for rmclust). Finally, PCA with
$\mathrm{RMCD}_{0.75}$ is detrimental to clustering performance, yielding a
lower ARI or larger variability than on the observed data in almost every
setting.

From Figures 14 and 15, we see that the main findings remain valid in the
settings with outliers. For all clustering methods, ICS with
$\mathrm{TCOV}-\mathrm{COV}$ and $\mathrm{LCOV}-\mathrm{COV}$ clearly yields
the best performance except for rmclust, where performance improves for some
settings but deteriorates in others. The simple clustering methods $k$means
and PAM after ICS even outperform rmclust on observed data in the settings
with $q=5$ clusters. Finally, PCA with $\mathrm{RMCD}_{0.75}$ still has a
detrimental effect on clustering performance.

###  B.3 Illustration for the barrow wheel distribution

In the simulation study of Section 4, we generate the data from a Gaussian
mixture model due to its simplicity and popularity. To provide an illustration
of tandem clustering using ICS in a non-elliptical framework, we simulate a
data set from the barrow wheel distribution, which was introduced in Hampel et
al. [1986] and used by Stahel & Mächler [2009] as an example in the discussion
of Tyler et al. [2009]. A $d$-variate random vector $\mathbf{X}$ follows a
barrow wheel distribution if,

${\mathbf{X}}\sim\mathbf{O}((1-\epsilon){\cal
N}(\mathbf{0},\mathrm{diag}(\sigma_{1}^{2},1,\dots,1))+\epsilon H),$

with $H$ so that $h=(h_{1},\mathbf{h}_{2}^{\prime})^{\prime}\sim H$ implies
$h_{1}\sim{\chi}_{d-1}$ and $h_{2}\sim{\cal N}(0,\sigma_{2}^{2}I_{d-1})$. The
$d\times d$ matrix $\mathbf{O}$ can be used to rotate and rescale the vector.

Figure 16: Scatter plot matrix of a simulated data set from the barrow wheel
distribution.

Figure 16 shows our simulated data set of size $n=1000$ from the barrow wheel
distribution with $d=3$, $\sigma_{1}=0.1$, $\sigma_{2}=0.2$ and
$\epsilon=0.2$. Table 2 reports the ARI values of $k$means, mclust, and
rmclust on this data set, as well as after performing ICS with
$\mathrm{TCOV}-\mathrm{COV}$. The approaches based on a Gaussian mixture model
(mclust and rmclust) fail to detect the clusters in the simulated data set
with ARI values around 0.34, while $k$means performs even worse with an ARI
close to 0. However, using tandem clustering with ICS and TCOV-COV recovers
the cluster structure in the first coordinate, as illustrated in Figure 17.
Applying the clustering methods to the first invariant coordinate greatly
improves the identification of the three clusters and yields ARI values
ranging from 0.835 for $k$means to about 0.96 for mclust an rmclust.

|  |  |  | ICS + | ICS + | ICS +
---|---|---|---|---|---|---
|  $k$means | mclust | rmclust |  $k$means | mclust | rmclust
ARI | 0.030 | 0.346 | 0.339 | 0.835 | 0.958 | 0.955
Table 2: Adjusted Rand Index for different clustering techniques on the
simulated data from the barrow wheel distribution. Here, ICS is applied using
the scatter pair $\mathrm{TCOV}-\mathrm{COV}$. Figure 17: Scatter plot
matrices of the invariant coordinates obtained by ICS with
$\mathrm{TCOV}-\mathrm{COV}$ on the simulated data set from the barrow wheel
distribution.

##  Appendix C Additional results for the empirical applications

For the three empirical applications discussed in Section 5, Figure 18
compares the five considered clustering methods to two variants of $k$means
with a built-in approach for subspace estimation, namely reduced $k$means
[Soete & Carroll, 1994] and factorial $k$means [Vichi & Kiers, 2001]. We focus
on applying the five previously considered methods to the the observed data,
as well as after ICS with the med criterion for $\mathrm{TCOV}-\mathrm{COV}$
and $\mathrm{MCD}_{0.5}-\mathrm{COV}$. Clearly, reduced $k$means and the
factorial $k$means do not improve the results and are always outperformed by
at least one version of tandem clustering with ICS.

Figure 18: Adjusted Rand Index (ARI) in the empirical applications for
different clustering methods on the observed data, as well as after ICS with
the med criterion for $\mathrm{TCOV}-\mathrm{COV}$ and
$\mathrm{MCD}_{0.5}-\mathrm{COV}$. Results for different data sets are shown
in separate columns.