# Generalized Spherical Principal Component Analysis

Sarah Leyder Corresponding author. Email<EMAIL_ADDRESS>Department of Mathematics, University of Antwerp, Belgium  Jakob Raymaekers
Department of Quantitative Economics, Maastricht University, The Netherlands
Tim Verdonck  Department of Mathematics, University of Antwerp, Belgium
IDLab, Imec Antwerp, Belgium

###### Abstract

Outliers contaminating data sets are a challenge to statistical estimators.
Even a small fraction of outlying observations can heavily influence most
classical statistical methods. In this paper we propose generalized spherical
principal component analysis, a new robust version of principal component
analysis that is based on the generalized spatial sign covariance matrix.
Supporting theoretical properties of the proposed method including influence
functions, breakdown values and asymptotic efficiencies are studied, and a
simulation study is conducted to compare our new method to existing methods.
We also propose an adjustment of the generalized spatial sign covariance
matrix to achieve better Fisher consistency properties. We illustrate that
generalized spherical principal component analysis, depending on a chosen
radial function, has both great robustness and efficiency properties in
addition to a low computational cost.

Keywords: Principal component analysis, Robustness, Influence functions,
Efficiency, Breakdown value

##  1 Introduction

A well-known and frequently used technique to analyze the structure of data
sets is principal component analysis (PCA). The objective of this technique is
to obtain a new set of uncorrelated variables while preserving as much as
possible of the variation present in the original data. These new variables
are constructed as linear combinations of the original variables.
Equivalently, they are obtained by projecting the original data on the PCA
loading vectors of the data set, i.e. the directions in which the data has the
greatest variability. PCA is a key building block in statistical data analysis
and is widely used as a first step in clustering, discriminant analysis and
regression.

In classical PCA (CPCA), the principal components can be calculated through a
spectral decomposition of the covariance matrix. However, it is widely known
that this matrix is very sensitive to outliers and potentially heavily
influenced by anomalous observations. As a result, the directions of greatest
variability are easily attracted towards these outliers, distorting the output
of PCA. In order to avoid this, robust PCA methods have been developed which
are resistant to such outlying observations. There exist many different
approaches to robust principal component analysis.

One approach is to use a spectral decomposition of a robust estimate of the
covariance matrix. For affine equivariant covariance matrices, this approach
was studied by Campbell (1980) and Boente (1987) who used M-estimators, which
unfortunately cannot withstand many outliers. Croux and Haesbroeck (2000)
revisited this approach, suggesting instead to use high-breakdown estimators
of location and scatter such as S-estimators (Davies, 1987; Rousseeuw and
Yohai, 1984) or the MCD (Rousseeuw, 1984). In particular, they derive general
expressions for the influence functions and efficiencies of the resulting
eigenvector and eigenvalue estimates. A drawback of the approach is that many
highly robust covariance estimators are computationally demanding.
Additionally, some of them can only be computed when the number of samples is
(substantially) larger than the dimension.

A second approach works incrementally starting from the principal component
corresponding with the largest projected variance. Each new principal
component is then estimated as a maximizer of the projected robust variance
conditional on being orthogonal to the already estimated components. This
approach is detailed in Li and Chen (1985); Croux and Ruiz-Gazen (2005); Croux
et al. (2007). The approach works well when a relatively small number of
principal components is required. However, it can also be computationally
demanding as the number of projections needed should increase rapidly with the
dimension of the data to guarantee a stable performance. A combination of the
two approaches mentioned above was used by Hubert et al. (2005) to develop the
ROBPCA algorithm. ROBPCA often outperforms the projection-pursuit based
methods as well as the covariance-based methods. It remains fairly slow to
compute, especially on larger data sets.

To mitigate the computational burden of the previously mentioned approaches to
robust principal component analysis, one elegant and popular approach for
robust PCA is spherical principal component analysis (SPCA), introduced
independently by Marden (1999) and Locantore et al. (1999). Spherical PCA
starts by projecting the centered data onto a unit sphere before performing
classical principal component analysis on this transformed data set. It was
studied by, among others, Visuri et al. (2001); Taskinen et al. (2012); Croux
et al. (2002). This procedure is equivalent to performing PCA on the spatial
sign covariance matrix (SSCM):

$\Sigma_{\text{SSCM}}(X)=\mathbb{E}\left[\frac{(X-\mu)(X-\mu)^{T}}{\left\lVert
X-\mu\right\rVert^{2}}\right]$

with $X$ a $p$-variate random variable, $\mu$ the location of the distribution
of $X$ and $\left\lVert.\right\rVert$ the Euclidean norm. Under mild
assumptions on the underlying distribution, the SSCM is a Fisher consistent
estimator of the eigenvectors and it preserves the order of the eigenvalues.
The SSCM was studied in detail in Magyar and Tyler (2014); Dürre et al. (2014,
2016); Boente et al. (2019).

In Raymaekers and Rousseeuw (2019) a generalisation to the SSCM was
introduced, namely the generalized spatial sign covariance matrix (GSSCM).
They found that the SSCM is part of a large class of orthogonally equivariant
scatter estimates, namely the generalized spatial sign covariance matrices.
Whereas in the SSCM all centered data points $x_{i}$ are given the weight
$1/\left\lVert x_{i}\right\rVert$, the generalized SSCM assigns different
weights, depending on the distribution of the random variable:

$\displaystyle\Sigma_{g_{X}}(X)$
$\displaystyle=\mathbb{E}\left[g_{X}(X-\mu)g_{X}(X-\mu)^{T}\right]$ (1)
$\displaystyle\text{with}\quad g_{X}(t)=t\ \xi_{X}(\left\lVert t\right\rVert)$

where $\xi_{X}:\mathbb{R}^{+}\to\mathbb{R}^{+}$ is the radial function. By
using the Euclidean norm, the GSSCM becomes an orthogonally equivariant
scatter estimator. In Raymaekers and Rousseeuw (2019), it is shown that the
GSSCM inherits the consistency properties of the SSCM in that it is a Fisher
consistent estimator of the eigenvectors and preserves the ranks of the
eigenvalues under the same assumptions.

A new robust method for principal component analysis emerges when we combine
the idea of spherical PCA with the GSSCM. Instead of computing the principal
components from the SSCM, we can compute them from the generalized SSCM. We
refer to this method as generalized spherical principal component analysis
(GSPCA).

In this paper, we introduce and investigate the GSPCA method from a robustness
perspective. The rest of the paper is organized as follows. Section 2 formally
introduces GSPCA. Section 3 covers several theoretical properties of the new
method, such as breakdown values, influence functions and asymptotic
variances. Additionally, we further improve our approach by adapting the GSSCM
estimator to get Fisher consistency for the eigenvalues. In Section 4 an
extensive simulation study compares our method to SPCA (Locantore et al.,
1999) and ROBPCA (Hubert et al., 2005), the state-of-the-art method. Section 5
finalizes by applying GSPCA to real life data sets.

##  2 Generalized spherical principal component analysis

As discussed in the introduction, the GSPCA method computes the principal
component directions as the eigenvectors of the GSSCM given in Equation (1).
For a finite, $p$-variate data set $\mathbf{X}$ containing
$\\{x_{1},\dots,x_{n}\\}$ the sample GSSCM becomes

$\displaystyle S$
${}_{g_{\mathbf{X}}}(\mathbf{X})=\frac{1}{n}\sum_{i=1}^{n}\xi_{\mathbf{X}}^{2}(\left\lVert
x_{i}-T(\mathbf{X})\right\rVert)\
(x_{i}-T(\mathbf{X}))(x_{i}-T(\mathbf{X}))^{T}=\frac{1}{n}\
g_{\mathbf{X}}(\mathbf{X})^{T}g_{\mathbf{X}}(\mathbf{X})$

where $T$ is a (orthogonally equivariant) location estimator for the center of
the data matrix $\mathbf{X}\in\mathbb{R}^{n\times p}$. In Raymaekers and
Rousseeuw (2019), it is suggested to use the $k$-step least trimmed squares
(LTS) estimator as a robust location. This estimator starts from the spatial
median, but adds some extra successive steps to improve robustness against
outliers. The $k$-step LTS estimator has a breakdown value of
$\left\lfloor{(n+1)/2}\right\rfloor/n$.
In this paper we will consider the following five radial functions suggested
in Raymaekers and Rousseeuw (2019), apart from the evident radial functions
$\xi(r)=1$ and $\xi(r)=1/r$, respectively corresponding to the classical
covariance matrix and the SSCM.

  1. 1.

Winsorizing (Winsor):

$\displaystyle\qquad\xi_{X}(r)$ $\displaystyle=\begin{cases}1&\quad\text{if
}r\leq Q_{2}\\\ Q_{2}/r&\quad\text{if }Q_{2}<r\end{cases}$ (2)

  2. 2.

Quadratic Winsor (Quad):

$\displaystyle\qquad\xi_{X}(r)$ $\displaystyle=\begin{cases}1&\quad\text{if
}r\leq Q_{2}\\\ Q_{2}^{2}/r^{2}&\quad\text{if }Q_{2}<r\end{cases}$ (3)

  3. 3.

Ball:

$\displaystyle\qquad\xi_{X}(r)$ $\displaystyle=\begin{cases}1&\quad\text{if
}r\leq Q_{2}\\\ 0&\quad\text{if }Q_{2}<r\end{cases}$ (4)

  4. 4.

Shell

$\displaystyle\qquad\xi_{X}(r)$ $\displaystyle=\begin{cases}0&\quad\text{if
}r<Q_{1}\\\ 1&\quad\text{if }Q_{1}\leq r\leq Q_{3}\\\ 0&\quad\text{if
}Q_{3}<r\end{cases}$ (5)

  5. 5.

Linearly Redescending (LR):

$\displaystyle\qquad\xi_{X}(r)$ $\displaystyle=\begin{cases}1&\ \text{if
}r\leq Q_{2}\\\ \frac{Q_{3}^{*}-r}{Q_{3}^{*}-Q_{2}}&\ \text{if }Q_{2}<r\leq
Q_{3}^{*}\\\ 0&\ \text{if }Q_{3}^{*}<r\end{cases}$ (6)

The cutoffs $Q_{1},Q_{2},Q_{3}\text{ and }Q_{3}^{*}$ are robust estimates of
the first, second and third quartile of the distribution of the Euclidean
distances. We will come back to these in the next section.

For a given data set $\mathbf{X}$ we can now use the spectral decomposition of
the sample GSSCM to obtain the GSPCA loading vectors:

$S_{g_{\mathbf{X}}}(\mathbf{X})=V_{g}\Lambda_{g}V_{g}^{T}$

Here the matrix $\Lambda_{g}$ is the diagonal matrix containing the
eigenvalues $\lambda_{g,i}$ of $S_{g_{\mathbf{X}}}(\mathbf{X})$ in descending
order ($\lambda_{g,1}>\dots>\lambda_{g,p})$ and the columns of the matrix
$V_{g}$ consist of the corresponding eigenvectors $v_{g,i}$. The GSPCA loading
vectors then coincide with these eigenvectors and the principal components
correspond to the data projected onto these vectors.

Alternatively we can also state that the principal components are the
uncorrelated vectors maximizing the variance of the transformed data. In
particular, the $i$-th loading vector $a_{i}$ can be defined as

$\displaystyle a_{i}=\operatorname*{arg\,max}_{a\in
A}\left\\{\text{var}(g_{\mathbf{X}}(\mathbf{X})a)\right\\}$
$\displaystyle\text{ where }A$ $\displaystyle=\\{a\in\mathbb{R}^{p}\
|\left\lVert a\right\rVert_{2}=1\text{ and }a\cdot a_{j}=0\text{ for all
}j=1,\ldots,{i-1}\\}$

It should be clear from the definition that GSPCA generally has a fairly low
computational cost. The only uncertain factor is the computation and
evaluation of the radial functions $g$. In the cases considered above however,
it requires the computation of the Euclidean norms of the observations
followed by robust estimates of the quantiles of these norms. The first can be
done in $\mathcal{O}(np)$ time, and the latter in $\mathcal{O}(n)$. Therefore,
GSPCA has a computational complexity given by
$\mathcal{O}(n+np+np^{2}+p^{3})=\mathcal{O}(np^{2}+p^{3})$, the same as
classical PCA.

Additionally, like CPCA, GSPCA has the property that it can be computed
through the singular value decomposition (SVD) as well as the spectral
composition of the GSSCM. In order to use the SVD for GSPCA, we first need to
transform the observations using the function $g_{X}$ of Equation (1) after
which we can apply SVD to the transformed data $g_{\mathbf{X}}(\mathbf{X})$.
More precisely, we first compute the singular value decomposition

$g_{\mathbf{X}}(\mathbf{X})=UDV^{T}$

where the orthogonal matrices $U$ and $V$ respectively contain the left and
right singular vectors of $g_{\mathbf{X}}(\mathbf{X})$ and the diagonal matrix
$D$ contains the corresponding singular values $\delta_{g,i}$. Next it can
easily be shown that the columns $v_{g,i}$ of $V$ correspond to the GSPCA
loading vectors and the values $\delta_{g,i}^{2}/n$ correspond to the GSPCA
eigenvalues:

$\displaystyle S_{g_{\mathbf{X}}}(\mathbf{X})$ $\displaystyle=\frac{1}{n}\
g_{\mathbf{X}}(\mathbf{X})^{T}g_{\mathbf{X}}(\mathbf{X})=\frac{1}{n}\
VDU^{T}UDV^{T}$ $\displaystyle=\frac{1}{n}\ VDDV^{T}$
$\displaystyle=V\text{diag}(\delta_{g,1}^{2}/n,\dots,\delta_{g,p}^{2}/n)V^{T}=V_{g}\Lambda_{g}V_{g}^{T}$

The SVD procedure has the advantage that the computational cost of GSPCA can
be lowered if the required number of principal components is known in advance.
In this scenario we can use truncated singular value decomposition on the
GSSCM to calculate only the first $k$ singular vectors, which correspond with
the first $k$ eigenvectors of the matrix $S_{g_{\mathbf{X}}}(\mathbf{X})$.
This not only avoids the calculation of the remaining principal components,
but also does not require the explicit construction of the $p\times p$
covariance estimate $S_{g_{\mathbf{X}}}(\mathbf{X})$.

##  3 Theoretical properties

In this section we study the robustness and efficiency of GSPCA. To evaluate
robustness, the breakdown value of our method is calculated and influence
functions are derived. The latter are then used to evaluate the efficiency
through the computation of the asymptotic variances. We use the obtained
theoretical properties to objectively compare different radial functions and
explain their statistical behavior.

###  3.1 Breakdown value

The breakdown value is a global measure for robustness defined by Hampel based
on an idea of Hodges (Hampel et al., 1986). For a given estimator, it is
defined as the smallest fraction of observations in the data set that needs to
be changed to carry the estimate arbitrarily far. More specifically, for a
scatter estimator $S$, it is defined as the minimal amount of contamination
required to make the largest eigenvalue $\lambda_{1}$ arbitrarily large
(explosion) or the smallest eigenvalue $\lambda_{p}$ arbitrarily close to zero
(implosion):

$\displaystyle\varepsilon(S,\mathbf{X})=\min$
$\displaystyle\left\\{\frac{m}{n}:\sup_{\mathbf{X}_{m}^{*}}\
\max[\lambda_{1}(S(\mathbf{X}_{m}^{*})),\lambda_{p}^{-1}(S(\mathbf{X}_{m}^{*}))]=\infty\right\\}$

with $\mathbf{X}_{m}^{*}$ the data set $\mathbf{X}$ where $m$ observations
have been replaced.

While in many applications, implosion of the covariance matrix is undesirable,
in the context of PCA, it is the explosion breakdown that is more relevant.
After all, if the dimension is equal to or larger than the sample size, we
will have that the smallest eigenvalue is zero and would thus have “breakdown”
if we were to use the definition above. In Raymaekers and Rousseeuw (2019) it
is shown that the breakdown value of the GSSCM is
$\left\lfloor{(n-p+1)/2}\right\rfloor/n$. This seems unsatisfactory, since it
is in contrast with the SSCM which has a higher breakdown value, namely
$\left\lfloor{(n+1)/2}\right\rfloor/n$. Fortunately, the lower breakdown value
of the GSSCM is due to the implosion breakdown value of the “hard
redescending” radial functions such as Ball, Shell and LR, in combination with
the estimation of the cutoffs $Q_{1},Q_{2},Q_{3}\text{ and }Q_{3}^{*}$ used in
these radial functions.

If we are no longer concerned with the implosion breakdown value, we can take
a slightly different approach. First, we estimate the cutoffs using the
estimators

$\displaystyle Q_{1}$ $\displaystyle=\Big{[}\text{med}_{i}\left(\left\lVert
x_{i}-T(\mathbf{X})\right\rVert^{2/3}\right)-\text{MAD}_{i}\left(\left\lVert
x_{i}-T(\mathbf{X})\right\rVert^{2/3}\right)\Big{]}^{3/2}$ $\displaystyle
Q_{2}$ $\displaystyle=\text{med}_{i}(\left\lVert
x_{i}-T(\mathbf{X})\right\rVert)$ $\displaystyle Q_{3}$
$\displaystyle=\Big{[}\text{med}_{i}\left(\left\lVert
x_{i}-T(\mathbf{X})\right\rVert^{2/3}\right)+\text{MAD}_{i}\left(\left\lVert
x_{i}-T(\mathbf{X})\right\rVert^{2/3}\right)\Big{]}^{3/2}$ $\displaystyle
Q_{3}^{*}$ $\displaystyle=\Big{[}\text{med}_{i}\left(\left\lVert
x_{i}-T(\mathbf{X})\right\rVert^{2/3}\right)+1.4826\cdot\text{MAD}_{i}\left(\left\lVert
x_{i}-T(\mathbf{X})\right\rVert^{2/3}\right)\Big{]}^{3/2}$

In these definitions med and MAD are the median and the median absolute
deviation. Unlike in Raymaekers and Rousseeuw (2019), these estimators for the
cutoffs no longer depend on the order statistic
$h=\left\lfloor{(n+p+1)/2}\right\rfloor$, which was required to prove the
implosion breakdown value of the GSSCM. The dependence on the $h$-th order
statistic is undesirable in the context of PCA, as the value of $p$ is often
close to $n$ or even considerably larger than $n$ which would lead to $h>n$,
making the cutoffs meaningless.

The theorem below states that for our adapted GSSCM, with cutoffs based on the
median and MAD, we obtain an explosion breakdown value independent of $p$,
equal to that of the SSCM. The proof can be found in the Supplementary
material.

######  Theorem 1 (Explosion breakdown value).

Given $\mathbf{X}=\\{x_{1},...,x_{n}\\}$ an $n\times p$ dimensional data set
and a location estimator $T(\mathbf{X})$ with a breakdown value of at least
$\left\lfloor{(n+1)/2}\right\rfloor/n$. Suppose that

  1. 1.

The radial function takes values in $[0,1]$

  2. 2.

$\forall\text{ data sets }\mathbf{X}:\\#\\{x_{i}:\xi(\left\lVert
x_{i}-T(\mathbf{X})\right\rVert)=1\\}\geq\left\lfloor{(n+1)/2}\right\rfloor$

  3. 3.

$\forall t:\ \left\lVert g(t)\right\rVert=\left\lVert
t\right\rVert\xi(\left\lVert
t\right\rVert)\leq\text{med}_{i}(d_{i})+1.4826\cdot\text{MAD}_{i}(d_{i})$

Then the explosion breakdown value $\varepsilon$ of the GSSCM is
$\left\lfloor{(n+1)/2}\right\rfloor/n$.

Note that the conditions in the theorem are all satisfied for the radial
functions in Equations (2) – (6). In addition, the k-step LTS estimator has a
breakdown value of $\left\lfloor{(n+1)/2}\right\rfloor/n$. Hence, we conclude
that GSPCA based on the GSSCM is robust up to
$\frac{\left\lfloor{(n+1)/2}\right\rfloor}{n}\approx 50\%$ contamination.

###  3.2 Influence functions of the loading vectors

We now consider the influence functions relevant for GSPCA. The influence
function for an estimator $T$ at distribution $F$ is defined as (see Hampel et
al. (1986)):

$\displaystyle\text{IF}(x,T,F)$ $\displaystyle=\lim_{\varepsilon\to
0}\frac{T(F_{\varepsilon,x})-T(F)}{\varepsilon}=\frac{\partial}{\partial\varepsilon}T(F_{\varepsilon,x})\Big{|}_{\varepsilon=0}$
(7)

with $F_{\varepsilon,x}=(1-\varepsilon)F+\varepsilon\Delta_{x}$ the
distribution contaminated by $x$ where $\Delta_{x}$ is the distribution
putting all its mass in $x$.

In contrast to the breakdown value, the influence function is a local measure
of robustness instead of a global one. Informally, it can be interpreted as a
measure for the effect that an infinitesimal small amount of contamination has
on the estimate. Therefore it is a complementary measure and describes a
different aspect of our method.

In the context of PCA, the relevant functionals $T$ in Equation (7) are the
eigenvector and eigenvalue functionals, which we denote by $v_{g,j}$ and
$\lambda_{g,j}$ for $j=1,\dots,p$. In this section specifically, we will study
the eigenvector functionals, therefore the relevant influence function is
given by

$\text{IF}(x,v_{g,j},F)=\lim_{\varepsilon\to
0}\frac{v_{g,j}(F_{\varepsilon,x})-v_{g,j}(F)}{\varepsilon}$

In what follows, we assume that $F$ is a distribution with the center at the
origin and covariance matrix $\Sigma$ which has distinct eigenvalues.
Furthermore, we denote the spectral decomposition of $\Sigma$ by
$\Sigma=V\Lambda V^{T}$ where $V$ is a matrix with the eigenvectors $v_{i}$ in
its columns and $\Lambda$ is a diagonal matrix with the eigenvalues
$\lambda_{1}>\ldots>\lambda_{p}$ on its diagonal.

Note that the functionals $v_{g,j}$ inherit the Fisher consistency properties
of the GSSCM in case of elliptically symmetric distributions, see Raymaekers
and Rousseeuw (2019). We say that a distribution $F_{X}$ has an elliptically
symmetric density if the density can be written as

$f_{X}(x)=det(\Sigma)^{-1/2}\cdot h((x-\mu)^{T}\Sigma^{-1}(x-\mu))$

In this case one thus has that $v_{g,j}(F)=v_{j}$.
Additionally, $v_{g,j}$ is orthogonally equivariant as the GSSCM is
orthogonally equivariant:

$\displaystyle\forall\ U$ $\displaystyle\in\mathbb{R}^{p\times p}\text{
orthogonal matrix:}v_{g,j}(F_{U})=Uv_{g,j}(F)U^{T}$

where $F_{U}$ is the distribution of $UX$ if $X\sim F$.

We now state the expression for the influence functions corresponding to the
loading vectors of GSPCA, for which the proof can be found in Section A.2 of
the Supplementary material:

######  Theorem 2 (Influence functions of the eigenvectors of the GSSCM).

For the influence functions of the eigenvectors resulting from GSPCA we have
the following analytical expression:

$\displaystyle\text{IF}(x,v_{g,j},F)=\sum_{k=1,k\neq
j}^{p}\frac{1}{\lambda_{g,j}-\lambda_{g,k}}\Big{[}$ $\displaystyle
v_{g,k}^{T}\int\left(dge(X)g(X)^{T}\right.+\left.g(X)dge(X)^{T}\right)dF(X)\
v_{g,j}$ $\displaystyle+(v_{g,k}^{T}g(x))(v_{g,j}^{T}g(x))\Big{]}v_{g,k}$ (8)

where
$dge(x)=\frac{\partial}{\partial\varepsilon}g_{\varepsilon}(X)\Big{|}_{\varepsilon=0}$

While the above expression looks rather involved, we can simplify it further
without (much) loss of generality. First, note that we can assume that $F$ has
a diagonal covariance matrix since the GSSCM is orthogonally equivariant. In
that case, the eigenvectors of the covariance matrix are the standard unit
vectors, i.e. $v_{j}=e_{j}$. Additionally, we can assume that $F$ is
elliptically symmetric. This assumption is very natural since it is required
for the Fisher consistency of the eigenvector functionals $v_{g,j}$ of the
GSSCM, yielding $v_{g,j}(F)=v_{j}$. We then obtain the result below, the proof
of which is in Section A.3 of the Supplementary material.

######  Corollary 1.

Assume that $F$ is a centered elliptically symmetric distribution, that it has
a density function and a diagonal covariance matrix $\Sigma$, then the
following holds:

$\text{IF}(x,v_{g,j},F)=\sum_{k=1,k\neq
j}^{p}\frac{1}{\lambda_{g,j}-\lambda_{g,k}}\big{(}g(x)_{k}g(x)_{j}\big{)}v_{k}$

We thus see that, under reasonable assumptions, the influence function of the
eigenvectors of GSPCA reduces to a fairly simple expression. In particular, if
we plug in the identity function for $g$, we obtain the influence function of
the loadings of classical PCA. In particular, we obtain $\sum_{k=1,k\neq
j}^{p}\frac{x_{k}x_{j}}{\lambda_{j}-\lambda_{k}}v_{k}$, corresponding to the
result obtained in Croux and Haesbroeck (2000).

Having obtained analytical forms for the influence functions of the loading
vectors resulting from GSPCA, we will now visualize them to compare the
different radial functions in Equations (2) – (6) and illustrate their
robustness. We consider a bivariate normal distribution,
$F=\mathcal{N}(0_{2},\Sigma)$ with $\Sigma=\text{diag}(1,0.5)$, and plot the
norm of the influence function of the largest eigenvector for the different
radial functions resulting in the 3D-plots shown in Figure 1.

First, we observe that the eigenvectors of the classical covariance matrix are
the only ones with an unbounded influence function, all others are bounded.
Second, all GSSCM radial functions are redescending to zero, except for
Winsor, whose influence function looks like a smoothed version of the
influence function of the SSCM, suggesting that Winsor will attain higher
efficiency than the SSCM. Further, we see that the norms of the influence
functions of LR, Ball and Shell look quite similar, all three demonstrate four
large spikes whereafter each influence function becomes zero. This is due to
their cutoffs at the second or third quantile. Quad only descends to zero in
the limit, making it more robust than Winsor, but not quite as robust as the
redescending radial functions. Overall these results suggest high robustness
for the LR, Ball and Shell radial functions.

(a) classical covariance matrix

(b) SSCM

(c) Winsor

(d) Quad

(e) LR

(f) Ball

(g) Shell

Figure 1: Norm of the influence function of the largest eigenvector at
$\mathcal{N}(0_{2},\text{diag}(1,0.5))$ for different radial functions

Additionally, we study a robustness measure derived from the influence
function, being the gross-error sensitivity (GES). It is defined as follows
for an estimator $T$ at distribution $F$:

$\gamma^{*}(T,F)=\sup_{x}|\text{IF}(x,T,F)|$

The gross-error sensitivity measures the maximal influence an infinitesimally
small amount of contamination can have on an estimator $T$. Therefore it
should be finite for robust estimators and preferably small.

For the loading vectors of GSPCA, we now compute the gross-error sensitivities
per radial function numerically to compare robustness. We again consider the
bivariate normal distribution $\mathcal{N}(0_{2},\text{diag}(1,0.5))$ and
study the first loading vector. We compute the gross-error sensitivity as the
supremum of the norm of the corresponding influence function. This results in
Table 1.

| classical | SSCM | Winsor | Quad
---|---|---|---|---
$\gamma^{*}$ | $\infty$ | 2.914213 | 3.100523 | 6.569927
| LR | Ball | Shell |
| 10.74909 | 29.27686 | 12.42133 |
Table 1: Gross-error sensitivities of the largest eigenvector at
$\mathcal{N}(0_{2},\text{diag}(1,0.5))$

Once again, we observe that CPCA is not robust, its GES is unbounded. In
contrast, GSPCA has bounded GES for each radial function, implying robustness.
The smallest values are attained by SSCM and Winsor. While this may suggest
these are the superior methods, the plots of the influence functions do show
that the Quad, LR, Ball and Shell methods will have a smaller influence
function for most values of the contamination $x$.

###  3.3 Asymptotic variances and asymptotic relative efficiencies of the
loading vectors

We will now use the results on the influence functions to study the asymptotic
variances and efficiencies of GSPCA. This will allow us to compare precision
among the different radial functions. Following Hampel et al. (1986), we
obtain that for well-behaved functionals $T$, the corresponding estimator
$T_{n}=T(F_{n})$ is asymptotically normal:
$\sqrt{n}\left(T_{n}-T(F)\right)\xrightarrow[]{D}\mathcal{N}\left(0,\text{ASV}(T,F)\right)$
where

$\text{ASV}(T,F)=\mathbb{E}_{F}\left[\mbox{IF}(x,T,F)\mbox{IF}(x,T,F)^{T}\right].$

In order to simplify the exposition, we consider the case of a bivariate
normal distribution $\mathcal{N}(0_{2},\text{diag}(1,\gamma))$. Since we are
interested in eigenvectors for PCA, we first calculate the asymptotic variance
(ASV) of the second element of the largest eigenvector, given that it is
important that we estimate this second element close to zero to get a good
estimation of the first eigenvector. This is similar to the approach taken in
Croux et al. (2010). We denote the asymptotic variance of the second element
of the largest eigenvector simply as $\text{ASV}_{g}$ per radial function.

In Section 3.2 we found analytical expressions for the influence functions of
the eigenvectors. Hence, we can use them to calculate the asymptotic
variances.

For the second component of the influence function of the first eigenvector,
which we denote by $v_{g,(1,2)}$, we obtained the following expression:

$\text{IF}(x,v_{g,(1,2)},F)=\frac{1}{\lambda_{g,1}-\lambda_{g,2}}g(x)_{1}g(x)_{2}$

Using this, we obtain the asymptotic variance:

$\displaystyle\text{ASV}_{g}(F)=\frac{1}{(\lambda_{g,1}-\lambda_{g,2})^{2}}\int\left(g(X)_{1}g(X)_{2}\right)^{2}\
dF(X)$ (9)

When we assume that $F$ is distributed as
$\mathcal{N}(0_{2},\text{diag}(1,\gamma))$, we can evaluate the asymptotic
variance of Equation (9) for different values of $\gamma$. The results of this
computation are presented in Figure 2.

(a) $0\leq\gamma\leq 0.5$

(b) $0.55\leq\gamma\leq 0.90$

Figure 2: Asymptotic variance of the second component of the largest
eigenvector for different values of $\gamma$ at bivariate normal distribution
$\mathcal{N}(0_{2},\text{diag}(1,\gamma))$

From the plots in Figure 2, it is clear that Ball performs the poorest of all
radial functions by a large margin. LR and Shell perform somewhat average. The
lowest ASV is evidently obtained by classical PCA. However, Winsor’s ASV is
quite close to it, followed by SSCM. Quad also has a low asymptotic variance.

We now use these asymptotic variances to calculate the asymptotic relative
efficiency in comparison to the classical method as follows:

$\text{Eff}_{g}(F)=\frac{\text{ASV}_{\text{class}}(F)}{\text{ASV}_{g}(F)}$

where $g$ refers to the chosen radial function. The closer this value is to
one, the more efficient is our method based on the corresponding radial
function $g$. This yields Figure 3. It is immediately clear that the highest
efficiency is obtained by Winsor, it seems to converge to 0.7 for higher
values of $\gamma$. SSCM is inferior, but still achieves values around 0.5.
This is interesting, since it suggests that we can improve on the efficiency
of the popular SSCM without sacrificing robustness (in terms of having a
bounded influence function). Next we observe Quad, attaining values around
0.3. Shell, LR and Ball have low efficiencies at bivariate normal
distributions. Overall Winsor clearly outperforms the other radial functions
in terms of relative efficiency.

The high variance of Ball translates into the lowest efficiency. This does not
come as a surprise, as Ball bears similarities to the affine equivariant
robust covariance estimators which use only half of the data such as the
minimum covariance determinant and minimum volume ellipsoid estimators
(Rousseeuw, 1984). These are known to have fairly low efficiencies on Gaussian
data but strong robustness properties (Davies, 1992; Butler et al., 1993;
Croux and Haesbroeck, 1999; Cator and Lopuhaä, 2012). Ball is also similar to
the BACON algorithm (Billor et al., 2000), which also has strong robustness
but weak efficiency properties.

(a) $0\leq\gamma\leq 0.5$

(b) $0.55\leq\gamma\leq 0.90$

Figure 3: Asymptotic relative efficiencies of the second component of the
largest eigenvector for different values of $\gamma$ at bivariate normal
distribution $\mathcal{N}(0_{2},\text{diag}(1,\gamma))$

###  3.4 Influence functions of the eigenvalues

In addition to the eigenvectors, we also study the influence functions of the
eigenvalues obtained by GSPCA. Denoting these functionals by $\lambda_{g,j}$
for $j=1,\ldots,p$, we could be interested in

$\text{IF}(x,\lambda_{g,j},F)=\lim_{\varepsilon\to
0}\frac{\lambda_{g,j}(F_{\varepsilon,x})-\lambda_{g,j}(F)}{\varepsilon}.$

The problem with using the expression above however, is that the eigenvalue
functionals are not Fisher consistent for elliptical distributions unlike
their eigenvector counterparts. Otherwise stated, in general we have
$\lambda_{g,j}(F)\neq\lambda_{j}(F)$, as well as
$\lambda_{g_{1},j}(F)\neq\lambda_{g_{2},j}(F)$ for two different functions
$g_{1}$ and $g_{2}$. As a result, we cannot directly compare these influence
functions for different functions $g$ which makes them less interesting.

In this section, we will correct the estimated eigenvalues to make them Fisher
consistent. We will then calculate the influence functions of these corrected
eigenvalue estimators, which now can be compared for different radial
functions. Finally, we briefly discuss the covariance matrix obtained by
combining the eigenvectors with the re-estimated eigenvalues.

####  3.4.1 Modification of GSPCA for Fisher consistent eigenvalues

To adapt our method for Fisher consistency, we will use the approach suggested
in Croux and Ruiz-Gazen (2005). The idea proposed there is that we can correct
the eigenvalues by first projecting the data onto the eigenvectors before
calculating the eigenvalues. More specifically we proceed as follows:

  1. 1.

First, we calculate the eigenvectors $v_{g,j}$ through GSPCA as done before.
We know that these are Fisher consistent provided that we have an elliptical
distribution.

  2. 2.

Second, we project the data onto these eigenvectors after which we can use a
robust, equivariant scale estimator $S$ to calculate the eigenvalues:

$\lambda_{S,j}=S^{2}(v_{g,j}^{T}X)$

and for a finite sample:

$\lambda_{S_{n},j}=S_{n}^{2}(v_{g,j}^{T}x_{1},...,v_{g,j}^{T}x_{n})$

Provided that the estimator $S$ itself is a Fisher consistent estimator, this
procedure yields Fisher consistent eigenvalue estimates at elliptically
symmetric distributions (see, e.g., Croux and Ruiz-Gazen (2005)).

  3. 3.

Lastly, should we be interested in the covariance matrix, we can combine the
two steps above to obtain a new estimator $\Sigma_{\text{comb}}$ for the
covariance matrix:

$\displaystyle\Sigma_{\text{comb}}(F)$
$\displaystyle=\sum_{k=1}^{p}\lambda_{S,k}(F)v_{g,k}(F)v_{g,k}(F)^{T}=\sum_{k=1}^{p}\lambda_{k}(F)v_{k}(F)v_{k}(F)^{T}$

Here the last equality holds if $F$ is an elliptically symmetric distribution,
then we have Fisher consistency for $\Sigma$.

In summary, we can correct our method to obtain Fisher consistency for the
eigenvalues and the resulting covariance estimate at elliptically symmetric
distributions. Next, we will take a look at the corresponding influence
functions and asymptotic variances to evaluate the robustness and efficiency
of the adapted GSPCA.

####  3.4.2 Influence functions of the new eigenvalues

In order to obtain the influence function for the corrected eigenvalue
functionals, we need a general property of projections of elliptically
symmetric distributions. Denote the distribution of the random vector
projected on $a$ as $H^{a}$, i.e., if $X\sim H$ then $a^{T}X\sim H^{a}$. We
will make use of following lemma from Croux and Ruiz-Gazen (2005) on the
distribution of this projection $H^{a}$.

######  Lemma 1 (Projected distributions).

Assume $H$ is an elliptically symmetric distribution with location parameter
$\mu$ and covariance matrix $\Sigma$. Then there exists a univariate symmetric
distribution $F_{0}$ such that:

$H^{a}(z)=F_{0}\left(\frac{z-\mu^{T}a}{\sqrt{a^{T}\Sigma a}}\right)$

Lemma 1 ensures that all projections of a random vector with an elliptically
symmetric distribution follow the same symmetric distribution $F_{0}$ after
proper scaling and centering.

We can now obtain an expression for the influence functions of the corrected
eigenvalue functionals. The proof can be found in Section A.4 of the
Supplementary material.

######  Theorem 3 (Influence functions of the new eigenvalues).

For our modified eigenvalues $\lambda_{S,k}$ with $S$ an equivariant scale
functional, we have following expression for the influence function at
elliptically symmetric distributions:

$\text{IF}\left(x,\lambda_{S,k},H\right)=2\lambda_{k}\
\text{IF}\left(\frac{x^{T}v_{k}}{\sqrt{\lambda_{k}}},S,F_{0}\right)$ (10)

This is the same expression as the one obtained by Croux and Ruiz-Gazen (2005)
in the context of robust PCA based on projection-pursuit. In particular, note
that the expression does not depend on the influence function of any of the
eigenvectors.

From the expression of the influence function of the eigenvalues, we readily
obtain the asymptotic variance of the new eigenvalues at elliptically
symmetric distribution:

$\displaystyle\text{ASV}(\lambda_{S,k},H)$
$\displaystyle=\mathbb{E}[\text{IF}(X,\lambda_{S,k},H)^{2}]=\mathbb{E}\left[4\lambda_{k}^{2}\
\text{IF}\left(\frac{v_{k}^{T}X}{\sqrt{\lambda_{k}}},S,F_{0}\right)^{2}\right]=4\lambda_{k}^{2}\
\text{ASV}(S,F_{0})$

Note that the expression for the influence function of the eigenvalues is
proportional to the influence function of the scale functional used. This
means that known optimality properties of equivariant scale functionals carry
over to the estimation of these eigenvalues. In particular, we can leverage
the theory on M-estimation to find optimal M-estimators of scale for
estimating the eigenvalues. This is stated in Corollary 2, which follows from
Theorem 3 and Hampel et al. (1986) (Section 2.5e), and leads us to the use of
the the median absolute deviation (MAD) defined by

$\displaystyle\text{MAD}(X)$ $\displaystyle=c\cdot\text{med}|X-\text{med}(X)|$

Here $c=1/\Phi^{-1}\left(\frac{3}{4}\right)\approx 1.4826$ is a consistency
factor at normal distributions, hence we obtain Fisher consistency for the
estimated eigenvalues.

######  Corollary 2.

If $F_{0}=\mathcal{N}(0,1)$, the median absolute deviation is the most
B-robust M-estimator of the eigenvalues, i.e. it has the lowest gross error
sensitivity among all (well-behaved) Fisher consistent M-estimators of scale.
Similarly, Huber’s M-estimator is the optimal B-robust M-estimator of the
eigenvalues as it has the highest efficiency for a given bound on the gross
error sensitivity.

We visualize the influence function of the new eigenvalue for the multivariate
normal distribution below. Suppose
$H\sim\mathcal{N}(0_{2},\mbox{diag}(1,\rho))$ and thus
$F_{0}\sim\mathcal{N}(0,1)=\Phi$. The influence function of the MAD is given
by (Hampel et al. (1986), p. 107):

$\text{IF}(x,\text{MAD},\Phi)=\frac{\text{sign}\left(|x|-\Phi^{-1}\left(\frac{3}{4}\right)\right)}{4\Phi^{-1}(\frac{3}{4})\phi(\Phi^{-1}\left(\frac{3}{4}\right))}.$

The influence function of the largest eigenvalue $\lambda_{1}=1$ then becomes
$2\ \text{IF}(x_{1},\text{MAD},\Phi)$, which is shown in Figure 4.

Figure 4: Influence function of the largest (new) eigenvalue for $x\in[-3,3]$

####  3.4.3 Influence function of the new covariance estimator

Finally, as an additional result, we compute the influence function of the new
estimate for the covariance matrix $\Sigma_{\text{comb}}$ (see Section A.5 of
the Supplementary Material for the proof).

######  Theorem 4 (Influence function of the new covariance matrix
$\Sigma_{\text{comb}}$).

Let $H$ be an elliptically symmetric distribution with a diagonal covariance
matrix and $S$ an equivariant scale estimator. Then

$\displaystyle\text{IF}(x,\Sigma_{\text{comb}},H)=2\sum_{k=1}^{p}\lambda_{k}\
\text{IF}\left(\frac{v_{k}^{T}x}{\sqrt{\lambda_{k}}},S,F_{0}\right)v_{k}v_{k}^{T}$
$\displaystyle+\sum_{k=1}^{p}\lambda_{k}\sum_{j=1,j\neq
k}^{p}\frac{1}{\lambda_{g,k}-\lambda_{g,j}}(g(x)_{j}g(x)_{k})(v_{j}v_{k}^{T}+v_{k}v_{j}^{T})$

Next we compare the influence functions of the new covariance estimate using
MAD as scale estimator $S$ with the common GSSCM. For this we plot the
influence functions for $H\sim\mathcal{N}(0_{2},\text{diag}(1,\rho))$ and
focus on the off-diagonal element. For that, we obtain:

$\text{IF}(x,\Sigma_{\text{comb}},H)_{1,2}=\frac{1-\rho}{\lambda_{g,1}-\lambda_{g,2}}g(x)_{1}g(x)_{2}$

This corresponds to a rescaling of the influence function of the common GSSCM.
For $\rho=0.5$ and contamination in the direction of $(x,x)^{T}$ this yields
Figure 5.

Figure 5: Influence function of the uncorrected GSSCM $\Sigma$ (left) and the
corrected GSSCM $\Sigma_{\text{comb}}$ (right) for contamination in the
direction of $(x,x)^{T}$: the off-diagonal element

As we can see from the comparison of the raw GSSCM and our proposed
correction, the relative performances among different radial functions remain
relatively stable. However, the influence function of Winsor is now similar to
that of SSCM except in the center where it is lower. This again speaks in
favor of Winsor instead of the popular SSCM. Contamination in the direction of
$(x,0)^{T}$ results in a zero influence function for each radial function
since $g(x)_{2}=0$.

##  4 Simulation study

In this section we conduct a simulation study on synthetically sampled data to
investigate the performance of generalized spherical PCA using the five radial
functions defined in Equations (2) – (6). We will compare them with classical
PCA, spherical PCA based on the SSCM and ROBPCA, the state-of-the-art method
for robust PCA.

The ROBPCA method was published by Hubert et al. (2005). It is based on the
projection pursuit and the minimum covariance determinant estimator. Just as
GSPCA, ROBPCA has the favorable property of being orthogonally equivariant.
For our simulation, we will largely base ourselves on the methodology of the
study conducted in the prior mentioned ROBPCA paper (Hubert et al., 2005).

| $n$ | $p$ | CPCA | ROBPCA | SSCM | Ball | LR | Quad | Shell | Winsor
---|---|---|---|---|---|---|---|---|---|---
Normal | 100 | 4 | .100 | .138 | .137 | .257 | .163 | .132 | .149 | .113
| 50 | 100 | .216 | .245 | .275 | .332 | .272 | .243 | .311 | .227
$t_{5}$ | 100 | 4 | .134 | .138 | .131 | .243 | .163 | .137 | .170 | .113
| 50 | 100 | .303 | .281 | .274 | .395 | .313 | .277 | .378 | .254
$t_{3}$ | 100 | 4 | .185 | .140 | .122 | .222 | .159 | .135 | .174 | .113
| 50 | 100 | .419 | .305 | .272 | .416 | .332 | .294 | .407 | .267
$t_{2}$ | 100 | 4 | .249 | .156 | .129 | .235 | .162 | .142 | .186 | .119
| 50 | 100 | .554 | .354 | .280 | .455 | .358 | .318 | .445 | .285
$t_{1}$ | 100 | 4 | .403 | .183 | .120 | .236 | .173 | .155 | .204 |  .117
(Cauchy) | 50 | 100 | .710 | .495 | .277 | .494 | .403 | .366 | .497 | .304

Table 2: Maxsub measure at uncontaminated data

###  4.1 Setting

In the simulation, we will look at uncontaminated data ($\varepsilon=0$) and
data with different levels of contamination ($\varepsilon=0.1,\
\varepsilon=0.2$) obtained from the following $p$-variate Gaussian
distribution and the $p$-variate Student’s $t$-distribution with five degrees
of freedom:

$\displaystyle(1-\varepsilon)\ \mathcal{N}_{p}(0,\Sigma)+\varepsilon\
\mathcal{N}_{p}(\tilde{\mu},\tilde{\Sigma})$ $\displaystyle(1-\varepsilon)\
t_{5}(0,\Sigma)+\varepsilon\ t_{5}(\tilde{\mu},\tilde{\Sigma})$

From these distributions, we will repeatedly generate 500 data samples of size
$n$. As the GSSCM scatter estimator is orthogonally equivariant, we only need
to look at diagonal covariance matrices $\Sigma$. More in specific, we
consider the following two situations, just as in the ROBPCA paper (Hubert et
al., 2005):

  1. 1.

Low-dimensional data:
$n=100$, $p=4$, $\Sigma=\text{diag}(8,4,2,1)$, $\tilde{\mu}=f_{1}\cdot
e_{4}=(0,0,0,f_{1})^{T}$, $\tilde{\Sigma}=\Sigma/f_{2}$

  2. 2.

High-dimensional data:
$n=50$, $p=100$, $\Sigma=\text{diag}(17,13.5,8,3,1,0.095,...,0.002,0.001)$,
$\tilde{\mu}=f_{1}\cdot e_{6}$, $\tilde{\Sigma}=\Sigma/f_{2}$

In the low-dimensional case we compute $k=3$ principal components, whereas in
the high-dimensional case we compute $k=5$ components, as we want our
principal component analysis to explain at least 90% of the total variance.
Parameter $f_{1}\in\\{6,8,10,\dots,20\\}$ determines the location shift for
the contaminated data in the direction of the $k+1^{\text{th}}$ principal
component, being $e_{4}$ or $e_{6}$. Parameter $f_{2}\in\\{1,15\\}$ specifies
the concentration of the contaminated data.

###  4.2 Maxsub measure

To evaluate performance in our simulation, we will compute the maxsub measure,
which calculates the maximal angle between the estimated PCA subspace and the
space spanned by $\\{e_{1},e_{2},...,e_{k}\\}$. This can be computed as
follows (Hubert et al., 2005):

$\text{maxsub}=\text{arccos}(\sqrt{\lambda_{k}})/(\pi/2)$

where $\lambda_{k}$ represents the smallest eigenvalue of
$I_{p,k}^{T}V_{p,k}V_{p,k}^{T}I_{p,k}$ with $I_{p,k}=(e_{1}\ e_{2}\dots
e_{k})$ and $V_{p,k}=(v_{1}\ v_{2}\dots v_{k})$. We divide by $\pi/2$ to
standardize the value of the angle. The ideal value of the maxsub measure is
zero, the worst value is one.

####  4.2.1 Uncontaminated data

In a first step, we simulate uncontaminated data ($\varepsilon=0$). As
discussed, we sample from the multivariate normal distribution and the
Student’s $t$-distribution. We however do not limit ourselves here to five
degrees of freedom for the $t$-distribution, we also consider $t_{3},\ t_{2}\
\text{and}\ t_{1}$. The mean of the maxsub measure for 500 samples is shown in
Table 2.

When there is no contamination, we see that for the normal and for the
$t_{5}$-distribution the best values are attained by classical PCA (evidently)
and GSPCA with the Winsor radial function. This result is in line with the
computed efficiencies in Section 3.3, where Winsor attained the highest
efficiency of all radial functions. ROBPCA, SPCA and Quad perform almost as
good. The worst results are those of GSPCA with the Ball function, making this
method least efficient at no contamination. Shell also performs poorly,
especially in the high dimensional case. The low efficiency at no
contamination for Ball and Shell is due to the fact that a lot of data points
are given a zero weight, see Equation (4) and (5).

For the $t_{1}$-, $t_{2}$\- and $t_{3}$-distribution the results are
different. These distributions are heavy tailed, meaning that some sampled
data points are located far away from the sample mean. These points at the
tails typically have a significant effect on the classical covariance matrix,
which is clear from Table 2. As the degrees of freedom in the $t$-distribution
decrease, making it more heavy tailed, CPCA attains high values for the maxsub
measure, meaning that it fails. GSPCA with Winsor and SPCA on the other hand,
are not affected by the heavy tailed distribution. Hence in this scenario, the
use of Winsor or the SSCM is suggested over the use of CPCA, even if there is
no contamination.

|
---|---

|
---|---

Figure 6: Maxsub measure for low-dimensional, multivariate normal data

####  4.2.2 Contaminated data

In a second step, we look at the mean of the maxsub measure for 500 samples of
data with different levels of contamination ($\varepsilon=0.1$ and
$\varepsilon=0.2$). The results for multivariate normal and $t_{5}$-data, for
the low- and high-dimensional case, are shown in Figures 6 to 9.

Low-dimensional, normal data:
For the low-dimensional, multivariate normal data (Figure 6) we can clearly
observe that CPCA, SPCA and Winsor fail. Their maxsub measure is close to one,
implying that outliers influenced the estimated PCA subspace to the extent
that one of the estimated principal components is orthogonal to
$\text{span}\\{e_{1},\dots,e_{k}\\}$.

The poor result for these radial functions can be explained by their influence
functions in Section 3.2, where we saw that the influence functions of the
eigenvectors of the classical covariance matrix, the SSCM and the GSSCM of
Winsor did not redescend to zero in contrast with the others. The maxsub
measure of Quad also attains some high values, especially for
$\varepsilon=0.2$ and with $f_{1}$ rather small. When $f_{1}$ is larger,
indicating that the outliers are shifted far enough from the regular
observations, Quad performs better. In contrast, we attain very good maxsub
measure results for LR, Shell and Ball, whose values are as good as the values
for ROBPCA. These three radial functions had influence functions that became
zero after a certain value.

High-dimensional, normal data:
In the high-dimensional normal case (Figure 7), Quad performs worse and fails,
just as CPCA, SPCA and Winsor. LR, Shell and Ball are able to distinguish the
outliers when $f_{1}\geq 10$, in which case their results are comparable to
ROBPCA. However, for $f_{1}<10$, only ROBPCA achieves low values.

|
---|---

|
---|---

Figure 7: Maxsub measure for high-dimensional, multivariate normal data

Low-dimensional, $t_{5}$-data:
For $t_{5}$-data in the low-dimensional case (Figure 8), the results are
similar to low-dimensional normal data. The only difference is that for
$\varepsilon=0.2$ and $f_{1}=6$, all methods perform poorly, except for Ball
who still attains a low value for the maxsub measure. So in this situation,
when the outliers are very close to the regular observations, Ball is the only
one who can distinguish them.

|
---|---

|
---|---

Figure 8: Maxsub measure for low-dimensional, multivariate $t_{5}$-data

High-dimensional, $t_{5}$-data:
Lastly, for the high-dimensional multivariate $t_{5}$-data (Figure 9), we
observe that once again CPCA, SPCA, Winsor and Quad fail. However, Ball and
ROBPCA perform well once $f_{1}$ is greater than or equal to 10. LR and Shell
attain higher values, they only achieve good results when $f_{1}$ exceeds 12
or 14. So for the high-dimensional $t_{5}$ case, the performance of the
methods depends on whether the outliers are separated enough from the regular
observations.
In summary, we can conclude that for the multivariate normal data GSPCA with
radial functions LR, Shell and Ball performs comparable to ROBPCA. For the
$t_{5}$-data, Ball achieves the same results as ROBPCA. LR and Shell perform
somewhat worse, but can still be considered as good alternatives. The fact
that GSPCA with the Ball radial function achieves such great results is
somewhat surprising, since from its definition in Equation (4) it follows that
half of the data points get weight zero and therefore do not contribute
anymore. However, this cutoff makes Ball very robust, as can be seen by the
influence functions of its eigenvectors in Section 3.2, Ball’s influence
function is the first one that becomes zero among all radial functions.

|
---|---

|
---|---

Figure 9: Maxsub measure for high-dimensional, multivariate $t_{5}$-data

###  4.3 Computational time

One of the biggest advantages of our new PCA method is its computational
speed. It is considered to be as fast as classical PCA, while ROBPCA is a
significantly slower method. To illustrate this, we will compare the
computational time of CPCA and ROBPCA with GSPCA using the radial functions
Ball and LR.

We will use the same setting as previously and sample multivariate normal
data, from the low- and high-dimensional case, with the mean vector and
covariance matrix as specified in Section 4.1. We set $\varepsilon=0$ (no
contamination) and measure the total computational time for 100 runs. The
computations were done on a 2.50 GHz core i5 processor ($7^{th}$ gen) and
results are shown in Table 3 and 4.

From both tables it is clear that GSPCA outperforms ROBPCA when we consider
the computational time. The time required for GSPCA is comparable to that of
classical PCA for both the low- and high-dimensional situation. However, the
time needed for ROBPCA increases significantly when the sample size $n$
becomes larger. We conclude that a lot of time can be saved by using GSPCA
instead of ROBPCA.

$n$ | 50 | 100 | 150 | 200 | 250 | 300 | 350 | 400 | 450 | 500
---|---|---|---|---|---|---|---|---|---|---
CPCA | 0.38 | 0.25 | 0.36 | 0.32 | 0.33 | 0.35 | 0.34 | 0.37 | 0.47 | 0.41
ROBPCA | 1.66 | 4.74 | 13.18 | 28.61 | 54.52 | 92.37 | 146.24 |  |  |
GSPCA Ball | 0.23 | 0.23 | 0.38 | 0.37 | 0.41 | 0.45 | 0.52 | 0.58 | 0.62 | 0.71
GSPCA LR | 0.25 | 0.22 | 0.5 | 0.47 | 0.46 | 0.6 | 0.55 | 0.58 | 0.61 | 0.73
Table 3: Total computational time for 100 runs (in seconds) for low-dimensional, multivariate normal data ($p=4$, $k=3$) for various values of $n$ $n$ | 50 | 100 | 150 | 200 | 250 | 300 | 350 | 400 | 450 | 500
---|---|---|---|---|---|---|---|---|---|---
CPCA | 0.89 | 1.34 | 1.88 | 2.28 | 2.86 | 2.93 | 3.58 | 3.88 | 4.18 | 4.6
ROBPCA | 2.78 | 9.34 | 23.78 | 50.18 | 95.04 | 154.86 | 240.79 |  |  |
GSPCA Ball | 1.14 | 1.4 | 1.58 | 2 | 2.19 | 2.67 | 3.11 | 3.45 | 4.11 | 4.35
GSPCA LR | 1.11 | 1.38 | 1.65 | 2 | 2.3 | 2.56 | 3.41 | 3.66 | 3.76 | 4.56
Table 4: Total computational time for 100 runs (in seconds) for high-
dimensional, multivariate normal data ($p=100$, $k=5$) for various values of
$n$

##  5 Data Examples

In this section we illustrate the performance of GSPCA on two real data sets:
the Top Gear car data set and a surveillance video of a beach, previously
studied in a different context in Rousseeuw et al. (2018).

###  5.1 Top Gear data

The first data set comprises the Top Gear car data set from the R package
RobustHD, which has been studied frequently in the context of robustness and
PCA. We consider the numeric variables Price, Displacement, BHP, Torque,
Acceleration, TopSpeed, MPG, Weight, Length, Width and Height and remove
incomplete observations, retaining 245 observations. Next we scale the data
and perform CPCA and GSPCA combined with the LR radial function, retaining 3
principal components to explain at least 85% of the total variance. To
illustrate the robustness of GSPCA, we use the diagnostic plot defined by
Hubert et al. (2005) to classify the flagged outliers, resulting in Figure 10.

We observe that GSPCA identifies many bad leverage points, outlying in score
distance and orthogonal distance, while CPCA only flags 2 of them as bad
leverage points. This suggests that CPCA was heavily influenced by the bad
leverage points identified by GSPCA. We also ran the state-of-the-art method
ROBPCA on the data, its results were very similar to GSPCA.

(a)

(b)

Figure 10: Diagnostic plots for CPCA and GSPCA on the Top Gear data

###  5.2 Video data

|  |  |
---|---|---|---
Figure 11: Video frames 108, 487, 491 and 564

The second data set is comprised of a surveillance video of a beach,
consisting of 633 frames of 160×128 pixels in the RGB color model. The data
was previously studied in the context of outlier detection for functional data
analysis (Rousseeuw et al., 2018) and originates from Li et al. (2004). In the
video we see a beach scenery, where in frame 483 a man comes into view, in
frame 489 he disappears behind a tree and in frames 493 to 633 he remains in
view, see for example Figure 11. Hence from a robustness perspective, we
consider frames 489 to 633 as outliers, while frames 1 to 488 are considered
normal observations.

|  |
---|---|---

|  |
---|---|---

Figure 12: Diagnostic plots for CPCA and GSPCA on the video data in RGB color

To test GSPCA, we run our new method using the LR radial function and the
classical method on the high dimensional, video data. As the video is color
coded in the RGB model, we have 3 data sets ($\mathbf{X1}$, $\mathbf{X2}$ and
$\mathbf{X3}$) of 633 observations with 20480 variables. We select the number
of principal components to explain at least 85% of the variance. The resulting
diagnostic plots of the principal component analyses are shown in Figure 12
for the 3 colors.

Figure 12 shows that GSPCA distinguishes the normal observations (no man,
black points) from the outliers (man in view, colored points). The method
flags all outliers as orthogonal outliers or bad leverage points. It also
separates the frames where the man is behind the tree (green). For CPCA, the
distinction between the outliers and the normal observations is less clear.
Many of the frames containing the man are flagged as good leverage points,
indicating that CPCA was strongly influenced by these points. CPCA also
doesn’t separate the frames where the man is behind the tree.

Also worth noting, GSPCA distinguishes one point with a large score distance
of 20. This observation corresponds to the first frame which is a bit lighter
compared to the others, probably due to starting the video recording. CPCA
also separates this point, but as a point with a large orthogonal distance
around 1000.

To illustrate the effect of the outliers on the analysis, we calculate the
predicted values $\mathbf{\hat{X}}_{(k)}$ of the principal component analysis:

$\mathbf{\hat{X}}_{(k)}=(\mathbf{X}_{n,p}-1_{n}T(\mathbf{X})^{T})V_{p,k}V_{p,k}^{T}+1_{n}T(\mathbf{X})^{T}$

with $k$ the number of principal components to explain 85% of the variance.
From these, we compute the residuals from CPCA and GSPCA to detect outlying
pixels and to study the difference between the original frames and the
predictions:

$r=\mathbf{X}-\mathbf{\hat{X}}_{(k)}$

We also standardize the residuals per frame using the mean and standard
deviation for CPCA, and the median and MAD for GSPCA. In Figure 13 we show
them for frames 108 (no man), 487 (man left), 491 (man behind tree) and 564
(man right). Here all scaled residuals below outlier cutoff 2.5 are given the
same color.

From the residual plot we observe that CPCA only has few mild residuals,
consisting mostly of noise around the man. In GSPCA, on the contrary, the
entire man has very high residuals and is clearly detected by the analysis.
The waves of the sea also have significant residuals as they move throughout
the video.

---
Figure 13: Scaled residuals of data set $\mathbf{X2}$ for frames 108, 487, 491
and 564

##  6 Conclusion

In this paper we developed a new robust version of PCA by calculating the
loading vectors on the generalized spatial sign covariance matrix (GSSCM)
instead of on the classical covariance matrix. The GSSCM was subsequently
combined with different radial functions, five were discussed in this work:
Winsor, Quad, LR, Ball and Shell. This resulted in a new method we named
generalized spherical principal component analysis (GSPCA).

To evaluate GSPCA, three theoretical properties were studied. First, we
calculated a breakdown value of $\left\lfloor{(n+1)/2}\right\rfloor/n$,
meaning that GSPCA can resist upto 50% contamination in a data set. Second,
analytical forms for the influence functions of the eigenvectors of the GSSCM
were derived and studied for the bivariate normal distribution. We
demonstrated that all radial functions except Winsor, had bounded influence
functions which redescended to zero, implying robustness for GSPCA. Last, the
asymptotic relative efficiencies at bivariate normal distribution were
computed, it was shown that GSPCA combined with Winsor is more efficient than
SPCA. Ball, LR and Shell on the contrary had lower efficiencies.

To further examine GSPCA, a simulation study was conducted where the maxsub
measure was analyzed for the different radial functions. First we studied data
sets without contamination to examine efficiency. For the normal distribution
and the Student’s $t$-distribution with five degrees of freedom, classical PCA
and Winsor performed the best. However, for distributions with heavier tails
such as the Student’s $t$-distribution with one or two degrees of freedom,
classical PCA failed. GSPCA with the Winsor radial function, on the contrary,
kept achieving great results.

Subsequently, contamination was introduced to the data sets to investigate
robustness. In the presence of this contamination, Winsor and Quad performed
bad. The other three radial functions LR, Shell and Ball achieved good
results. Their outcomes were comparable to those of ROBPCA, the state-of-the-
art method. However, GSPCA requires significantly less computational time.

Thereafter, our method was tested on real life data. The first data set was
the well-studied, moderately sized Top Gear data set. The second data set
consisted of high-dimensional video data. For both examples, GSPCA
demonstrated favorable robustness properties, whereas CPCA was clearly
influenced by outliers.

We also proposed an adjustment to the GSSCM to achieve Fisher consistency for
the eigenvalues and the covariance estimate, this was done using projections.

In summary, we conclude that the newly developed GSPCA method can be
considered a very interesting alternative for robust PCA based on the findings
stated above. It is however important to combine GSPCA with the appropriate
radial function. When aiming for efficiency, using GSPCA with the Winsor
radial function can be a good alternative to classical PCA. When robustness is
the main priority, GSPCA combined with the LR or Ball radial function are the
most reliable.

## Acknowledgments

First author is supported by Fonds Wetenschappelijk onderzoek - Vlaanderen
(FWO) as a PhD fellow Fundamental Research (PhD fellowship 11K5523N).

## References

  * Billor et al. (2000) N. Billor, A. S. Hadi, and P. F. Velleman.  Bacon: blocked adaptive computationally efficient outlier nominators.  _Computational statistics & data analysis_, 34(3):279–298, 2000.  doi: https://doi.org/10.1016/S0167-9473(99)00101-2.
  * Boente (1987) G. Boente.  Asymptotic theory for robust principal components.  _Journal of Multivariate Analysis_ , 21(1):67–78, 1987.  doi: https://doi.org/10.1016/0047-259X(87)90099-6.
  * Boente et al. (2019) G. Boente, D. Rodriguez, and M. Sued.  The spatial sign covariance operator: Asymptotic results and applications.  _Journal of Multivariate Analysis_ , 170:115–128, 2019\.  doi: https://doi.org/10.1016/j.jmva.2018.10.002.
  * Butler et al. (1993) R. Butler, P. Davies, and M. Jhun.  Asymptotics for the minimum covariance determinant estimator.  _The Annals of Statistics_ , pages 1385–1400, 1993.  doi: https://doi.org/10.1214/aos/1176349264.
  * Campbell (1980) N. A. Campbell.  Robust procedures in multivariate analysis i: Robust covariance estimation.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 29(3):231–237, 1980.  doi: https://doi.org/10.2307/2346896.
  * Cator and Lopuhaä (2012) E. A. Cator and H. P. Lopuhaä.  Central limit theorem and influence function for the mcd estimators at general multivariate distributions.  _Bernoulli_ , 18(2):520–551, 2012.  doi: https://www.jstor.org/stable/23238609.
  * Croux and Haesbroeck (1999) C. Croux and G. Haesbroeck.  Influence function and efficiency of the minimum covariance determinant scatter matrix estimator.  _Journal of Multivariate Analysis_ , 71(2):161–190, 1999.  doi: https://doi.org/10.1006/jmva.1999.1839.
  * Croux and Haesbroeck (2000) C. Croux and G. Haesbroeck.  Principal component analysis based on robust estimators of the covariance or correlation matrix: Influence functions and efficiencies.  _Biometrika_ , 87(3):603–618, 09 2000.  doi: https://doi.org/10.1093/biomet/87.3.603.
  * Croux and Ruiz-Gazen (2005) C. Croux and A. Ruiz-Gazen.  High breakdown estimators for principal components: the projection-pursuit approach revisited.  _Journal of Multivariate Analysis_ , 95:206–226, 2005.  doi: https://doi.org/10.1016/j.jmva.2004.08.002.
  * Croux et al. (2002) C. Croux, E. Ollila, and H. Oja.  Sign and rank covariance matrices: statistical properties and application to principal components analysis.  In _Statistical data analysis based on the L1-norm and related methods_ , pages 257–269. Springer, 2002.  doi: 10.1007/978-3-0348-8201-9˙22.
  * Croux et al. (2007) C. Croux, P. Filzmoser, and M. R. Oliveira.  Algorithms for projection–pursuit robust principal component analysis.  _Chemometrics and Intelligent Laboratory Systems_ , 87(2):218–225, 2007.  doi: https://doi.org/10.1016/j.chemolab.2007.01.004.
  * Croux et al. (2010) C. Croux, C. Dehon, and A. Yadine.  The k-step spatial sign covariance matrix.  _Advances in Data Analysis and Classification_ , 4:137–150, 2010.  doi: https://doi.org/10.1007/s11634-010-0062-7.
  * Davies (1992) L. Davies.  The asymptotics of rousseeuw’s minimum volume ellipsoid estimator.  _The Annals of Statistics_ , 20(4):1828–1843, 1992.  doi: https://doi.org/10.1214/aos/1176348891.
  * Davies (1987) P. L. Davies.  Asymptotic behaviour of s-estimates of multivariate location parameters and dispersion matrices.  _The Annals of Statistics_ , pages 1269–1292, 1987.  doi: https://www.jstor.org/stable/2241828.
  * Dürre et al. (2014) A. Dürre, D. Vogel, and D. E. Tyler.  The spatial sign covariance matrix with unknown location.  _Journal of Multivariate Analysis_ , 130:107–117, 2014\.  doi: https://doi.org/10.1016/j.jmva.2014.05.004.
  * Dürre et al. (2016) A. Dürre, D. E. Tyler, and D. Vogel.  On the eigenvalues of the spatial sign covariance matrix in more than two dimensions.  _Statistics & Probability Letters_, 111:80–85, 2016.  doi: https://doi.org/10.1016/j.spl.2016.01.009.
  * Hampel et al. (1986) F. R. Hampel, E. M. Ronchetti, P. J. Rousseeuw, and W. A. Stahel.  _Robust Statistics: The Approach Based on Influence Functions_.  John Wiley & Sons, 1986.
  * Hubert et al. (2005) M. Hubert, P. Rousseeuw, and K. Vanden Branden.  Robpca: A new approach to robust principal component analysis.  _Technometrics_ , 47(1):64–79, 02 2005.  doi: 10.1198/004017004000000563.
  * Li and Chen (1985) G. Li and Z. Chen.  Projection-pursuit approach to robust dispersion matrices and principal components: primary theory and monte carlo.  _Journal of the American Statistical Association_ , 80(391):759–766, 1985.  doi: https://doi.org/10.2307/2288497.
  * Li et al. (2004) L. Li, W. Huang, I. Y.-H. Gu, and Q. Tian.  Statistical modeling of complex backgrounds for foreground object detection.  _IEEE Transactions on Image Processing_ , 13(11):1459–1472, 2004.  doi: 10.1109/TIP.2004.836169.
  * Locantore et al. (1999) N. Locantore, J. Marron, D. Simpson, N. Tripoli, J. Zhang, and K. Cohen.  Robust principal component analysis for functional data.  _Test_ , 8(1):1–73, 06 1999.  doi: https://doi.org/10.1007/BF02595862.
  * Magyar and Tyler (2014) A. F. Magyar and D. E. Tyler.  The asymptotic inadmissibility of the spatial sign covariance matrix for elliptically symmetric distributions.  _Biometrika_ , 101(3):673–688, 2014.  doi: https://doi.org/10.1093/biomet/asu020.
  * Marden (1999) J. I. Marden.  Some robust estimates of principal components.  _Statistics & probability letters_, 43(4):349–359, 1999.  doi: https://doi.org/10.1016/S0167-7152(98)00272-7.
  * Raymaekers and Rousseeuw (2019) J. Raymaekers and P. Rousseeuw.  A generalized spatial sign covariance matrix.  _Journal of Multivariate Analysis_ , 171:94–111, 2019.  doi: https://doi.org/10.1016/j.jmva.2018.11.010.
  * Rousseeuw and Yohai (1984) P. Rousseeuw and V. Yohai.  Robust regression by means of s-estimators.  In _Robust and nonlinear time series analysis_ , pages 256–272. Springer, 1984.  doi: https://doi.org/10.1007/978-1-4615-7821-5˙15.
  * Rousseeuw (1984) P. J. Rousseeuw.  Least median of squares regression.  _Journal of the American statistical association_ , 79(388):871–880, 1984.  doi: –10.1080/01621459.1984.10477105˝.
  * Rousseeuw et al. (2018) P. J. Rousseeuw, J. Raymaekers, and M. Hubert.  A measure of directional outlyingness with applications to image data and video.  _Journal of Computational and Graphical Statistics_ , 27(2):345–359, 2018.  doi: 10.1080/10618600.2017.1366912.
  * Sibson (1979) R. Sibson.  Studies in the robustness of multidimensional scaling: Perturbational analysis of classical scaling.  _Journal of the royal statistical society series b-methodological_ , 41(2):217–229, 1979.  doi: https://doi.org/10.1111/j.2517-6161.1979.tb01076.x.
  * Taskinen et al. (2012) S. Taskinen, I. Koch, and H. Oja.  Robustifying principal component analysis with spatial sign vectors.  _Statistics & Probability Letters_, 82(4):765–774, 2012.  doi: https://doi.org/10.1016/j.spl.2012.01.001.
  * Visuri et al. (2001) S. Visuri, H. Oja, and V. Koivunen.  Subspace-based direction-of-arrival estimation using nonparametric statistics.  _IEEE Transactions on Signal Processing_ , 49(9):2060–2073, 2001.  doi: 10.1109/78.942634.

##  Appendix A

###  A.1 Proof of Theorem 1

We follow the train of thought of the proof in Raymaekers and Rousseeuw
(2019):

###### Proof.

Part 1: $\varepsilon\geq\left\lfloor{(n+1)/2}\right\rfloor/n$:

Given an $m<\left\lfloor{(n+1)/2}\right\rfloor$, change $m$ of the
observations in $\mathbf{X}$ obtaining the contaminated data set
$\mathbf{X_{m}^{*}}$. As location estimates we have $T(\mathbf{X})$ and
$T(\mathbf{X_{m}^{*}})$ and we define $c_{1}=\max_{i}\left\lVert
x_{i}-T(\mathbf{X})\right\rVert<\infty$. Since $m/n$ is smaller than the
breakdown value of the location estimator, we have that there exist a constant
$c_{2}$ such that $\left\lVert
T(\mathbf{X})-T(\mathbf{X_{m}^{*}})\right\rVert\leq c_{2}<\infty$. Hence by
the triangle inequality, we get $d_{i}^{*}\coloneqq\left\lVert
x_{i}-T(\mathbf{X_{m}^{*}})\right\rVert\leq c_{1}+c_{2}<\infty$. Therefore we
get $\text{med}_{i}(d_{i}^{*})\leq c_{1}+c_{2}$ and hence
$\text{med}_{i}(d_{i}^{*})+1.4826\cdot\text{MAD}_{i}(d_{i}^{*})\leq
2.4826\cdot\text{med}_{i}(d_{i}^{*})\leq 2.4826\cdot(c_{1}+c_{2})$. Using
condition 3, this yields $\left\lVert g(t)\right\rVert\leq
2.4826\cdot(c_{1}+c_{2})$. Now we compute:

$\displaystyle\lambda$ ${}_{\text{max}}=\underset{\left\lVert
u\right\rVert=1}{\sup}u^{T}S_{g}(\mathbf{X_{m}^{*}})u$
$\displaystyle=\underset{\left\lVert
u\right\rVert=1}{\sup}\frac{1}{n}\sum_{i=1}^{n}u^{T}g(x_{i}^{*}-T(\mathbf{X_{m}^{*}}))g(x_{i}^{*}-T(\mathbf{X_{m}^{*}}))^{T}u$
$\displaystyle=\underset{\left\lVert
u\right\rVert=1}{\sup}\frac{1}{n}\sum_{i=1}^{n}\left[u^{T}g(x_{i}^{*}-T(\mathbf{X_{m}^{*}}))\right]^{2}$
$\displaystyle\leq\underset{\left\lVert
u\right\rVert=1}{\sup}\frac{1}{n}\sum_{i=1}^{n}\left\lVert
u\right\rVert^{2}\left\lVert
g(x_{i}^{*}-T(\mathbf{X_{m}^{*}}))\right\rVert^{2}$
$\displaystyle\leq(2.4826\cdot(c_{1}+c_{2}))^{2}<\infty$

Hence we have shown that the largest eigenvalue of $S_{g}(\mathbf{X_{m}^{*}})$
is bounded.

Part 2: $\varepsilon\leq\left\lfloor{(n+1)/2}\right\rfloor/n$:

Given an $m>\left\lfloor{(n+1)/2}\right\rfloor$, replace the last $m$
observations from $\mathbf{X}$ yielding the contaminated data set
$\mathbf{X_{m}^{*}}=\\{x_{1},...,x_{n-m},x_{n-m+1}^{*},...,x_{n}^{*}\\}$.
WLOG, by using location equivariance, we can assume that the mean of
$x_{1},...,x_{n-m}$ is zero. For the other data points,
$j\in\\{n-m+1,...,n\\}$, we put $x_{j}^{*}=\lambda a_{j}$, with $a_{j}$ such
that $\min_{i\in\\{n-m+1,...,n\\}}\left\lVert a_{j}-a_{i}\right\rVert\geq 1$
and such that for all $\lambda>1$:
$\min_{i\in\\{1,...,n-m\\}}\left\lVert\lambda
a_{j}-x_{i}\right\rVert\geq\lambda$. We can do this by placing the $a_{j}$’s
outside of the convex hull of $\mathbf{X}$, far enough apart from each other.

Further, we consider an increasing sequence $(\lambda_{k})_{k}>1$. Then for
every $\lambda_{k}$ the set $\\{x_{n-m+1}^{*},...,x_{n}^{*}\\}$ must contain
at least one point such that $\xi(\left\lVert
x_{i}^{*}-T(\mathbf{X_{m}^{*}})\right\rVert)=1$ by condition 2, say
$x_{b}^{*}$. The set $\mathbf{X_{m}^{*}}$ contains other points with weight 1,
take one arbitrarily and call it $x_{c}^{*}$. By the previous paragraph, we
then have $\left\lVert x_{b}^{*}-x_{c}^{*}\right\rVert\geq\lambda$ and hence
$\left\lVert x_{b}^{*}-T(\mathbf{X_{m}^{*}})\right\rVert+\left\lVert
x_{c}^{*}-T(\mathbf{X_{m}^{*}})\right\rVert\geq\lambda$ and further
$\left\lVert x_{b}^{*}-T(\mathbf{X_{m}^{*}})\right\rVert^{2}+\left\lVert
x_{c}^{*}-T(\mathbf{X_{m}^{*}})\right\rVert^{2}\geq\lambda^{2}/2$. We can then
compute:

$\displaystyle\sum$
${}_{j=1}^{p}\lambda_{j}(S_{g}(\mathbf{X_{m}^{*}}))=\text{trace}(S_{g}(\mathbf{X_{m}^{*}}))$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\text{trace}\left[g(x_{i}^{*}-T(\mathbf{X_{m}^{*}}))g(x_{i}^{*}-T(\mathbf{X_{m}^{*}}))^{T}\right]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left\lVert
g(x_{i}^{*}-T(\mathbf{X_{m}^{*}}))\right\rVert^{2}$
$\displaystyle\geq\frac{1}{n}\left(\left\lVert
x_{b}^{*}-T(\mathbf{X_{m}^{*}})\right\rVert^{2}+\left\lVert
x_{c}*-T(\mathbf{X_{m}^{*}})\right\rVert^{2}\right)$
$\displaystyle\geq\lambda^{2}/(2n)$

For an unbounded increasing sequence of $\lambda$’s this becomes arbitrarily
large and hence $\varepsilon\leq\left\lfloor{(n+1)/2}\right\rfloor/n$. ∎

###  A.2 Proof of Theorem 2

To proof the theorem, we use the following expression for the influence
function of the GSSCM found in the paper of Raymaekers and Rousseeuw (2019):

$\displaystyle\text{IF}(x,\Sigma_{g}$
$\displaystyle,F)=g(x)g(x)^{T}-\Sigma_{g}(F)$
$\displaystyle+\frac{\partial}{\partial\varepsilon}\int
g_{\varepsilon}(X)g_{\varepsilon}(X)^{T}dF(X)\Big{|}_{\varepsilon=0}$ (11)

Here the last term can be expanded as:

$\displaystyle\frac{\partial}{\partial\varepsilon}\int
g_{\varepsilon}(X)g_{\varepsilon}(X)^{T}dF(X)\Big{|}_{\varepsilon=0}=$
$\displaystyle\int\Big{(}\left\\{\frac{\partial}{\partial\varepsilon}g_{\varepsilon}(X)\Big{|}_{\varepsilon=0}\right\\}g(X)^{T}+g(X)\left\\{\frac{\partial}{\partial\varepsilon}g_{\varepsilon}(X)^{T}\Big{|}_{\varepsilon=0}\right\\}\Big{)}dF(X)$

To ease the notational burden, we write $dge(x)$ for the derivative to
$\varepsilon$ of $g_{\varepsilon}(x)$ in $\varepsilon=0$, for which we have:

$\displaystyle dge(x)$
$\displaystyle=\frac{\partial}{\partial\varepsilon}g_{\varepsilon}(X)\Big{|}_{\varepsilon=0}$
$\displaystyle=X\frac{\partial}{\partial\varepsilon}\xi_{\varepsilon}(\left\lVert
X\right\rVert)\Big{|}_{\varepsilon=0}$ (12)

Explicit expressions for Equation (12) per radial function can be found in
appendix A.3 of the GSSCM paper (Raymaekers and Rousseeuw, 2019).
Additionally, we make use of following lemma published in Croux and Haesbroeck
(2000) to compute the influence functions of the eigenvectors of the GSSCM:

######  Lemma 2 (Croux & Haesbroeck 2000).

Given $\Sigma:F\to\text{SPD}(p)$ a statistical functional such that
$\text{IF}(x,\Sigma,F)$ exists. Let $v_{g,1},...,v_{g,p}$ and
$\lambda_{g,1},...,\lambda_{g,p}$ be the eigenvectors and eigenvalues of
$\Sigma(F)$. Then the influence function of $v_{g,j}$ is given by:

$\text{IF}(x,v_{g,j},F)=\sum_{k=1,k\neq
j}^{p}\frac{1}{\lambda_{g,j}-\lambda_{g,k}}\left(v_{g,k}^{T}\text{IF}(x,\Sigma,F)v_{g,j}\right)v_{g,k}$

The proof for this lemma follows from Lemma 2.1 in Sibson (1979). We are now
ready to proof the theorem.

###### Proof.

Making use of Lemma 2 and Equation (11) for the influence function of the
GSSCM, we note:

$\displaystyle\text{IF}(x,v_{g,j},F)$ $\displaystyle=\sum_{k=1,k\neq
j}^{p}\frac{1}{\lambda_{g,j}-\lambda_{g,k}}\left(v_{g,k}^{T}\text{IF}(x,\Sigma_{g},F)v_{g,j}\right)v_{g,k}$
$\displaystyle=\sum_{k=1,k\neq
j}^{p}\frac{1}{\lambda_{g,j}-\lambda_{g,k}}\Big{(}v_{g,k}^{T}g(x)g(x)^{T}v_{g,j}-v_{g,k}^{T}\Sigma_{g}(F)v_{g,j}$
$\displaystyle\qquad\quad+v_{g,k}^{T}\frac{\partial}{\partial\varepsilon}\int
g_{\varepsilon}(X)g_{\varepsilon}(X)^{T}dF(X)\Big{|}_{\varepsilon=0}v_{g,j}\Big{)}v_{g,k}$
$\displaystyle=\sum_{k=1,k\neq
j}^{p}\frac{1}{\lambda_{g,j}-\lambda_{g,k}}\Big{(}(v_{g,k}^{T}g(x))(v_{g,j}^{T}g(x))-v_{g,k}^{T}\lambda_{g,j}v_{g,j}$
$\displaystyle\qquad\quad+v_{g,k}^{T}\frac{\partial}{\partial\varepsilon}\int
g_{\varepsilon}(X)g_{\varepsilon}(X)^{T}dF(X)\Big{|}_{\varepsilon=0}v_{g,j}\Big{)}v_{g,k}$
$\displaystyle=\sum_{k=1,k\neq
j}^{p}\frac{1}{\lambda_{g,j}-\lambda_{g,k}}\Bigg{(}(v_{g,k}^{T}g(x))(v_{g,j}^{T}g(x))\
-\ 0$
$\displaystyle\qquad\quad+v_{g,k}^{T}\int\Big{(}\left\\{\frac{\partial}{\partial\varepsilon}g_{\varepsilon}(X)\Big{|}_{\varepsilon=0}\right\\}g(X)^{T}$
$\displaystyle\qquad\quad+g(X)\left\\{\frac{\partial}{\partial\varepsilon}g_{\varepsilon}(X)^{T}\Big{|}_{\varepsilon=0}\right\\}\Big{)}dF(X)\
v_{g,j}\Bigg{)}v_{g,k}$ $\displaystyle=\sum_{k=1,k\neq
j}^{p}\frac{1}{\lambda_{g,j}-\lambda_{g,k}}\Big{(}(v_{g,k}^{T}g(x))(v_{g,j}^{T}g(x))$
$\displaystyle\qquad\quad+v_{g,k}^{T}\int\big{(}dge(X)g(X)^{T}$
$\displaystyle\qquad\quad+g(X)dge(X)^{T}\big{)}dF(X)\ v_{g,j}\Big{)}v_{g,k}$

∎

###  A.3 Proof of Corollary 1

###### Proof.

When we assume that $F$ is a centered elliptically symmetric distribution with
diagonal covariance matrix, we have that $v_{g,j}=v_{j}=e_{j}$. Hence Equation
(2) simplifies to:

$\displaystyle\text{IF}(x,v_{g,j},F)=\sum_{k=1,k\neq
j}^{p}\frac{1}{\lambda_{g,j}-\lambda_{g,k}}\Big{(}g(x)_{k}g(x)_{j}$
$\displaystyle+\int\\{dge(X)_{k}\ g(X)_{j}+g(X)_{k}\
dge(X)_{j}\\}dF(X)\Big{)}v_{k}$ (13)

Next we proof that the above integral in the second term equals zero. For
this, we first know that $g(x)=x\cdot\xi(\left\lVert x\right\rVert)$, hence
$g(x)_{i}=x_{i}\cdot\xi(\left\lVert x\right\rVert)$. Second, if we study
appendix A.3 of the GSSCM paper (Raymaekers and Rousseeuw, 2019), one can find
that for every suggested radial function $dge(x)$ is of the form
$dge(x)=\text{scalar}\cdot x\cdot f(\left\lVert x\right\rVert)$. Therefore, we
have that $dge(x)_{i}=\text{scalar}\cdot x_{i}\cdot f(\left\lVert
x\right\rVert)$. Last, the density function of a centered elliptically
symmetric distribution has the following form:
$f_{X}(x)=\text{det}(\Sigma)^{-1/2}\cdot h(x^{T}\Sigma^{-1}x)$ with $\Sigma$
diagonally assumed here, hence it is an even function in each of its variables
$x_{i}$. Putting all this together, the integral in the second term of
Equation (A.3) is of following form

$\displaystyle\int_{-\infty}^{+\infty}\dots\int_{-\infty}^{+\infty}x_{k}x_{j}\cdot
f(\left\lVert x\right\rVert)\xi(\left\lVert x\right\rVert)\cdot\omega(x)\cdot
dx_{1}dx_{2}\dots dx_{p}$ $\displaystyle\text{for }k\neq j$

where $\omega(x)$ is an even function in each variable. Then one has that this
integral is equal to zero because of symmetry reasons. Hence the second term
in Equation (A.3) becomes zero for centered elliptically symmetric
distributions. ∎

###  A.4 Proof of Theorem 3

Using Lemma 1, we can find an expression for the influence functions of our
new eigenvalues:

###### Proof.

$\displaystyle\lambda_{S,k}(H)=S^{2}(v_{g,k}^{T}X)=S^{2}(H^{v_{g,k}})$
$\displaystyle\implies$
$\displaystyle\text{IF}(x,\lambda_{S,k},H)=\frac{d}{d\varepsilon}\left(S^{2}\left(H_{\varepsilon,x}^{v_{g,k}(H_{\varepsilon,x})}\right)\right)\Big{|}_{\varepsilon=0}$
$\displaystyle=\frac{d}{d\varepsilon}\left(S^{2}(H_{\varepsilon,x}^{v_{k}})\right)\Big{|}_{\varepsilon=0}$
$\displaystyle\qquad+\left(\frac{d}{da}\left(S^{2}(H^{a})\right)\Big{|}_{a=v_{g,k}}\right)^{T}\cdot\text{IF}(x,v_{g,k},H)$

We can use that $S^{2}(H^{a})=a^{T}\Sigma a$, as $S$ is equivariant, which
yields:

$\displaystyle\frac{d}{da}\left(S^{2}(H^{a})\right)\Big{|}_{a=v_{g,k}}$
$\displaystyle=\left(\frac{d}{da}\left(a^{T}\Sigma
a\right)\Big{|}_{a=v_{g,k}}\right)$ $\displaystyle=2\Sigma v_{g,k}=2\Sigma
v_{k}=2\lambda_{k}v_{k}$

We can then continue our calculation:

$\displaystyle\text{IF}(x,\lambda_{S,k},H)$
$\displaystyle=\frac{d}{d\varepsilon}\left(S^{2}(H_{\varepsilon,x}^{v_{k}})\right)\Big{|}_{\varepsilon=0}+2\lambda_{k}v_{k}^{T}\cdot\text{IF}(x,v_{g,k},H)$
$\displaystyle=\text{IF}(v_{k}^{T}x,\text{S}^{2},H^{v_{k}})+2\lambda_{k}v_{k}^{T}\cdot\text{IF}(x,v_{g,k},H)$

When we calculated the influence function of the eigenvector $v_{g,k}$ of
GSSCM, we saw that this function had no component in the direction of $v_{k}$.
Therefore we can remove the second term. Hence we obtain:

$\displaystyle\text{IF}(x,\lambda_{S,k},H)$
$\displaystyle=\text{IF}(v_{k}^{T}x,S^{2},H^{v_{k}})$
$\displaystyle=v_{k}^{T}\Sigma v_{k}\
\text{IF}\left(\frac{v_{k}^{T}x}{\sqrt{v_{k}^{T}\Sigma
v_{k}}},S^{2},F_{0}\right)$ $\displaystyle=2\lambda_{k}\
\text{IF}\left(\frac{v_{k}^{T}x}{\sqrt{\lambda_{k}}},S,F_{0}\right)$

∎

###  A.5 Proof of Theorem 4

###### Proof.

IF $\displaystyle(x,\Sigma_{\text{comb}},H)$
$\displaystyle=\frac{d}{d\varepsilon}\left(\sum_{k=1}^{p}\lambda_{S,k}(H_{\varepsilon,x})v_{g,k}(H_{\varepsilon,x})v_{g,k}(H_{\varepsilon,x})^{T}\right)\Big{|}_{\varepsilon=0}$
$\displaystyle=\sum_{k=1}^{p}\left[\text{IF}(x,\lambda_{S,k},H)v_{k}v_{k}^{T}\right.$
$\displaystyle\left.\qquad+\ \lambda_{k}\ \text{IF}(x,v_{g,k},H)\
v_{k}^{T}+\lambda_{k}v_{k}\ \text{IF}(x,v_{g,k},H)^{T}\right]$

Assuming that the original covariance matrix is diagonal and that $H$ is
elliptically symmetric, we can use Corollary 1 and Equation (10) to obtain:

$\displaystyle\text{IF}($ $\displaystyle x,\Sigma_{\text{comb}},H)$
$\displaystyle=\sum_{k=1}^{p}\big{[}\text{IF}(x,\lambda_{S,k},H)v_{k}v_{k}^{T}+\lambda_{k}\
\text{IF}(x,v_{g,k},H)\ v_{k}^{T}$ $\displaystyle\quad+\lambda_{k}v_{k}\
\text{IF}(x,v_{g,k},H)^{T}\big{]}$ $\displaystyle=\sum_{k=1}^{p}2\lambda_{k}\
\text{IF}\left(\frac{v_{k}^{T}x}{\sqrt{\lambda_{k}}},S,F_{0}\right)v_{k}v_{k}^{T}$
$\displaystyle\quad+\sum_{k=1}^{p}\lambda_{k}\Big{[}\sum_{j=1,j\neq
k}^{p}\frac{1}{\lambda_{g,k}-\lambda_{g,j}}(g(x)_{j}g(x)_{k})v_{j}v_{k}^{T}$
$\displaystyle\quad+\sum_{j=1,j\neq
k}^{p}\frac{1}{\lambda_{g,k}-\lambda_{g,j}}(g(x)_{j}g(x)_{k})v_{k}v_{j}^{T}\Big{]}$
$\displaystyle=2\sum_{k=1}^{p}\lambda_{k}\
\text{IF}\left(\frac{v_{k}^{T}x}{\sqrt{\lambda_{k}}},S,F_{0}\right)v_{k}v_{k}^{T}$
$\displaystyle\quad+\sum_{k=1}^{p}\lambda_{k}\sum_{j=1,j\neq
k}^{p}\frac{1}{\lambda_{g,k}-\lambda_{g,j}}(g(x)_{j}g(x)_{k})(v_{j}v_{k}^{T}+v_{k}v_{j}^{T})$

∎