# Cauchy robust principal component analysis with applications to high-
deimensional data sets

Ayisha Fayomi 1 Yannis Pantazis 2 Michail Tsagris 3
and Andrew T.A. Wood 4

1 Department of Statistics  King Abdulaziz University
Jeddah  Saudi Arabia
<EMAIL_ADDRESS>
2 Institute of Applied & Computational Mathematics
Foundation for Research & Technology - Hellas  Heraklion  Greece
<EMAIL_ADDRESS>
3 Department of Economics  University of Crete
Gallos Campus  Rethymnon  Greece
<EMAIL_ADDRESS>
4 Research School of Finance  Actuarial Studies & Statistics
Australian National University  Canberra  Australia
<EMAIL_ADDRESS>
###### Abstract

Principal component analysis (PCA) is a standard dimensionality reduction
technique used in various research and applied fields. From an algorithmic
point of view, classical PCA can be formulated in terms of operations on a
multivariate Gaussian likelihood. As a consequence of the implied Gaussian
formulation, the principal components are not robust to outliers. In this
paper, we propose a modified formulation, based on the use of a multivariate
Cauchy likelihood instead of the Gaussian likelihood, which has the effect of
robustifying the principal components. We present an algorithm to compute
these robustified principal components. We additionally derive the relevant
influence function of the first component and examine its theoretical
properties. Simulation experiments on high-dimensional datasets demonstrate
that the estimated principal components based on the Cauchy likelihood
outperform or are on par with existing robust PCA techniques.

##  1 Introduction

In the analysis of multivariate data, it is frequently desirable to employ
statistical methods which are insensitive to the presence of outliers in the
sample. To address the problem of outliers, it is important to develop robust
statistical procedures. Most statistical procedures include explicit or
implicit prior assumptions about the distribution of the observations, but
often without taking into account the effect of outliers. The purpose of this
paper is to present a novel robust version of PCA which has some attractive
features.

Principal components analysis (PCA) is considered to be one of the most
important techniques in statistics. However, the classical version of PCA
depends on either a covariance or a correlation matrix, both of which are very
sensitive to outliers. We develop an alternative method to classical PCA,
which is far more robust, by using a multivariate Cauchy likelihood to
construct a robust principal components (PC) procedure. It is an adaptation of
the classic method of PCA obtained by replacing the Gaussian log-likelihood
function by the Cauchy log-likelihood function, in a sense that will be
explained in section 2.2. Although we do not claim that the interpretation of
standard PCA in terms of operations on a Gaussian likelihood is new, see
Bolton and Krzanowski, this fact does not appear to have been exploited in the
development of a robust PCA procedure, as we do in this paper. An important
reason for using the multivariate Cauchy likelihood is that this likelihood
has only one maximum point, but the single most important motivation is that
it leads to a robust procedure.

In the next section we review briefly some of the techniques employed for
estimating parameters and for directing a PCA in ways which are robust against
the presence of outliers. We also present robustness preliminaries that
include some important techniques which are necessary to assess whether the
method used is robust or not. In Section 3 we develop the Cauchy-PCA and
theoretically explore its robustness properties. Finally, in Section 4 we
present the numerical algorithms for creating Cauchy PCs, and also give the
results of a number of very high-dimensional real-data and simulated examples.
Our approach is seen to be competitive with, and often gives superior results
to, that of the projection pursuit algorithm of Croux et al. (2007, 2013).
Finally we conclude the paper in Section 5.

###  1.1 Literature review on robust PCA

It is well known that PCA is an important technique for high-dimensional data
reduction. PCA is based on the sample covariance matrix $\hat{{\bf\Sigma}}$
and it involves searching for a linear combination $y_{j}={\bf u}^{T}{\bf
x}_{j}$ of the ${\bf x}$ components of the vector that maximize the sample
variance of the components of $y$. According to Mardia et al., (1979), the
solution will be given by the equation

$\hat{\bf\Sigma}={\bf U\Lambda U}^{T},$

where ${\bf\Lambda}=\hbox{diag}\\{\lambda_{1},\ldots,\lambda_{p}\\}$ and its
diagonal elements $\lambda_{i}$ are the sample variances, while ${\bf U}$ is
an orthogonal matrix, i.e. ${\bf UU}^{T}={\bf U}^{T}{\bf U}={\bf I}_{p}$,
whose columns ${\bf u}_{i}$ are the corresponding eigenvectors which represent
the linear combinations. [[The principal components are efficiently estimated
in practice via Singular Value Decomposition (SVD) (cite Lanczos for an
efficient algorithm).]]

Classical PCA, unfortunately, is non-robust, since it based on the sample
covariance or sample correlation matrix which are very sensitive to outlying
observations; see section 2. However, this problem has been handled by two
different methods which result in robust versions of PCA by:

i.

replacing the standard covariance or correlation matrix with a robust
estimator; or

ii.

maximising (or minimising) a different objective function to obtain a robust
PCA.

Many different proposes had been developed to carry out robust PCA, such that
using projection pursuit PP, $M-$estimators and so on.

Despite maximum likelihood estimation, perhaps, being considered as the most
important statistical inference method, sometimes this approach can lead to
improper results when the underlying assumptions are not satisfied, for
instance, when data contain outliers or deviate slightly from the supposed
model. A generalization of maximum likelihood estimation proposed by Huber,
(1964) which is called $M$-estimation, aims to produce a robust statistic by
constructing approaches that are resistant to deviations from the underline
assumptions. $M$-estimators were also defined for the multivariate case by
Maronna, (1976).

Campbell, (1980) provided a procedure for robust PCA by examining the
estimates of means and covariances which are less affected by outlier
observations, and by exploring the observations which have a large effect on
the estimates. He replaced the sample covariance sample by an $M-$estimator.
Hubert and Verboven, (2003) introduced a new approach to create robust PCA. It
combines the advantages of two methods, the first one is based on replacing
the covariance or correlation matrix by its robust estimator, while the second
one is based on maximizing the objective function for this robust estimate.

A robust PCA based on the projection pursuit (PP) method was developed by Li
and Chen, (1985), using Huber’s $M$-estimator of dispersion as the projection
index. The objective of PP is to seek projections, of the high-dimensional
data set onto low-dimensional subspaces, that optimise a function of
”interestingness”. The function that should be optimised is called an index or
objective function and its choice depends on a feature that the researcher is
concerned about. This property gives the PP technique a flexibility to handle
many different statistical problems range from clustering to identifying
outliers in a multivariate data set.

Bolton and Krzanowski, (1999) characterized the PC’s for PP in terms of
maximum likelihood under the assumption of normality. PCA can be considered as
a special case of PP as well as many other methods of multivariate analysis.
Li and Chen, (1985) used Huber’s $M$-estimator of dispersion as projective
index to develop a robust PCA based on the PP approach. The sample median was
used as a projective index to develop a robust PCA by Xie et al., (1993). In
their simulation studies, Xie et al., (1993) observed a PCA resistant to
outliers and deviations from the normal distribution. Croux et al., (2007,
2013) also suggested a robust PCA using projection pursuit and we will
contrast our methodology against their algorithm.

##  2 Preliminaries on standard PCA

PCA is an orthogonal linear transformation that projects the data to a new
coordinate system according to the variance of each direction. Given a data
matrix ${\bf X}\in\mathbb{R}^{n\times p}$ with each row correspond to a
sample, the first direction ${\bf u}_{1}$ that maximizes the variance is
defined through

${\bf u}_{1}=\operatorname*{\arg\max}_{||{\bf u}||_{2}=1}||({\bf X}-{\bf
1}_{n}\bar{{\bf x}}^{T}){\bf u}||_{2}^{2},$

where ${\bf 1}_{n}$ is an $n$-dimensional vector whose elements are all set to
1 while $\bar{{\bf x}}=\frac{1}{n}\sum_{i=1}^{n}{\bf x}_{i}$ is the empirical
mean. The process is repeated $k$ times and at each iteration the to-be-
estimated principal direction has to be orthogonal to all previously-computed
principal directions. Thus, the $k$-th direction which has to be orthogonal to
the previous ones is defined by

${\bf u}_{k}=\operatorname*{\arg\max}_{||{\bf u}||_{2}=1}||({\bf X}-{\bf
1}_{n}\bar{{\bf x}}^{T}){\bf u}||_{2}^{2}\ \ \text{subject to}\ \ {\bf
u}_{k}\perp{\bf u}_{j}\ \text{with}\ j=1,...,k-1\ .$

###  2.1 Non-robustness of standard PCA

We will show that the influence function for the largest eigenvalue of the
covariance matrix and the respective eigenvector are unbounded with respect to
the norm of an outlier sample. Suppose that $\boldsymbol{\Sigma}$ is the
covariance matrix of a population with distribution function $F$, i.e.,

${\boldsymbol{\Sigma}}=\int_{\mathbb{R}^{p}}({\bf x}-\boldsymbol{\mu})({\bf
x}-\boldsymbol{\mu})^{T}dF({\bf x}),$ (1)

where $\boldsymbol{\mu}=\int_{\mathbb{R}^{p}}{\bf x}dF({\bf x})$ corresponds
to the mean vector. Assume that the leading eigenvalue of
$\boldsymbol{\Sigma}$ has multiplicity 1, then we denote it by $\lambda$ and
the leading eigenvector by $\hat{{\bf u}}$ (i.e., ${\bf u}_{1}=\hat{{\bf
u}}$).

Let $T$ be an arbitrary functional, $F$ a distribution and ${\bf
z}\in\mathbb{R}^{p}$ an arbitrary point in the relevant sample space. The
influence function is defined as

$IF_{T}({\bf z};F)=\lim_{\epsilon\to
0+}\frac{T((1-\epsilon)F+\epsilon\Delta_{{\bf z}})-T(F)}{\epsilon},$ (2)

where $\Delta_{{\bf z}}$ is a unit point mass located at ${\bf z}$.

A robust estimator for $T$ means that the influence function is bounded with
respect to the norm of the outlier $\bf z$.

######  Proposition 2.1.

The influence function for the leading eigenvector of $\boldsymbol{\Sigma}$ is
given by111We use ${\bf A}^{+}$ to denote the Moore-Penrose inverse of a
matrix $\bf A$.

$IF_{\hat{{\bf u}}}({\bf z},F)=-\big{(}({\bf z}-\boldsymbol{\mu})^{T}\hat{{\bf
u}}\big{)}(\boldsymbol{\Sigma}-\lambda{\bf I}_{p})^{+}({\bf
z}-\boldsymbol{\mu}).$ (3)

Similarly, the IF for the largest eigenvalue of ${\boldsymbol{\Sigma}}$ is

$IF_{\lambda}({\bf z},F)=\big{(}({\bf z}-\boldsymbol{\mu})^{T}\hat{{\bf
u}}\big{)}^{2}-\lambda.$ (4)

The detailed calculations are presented in Appendix A1. The following result
shows that outliers with unbounded influence function do exist.

######  Corollary 2.1.

Let ${\bf z}=\boldsymbol{\mu}+\gamma\hat{{\bf u}}+\eta{\bf v}$ where ${\bf v}$
is orthogonal to $\hat{{\bf u}}$ and does not belong to the null space of
$\boldsymbol{\Sigma}$ and $\gamma,\eta\neq 0$ then

$\lim_{{\bf z}:\,||{\bf z}||_{2}\rightarrow\infty}||IF_{\hat{{\bf u}}}({\bf
z},F)||_{2}=\infty,$

and similarly for $IF_{\lambda}({\bf z},F)$.

###### Proof.

Direct substitution of ${\bf z}$ into the influence function gives:

$IF_{\hat{{\bf u}}}({\bf z},F)=-((\gamma\hat{{\bf u}}+\eta{\bf
v})^{T}\hat{{\bf u}})(\boldsymbol{\Sigma}-\lambda{\bf
I}_{p})^{+}(\gamma\hat{{\bf u}}+\eta{\bf
v})=-\gamma\eta(\boldsymbol{\Sigma}-\lambda{\bf I}_{p})^{+}{\bf v}.$

Since ${\bf v}$ does not belong to the null space of $\boldsymbol{\Sigma}$, it
holds that $(\boldsymbol{\Sigma}-\lambda{\bf I}_{p})^{+}{\bf v}\neq{\bf 0}$
thus $||(\boldsymbol{\Sigma}-\lambda{\bf I}_{p})^{+}{\bf v}||_{2}=c\neq 0$.
Hence,

$||IF_{\hat{{\bf u}}}({\bf z},F)||_{2}=|\gamma||\eta|c.$

Given that $||{\bf
z}||_{2}^{2}=\gamma^{2}+\eta^{2}+||\boldsymbol{\mu}||_{2}^{2}+\gamma\boldsymbol{\mu}^{T}\hat{{\bf
u}}+\eta\boldsymbol{\mu}^{T}{\bf v}$, either sending $|\gamma|\to\infty$ or
$|\eta|\to\infty$ completes the proof.

Similarly,

$IF_{\lambda}({\bf z},F)=\gamma^{2}-\lambda\rightarrow\infty,$

as $|\gamma|\to\infty$. ∎

###  2.2 Generalizations of standard PCA

Standard PCA can be viewed as a special case of a more general optimization
problem. We present two such generalization: the first one leads to projection
pursuit algorithms while the second leads to a maximum likelihood formulation.
Let ${\bf u}$ be a unit vector and define the projection values

$c_{i}({\bf u})={\bf x}^{T}_{i}{\bf u},{\hskip 8.53581pt}i=1,\ldots,n,$

and a function $\Phi:\mathbb{R}^{n}\to\mathbb{R}$ acting on the projected
values. The first generalization of PCA is defined as the maximization of
$\Phi$:

${\bf u}_{1}=\operatorname*{\arg\max}_{||{\bf u}||_{2}=1}\Phi(c_{1}({\bf
u}),...,c_{n}({\bf u}))\ .$

As in the standard PCA, the following principal directions are obtained after
removing the contribution of the current principal component from the data.
When $\Phi$ is the sample variance then we recover the standard PCA.

The second generalization interprets the computation of the principal
component as a maximum likelihood estimation problem. By letting,

$l_{G}(\mu,\sigma^{2}|c_{1},\ldots,c_{n})=-\frac{n}{2}\log{\sigma}^{2}-\frac{n}{2{\sigma}^{2}}\sum_{i=1}^{n}(c_{i}-\mu)^{2}.$
(5)

be the Gaussian log-likelihood, the first principal direction can be obtained
by solving the minimax problem:

$\min_{||{\bf u}||_{2}=1}\max_{\mu,\sigma^{2}}\
l_{G}(\mu,\sigma^{2}|c_{1}({\bf u}),\ldots,c_{n}({\bf u})).$

Indeed, the inner maximization can be solved analytically which leads to the
optimal solution

$\hat{\mu}({\bf u})=\frac{1}{n}\sum_{i=1}^{n}c_{i}({\bf u})=:\bar{c}({\bf u})$

and

${\hat{\sigma}}^{2}({\bf u})=\frac{1}{n}\sum_{i=1}^{n}(c_{i}({\bf
u})-\bar{c}({\bf u}))^{2}.$

Unsurprisingly, the optimal values are the sample mean and the sample
variance. Using the above formulas it is straightforward to show that

$\displaystyle\operatorname*{\arg\min}_{||{\bf u}||_{2}=1}\
l_{G}\big{(}\hat{\mu}({\bf u}),{\hat{\sigma}}^{2}({{\bf u}})|c_{1}({{\bf
u}}),\ldots,c_{n}({{\bf u}})\big{)}=\operatorname*{\arg\max}_{||{\bf
u}||_{2}=1}\ \hat{\sigma}^{2}({{\bf u}})\ .$ (6)

Variations of PCA can be derived by changing the likelihood function and in
the next section we analyze the case of Cauchy distribution.

##  3 Cauchy PCA

The Cauchy log-likelihood function is given by

$l_{C}({\mu},{\sigma}|{c}_{1}({{\bf u}}),\ldots,{c}_{n}({{\bf
u}}))=n\log{\frac{\sigma}{\pi}}-\sum_{i=1}^{n}\log\left\\{{\sigma}^{2}+(c_{i}({\bf
u})-{\mu})^{2}\right\\}.$ (7)

where $\mu$ and $\sigma$ are the two parameters of the Cauchy distribution.
The first Cauchy principal direction is also obtained by solving the minimax
optimization problem:

$\min_{||{\bf u}||_{2}=1}\max_{\mu,\sigma}\ l_{C}(\mu,\sigma^{2}|c_{1}({\bf
u}),\ldots,c_{n}({\bf u})).$ (8)

In contrast to the Gaussian case, the inner maximization cannot be performed
analytically. Therefore an iterative approach needs to be utilized. Here, we
apply the Newton-Raphson method with initial values the median and half the
interquartile range for the location and scale parameters, respectively.
According to Copas, (1975), although the mean of the Cauchy distribution does
not exist and it has infinite variance, the Cauchy log-likelihood function
$l_{C}(\mu,\sigma)$ has a unique maximum likelihood estimate,
$(\hat{\mu},\hat{\sigma})$.

Fixing $\mu$ and $\sigma$, the outer minimization is also non-analytic and a
fixed point iteration is applied to calculate ${\bf u}$. The iteration is
given by

$\hat{{\bf u}}=\frac{\hat{{\bf u}}_{un}}{||\hat{{\bf u}}_{un}||_{2}},$ (9)

where $\hat{{\bf u}}_{un}$ is the unnormalized direction which is obtained
from the differentiation of the Lagrangian function with respect to ${\bf u}$
and it is given by

$\displaystyle\hat{{\bf u}}_{un}=\sum_{i=1}^{n}\frac{({{\bf
x}}_{i}^{T}{\hat{{\bf u}}}-\hat{\mu}){{\bf
x}}_{i}}{\hat{\sigma}^{2}+\left({{\bf x}}_{i}^{T}{\hat{{\bf
u}}}-\hat{\mu}\right)^{2}}\ .$ (10)

Once the first principal direction has been computed, its contribution from
the dataset ${\bf X}$ is removed and the same procedure to estimate the next
principal direction is repeated. This iterative process is repeated $k$ times.
The removal of the contribution makes the principal directions orthogonal to
each other. We summarize the estimation of $k$ Cauchy principal components in
the following pseudo-code (Algorithm 1).

Algorithm 1 Cauchy PCA

for $j=1,...,k$ do

$\bullet$ Initialize ${\hat{{\bf u}}_{un}}$ and normalize $\hat{{\bf
u}}=\hat{{\bf u}}_{un}/||\hat{{\bf u}}_{un}||_{2}$

while not converged do

$\bullet$ Fix $\hat{{\bf u}}$ and set  $c_{i}(\hat{{\bf u}})={\bf
x}_{i}^{T}\hat{{\bf u}},\ \ i=1,...,n.$

$\bullet$ Via Newton-Raphson algorithm find
$(\hat{\mu},\hat{\sigma})=\operatorname*{\arg\max}_{\mu,\sigma}\
l_{C}(\mu,\sigma;c_{1}(\hat{{\bf u}}),\ldots,c_{n}(\hat{{\bf u}})).$

$\bullet$ Fix $(\hat{\mu},\hat{\sigma})$ and using fixed point iteration
(i.e., (10) & (9)) find  $\hat{{\bf u}}=\operatorname*{\arg\min}_{{\bf u}}\
l_{C}(\hat{\mu},\hat{\sigma}|c_{1}({\bf u}),\ldots,c_{n}({\bf
u}))-\lambda(||{\bf u}||_{2}^{2}-1)$

end while

$\bullet$ Set the $j$-th Cauchy principal direction  ${\bf u}_{j}=\hat{{\bf
u}}.$

$\bullet$ Remove the contribution from the dataset  $\displaystyle{\bf X}={\bf
X}({\bf I}_{p}-{\bf u}_{j}{\bf u}^{T}_{j}),$

end for

###  3.1 Robustness of the Leading Cauchy Principal Direction

Let $\boldsymbol{\theta}=\left(\mu,\sigma\right)^{T}$ be the parameter vector
of the Cauchy distribution and consider the infinite-sample normalized Cauchy
log-likelihood function

$l({\bf u}|\boldsymbol{\theta})=\int_{{\bf x}\in\mathbb{R}^{p}}g(c({\bf
u}),\boldsymbol{\theta})\,dF({\bf x}),$ (11)

where $g(c,\boldsymbol{\theta})=\log(\sigma/\pi)-\log(\sigma^{2}+(c-\mu)^{2})$
and $c({\bf u})={\bf x}^{T}{\bf u}$. We will estimate the influence function
for the leading Cauchy principal direction

$\hat{{\bf u}}=\operatorname*{\arg\min}_{||{\bf u}||_{2}=1}\ l({\bf
u}|\boldsymbol{\theta}_{F}({\bf u})),$ (12)

where $\boldsymbol{\theta}_{F}({\bf
u})=\operatorname*{\arg\max}_{\boldsymbol{\theta}}l({\bf x}^{T}{\bf
u}|\boldsymbol{\theta})$ is the optimal Cauchy parameters for a given
direction ${\bf u}$.

Since $\hat{{\bf u}}$ is restricted to be a unit vector, the standard
condition for the minimum, i.e., $\left.\frac{\partial}{\partial{\bf u}}l({\bf
u}|\boldsymbol{\theta}_{F}({\bf u}))\right|_{{\bf u}=\hat{{\bf u}}}={\bf 0}$
is not valid. The proper condition is defined by

${\bf P}_{\hat{{\bf u}}}\left.\frac{\partial}{\partial{\bf u}}l({\bf
u}|\boldsymbol{\theta}_{F}({\bf u}))\right|_{{\bf u}=\hat{{\bf u}}}={\bf 0},$
(13)

where ${\bf P}_{{\bf u}}$ is the projection matrix given by ${\bf P}_{{\bf
u}}={\bf I}_{p}-{\bf u}{\bf u}^{T}$.

######  Remark 3.1.

An equivalent condition is to satisfy ${\bf
h}^{T}\left.\frac{\partial}{\partial{\bf u}}l({\bf
u}|\boldsymbol{\theta}_{F}({\bf u}))\right|_{{\bf u}=\hat{{\bf u}}}={\bf 0}$
for all ${\bf h}$ such that ${\bf h}^{T}\hat{{\bf u}}=0$ and $||{\bf
h}||_{2}=1$. Both derived conditions are essentially a consequence of the
Lagrangian formulation of the constraint optimization problem. Indeed, the
Lagrange condition implies that at the minimum the direction of the objective
function’s derivative should be parallel to the direction of the constraint’s
derivative which translates to $\left.\frac{\partial}{\partial{\bf u}}l({\bf
u}|\boldsymbol{\theta}_{F}({\bf u}))\right|_{{\bf u}=\hat{{\bf
u}}}=\lambda\hat{{\bf u}}$ where $\lambda\neq 0$ is the Lagrange multiplier.

Let $\bar{g}({\bf x};{\bf u})=\left.g({\bf x}^{T}{\bf
u}|\theta)\right|_{\theta=\theta_{F}({\bf u})}$ be the likelihood function
computed at $\theta=\theta_{F}({\bf u})$ and let denote its partial
derivatives as

$\bar{g}_{c}({\bf x};{\bf u})=\left.\frac{\partial}{\partial c}g({\bf
x}^{T}{\bf u}|\theta)\right|_{\theta=\theta_{F}({\bf u})}$

and

$\bar{g}_{\boldsymbol{\theta}}({\bf x};{\bf
u})=\left.\frac{\partial}{\partial\boldsymbol{\theta}}g({\bf x}^{T}{\bf
u}|\theta)\right|_{\theta=\theta_{F}({\bf u})}.$

Similarly, $\bar{g}_{cc}$, $\bar{g}_{c\theta}$ and $\bar{g}_{\theta\theta}$
denote the second order derivatives. The following proposition establishes the
expression for the influence function of the leading Cauchy principal
direction, $\hat{{\bf u}}$.

######  Proposition 3.1.

Under the assumption of ${{\bf I}}_{F}(\hat{{\bf u}})$ and ${\bf A}$ being
invertible matrices, the influence function of $\hat{{\bf u}}$ is

$IF_{\hat{{\bf u}}}({\bf z},F)={\bf A}^{-1}{\bf b},$ (14)

where

$\displaystyle{\bf A}$ $\displaystyle={\bf
I}_{p}\int_{\mathbb{R}^{p}}\bar{g}_{c\boldsymbol{\theta}}({\bf x};\hat{{\bf
u}}){\bf x}^{T}\hat{{\bf u}}dF({\bf x})-{\bf P}_{\hat{{\bf
u}}}\int_{\mathbb{R}^{p}}\bar{g}_{cc}({\bf x};\hat{{\bf u}}){\bf x}^{T}{\bf
x}dF({\bf x}){\bf P}_{\hat{{\bf u}}}$ $\displaystyle+{\bf P}_{\hat{{\bf
u}}}\int_{\mathbb{R}^{p}}{\bf x}\bar{g}_{c\boldsymbol{\theta}}({\bf
x};\hat{{\bf u}})dF({\bf x})\,{{\bf I}}_{F}(\hat{{\bf
u}})^{-1}\,\int_{\mathbb{R}^{p}}\bar{g}_{\boldsymbol{\theta}c}({\bf
x};\hat{{\bf u}}){\bf x}^{T}dF({\bf x}){\bf P}_{\hat{{\bf u}}}$

and

${\bf b}={\bf b}(z)=\bar{g}_{c}({\bf z},\hat{{\bf u}}){\bf
z}+\int_{\mathbb{R}^{p}}{\bf x}\bar{g}_{c\boldsymbol{\theta}}({\bf
x};\hat{{\bf u}})dF({\bf x})\,{{\bf I}}_{F}(\hat{{\bf
u}})^{-1}\,\bar{g}_{\boldsymbol{\theta}}({\bf z};\hat{{\bf u}}),$

while

${{\bf I}}_{F}(\hat{{\bf
u}})=\int_{\mathbb{R}^{p}}\bar{g}_{\boldsymbol{\theta}\boldsymbol{\theta}}({\bf
x};\hat{{\bf u}})dF({\bf x})$

is the expected Fisher information matrix under $F$ for the parameters of the
Cauchy distribution computed at $\hat{{\bf u}}$.

###### Proof.

The proof consists of several straightforward series expansions and implicit
function calculations. The complete proof is given in Appendix A2. ∎

The following boundedness result for the influence function states the
conditions under which Cauchy PCA is robust.

######  Corollary 3.1.

Let the assumptions of the proposition hold. If ${\bf z}\not\perp\hat{{\bf
u}}$ or if ${\bf z}\perp\hat{{\bf u}}=0$ but $\mu_{F}(\hat{{\bf u}})=0$ then
the influence function for $\hat{{\bf u}}$ is bounded.

###### Proof.

First, observe that matrix ${\bf A}$ does not depend on ${\bf z}$. It is only
${\bf b}$ that depends on ${\bf z}$ and our goal is to prove that ${\bf b}$ is
bounded with respect to ${\bf z}$. Second, we have to compute the partial
derivatives $\bar{g}_{c}({\bf z};\hat{{\bf u}})$ and
$\bar{g}_{\boldsymbol{\theta}}({\bf z};\hat{{\bf u}})$. Straightforward
calculations lead to

$\bar{g}_{c}({\bf z};\hat{{\bf u}})=-\frac{2({\bf z}^{T}\hat{{\bf
u}}-\mu_{F}(\hat{{\bf u}}))}{\sigma_{F}^{2}(\hat{{\bf u}})+({\bf
z}^{T}\hat{{\bf u}}-\mu_{F}(\hat{{\bf u}}))^{2}}$ $\bar{g}_{\mu}({\bf
z};\hat{{\bf u}})=\frac{2({\bf z}^{T}\hat{{\bf u}}-\mu_{F}(\hat{{\bf
u}}))}{\sigma_{F}^{2}(\hat{{\bf u}})+({\bf z}^{T}\hat{{\bf
u}}-\mu_{F}(\hat{{\bf u}}))^{2}}$

and

$\bar{g}_{\sigma}({\bf z};\hat{{\bf u}})=\frac{1}{\sigma_{F}(\hat{{\bf
u}})}-\frac{2\sigma_{F}(\hat{{\bf u}})}{\sigma_{F}^{2}(\hat{{\bf u}})+({\bf
z}^{T}\hat{{\bf u}}-\mu_{F}(\hat{{\bf u}}))^{2}}.$

Let us now define an arbitrary scaling of the outlier ${\bf
z}\rightarrow\alpha{\bf z}$ and prove boundedness of ${\bf b}$ as we send
$\alpha\to\infty$. We consider the first case where ${\bf z}\not\perp\hat{{\bf
u}}$. It holds that $\lim_{\alpha\to\infty}\bar{g}_{c}(\alpha{\bf z};\hat{{\bf
u}})\alpha{\bf z}=-({\bf z}^{T}\hat{{\bf u}})^{-1}{\bf z}$,
$\lim_{\alpha\to\infty}\bar{g}_{\mu}(\alpha{\bf z};\hat{{\bf u}})=0$ and
$\lim_{\alpha\to\infty}\bar{g}_{\sigma}(\alpha{\bf z};\hat{{\bf
u}})=\frac{1}{\sigma_{F}(\hat{{\bf u}})}$ therefore ${\bf b}$ is bounded with
respect to $\alpha$.

For the second case, we have

$\lim_{\alpha\to\infty}\bar{g}_{c}(\alpha{\bf z};\hat{{\bf u}})\alpha{\bf
z}=\lim_{\alpha\to\infty}\frac{2\mu_{F}(\hat{{\bf
u}})}{\sigma_{F}^{2}(\hat{{\bf u}})+\mu_{F}(\hat{{\bf u}})^{2}}\alpha{\bf
z}=0,$ $\lim_{\alpha\to\infty}\bar{g}_{\mu}(\alpha{\bf z};\hat{{\bf
u}})=\frac{2\mu_{F}(\hat{{\bf u}})}{\sigma_{F}^{2}(\hat{{\bf
u}})+\mu_{F}(\hat{{\bf u}})^{2}}=0$

and

$\lim_{\alpha\to\infty}\bar{g}_{\sigma}(\alpha{\bf z};\hat{{\bf
u}})=\frac{1}{\sigma_{F}(\hat{{\bf u}})}-\frac{2\sigma_{F}(\hat{{\bf
u}})}{\sigma_{F}^{2}(\hat{{\bf u}})+\mu_{F}(\hat{{\bf
u}})^{2}}=-\frac{1}{\sigma_{F}(\hat{{\bf u}})}$

since $\mu_{F}(\hat{{\bf u}})=0$ by assumption. Thus ${\bf b}$ is bounded with
respect to $\alpha$ for the second case, too. ∎

The only case not covered by the corollary is when ${\bf z}^{T}\hat{{\bf
u}}=0$ and $\mu(\hat{{\bf u}})\neq 0$. Our experiments presented in the
following section show that outliers that are orthogonal to the Cauchy
principal direction do sometimes influence the estimation of the Cauchy
principal direction yet not significantly.

###  3.2 Several Cauchy principal components

We briefly mention possibilities for estimating several Cauchy principal
components. There are two obvious approaches: one approach, the sequential
approach, is to repeat the algorithm described above on the subspace
orthogonal to $\hat{{\bf u}}=\hat{{\bf u}}_{1}$ to obtain $\hat{{\bf u}}_{2}$,
the second Cauchy principal component, where $\hat{{\bf u}}_{1}$ is the first
Cauchy principal component; then repeat the procedure on the subspace
orthogonal to $\hat{{\bf u}}_{1}$ and $\hat{{\bf u}}_{2}$ to obtain $\hat{{\bf
u}}_{3}$; and so on. A second approach, the simultaneous approach, is to
decide in advance how many principal components we wish to determine, $p$ say,
and then use a $p$-dimensional multivariate Cauchy likelihood, which has
$p+p(p+1)/2$ free parameters, to obtain $\hat{{\bf u}}_{1},\ldots,\hat{{\bf
u}}_{p}$. It turns out that these two approaches lead to equivalent results in
classical (Gaussian) PCA but when a Cauchy likelihood is used the two
approaches produce different sets of principal components. Our current
thinking is this: the sequential approach is easier to implement (essentially
the same software can be used at each step) and it is faster. However, the
simultaneous approach could potentially be preferable if we know in advance
how many principal components we wish to estimate. Further investigation is
required.

##  4 Numerical Results

###  4.1 Simulation studies

In this section we will empirically validate our proposed methodology, via
simulation studies. We searched for R packages that offer robust PCA in the
$n<<p$ case and came up with FastHCS (Vakili,, 2018), rrcovHD (Todorov,,
2016), rpca (Sykulski,, 2017) and pcaPP (Filzmoser et al.,, 2018). Out of
them, pcaPP (Projection Pursuit PCA) is the only one which does not require
hyper-parameter tuning, e.g. selection of the LASSO penalty $\lambda$ or
choice of the percentage of observations used to estimate a robust covariance
matrix.

####  4.1.1 Setup of the simulations

Initially, we created a $p\times p$ (orthonormal) basis $\bf B$ by using QR
decomposition on some randomly generated data. We then generated eigenvalues
$\lambda_{i}\sim Exp(0.4)$, where $i=1,\ldots,p$ and hence we obtained the
covariance matrix $\boldsymbol{\Sigma}={\bf B}\boldsymbol{\Lambda}{\bf
B}^{T}$, where $\boldsymbol{\Lambda}=\text{diag}(\lambda_{i})$. The first
column of $\bf B$ served as the first “clean” eigenvector, and was the
benchmark in our comparative evaluations. Following this step, we simulated
$n$ random vectors ${\bf X}\sim N_{p}\left({\bf 0},\boldsymbol{\Sigma}\right)$
and in order to check the robustness of the results to the center of the data,
all observations were shifted right by adding $50$ everywhere. A number of
outliers equal to 2$\%$ of the sample size were introduced. These outliers
were $\bar{\bf x}+e^{\kappa}{\bf z}\in{\mathbb{R}}^{p}$, where $\bar{\bf x}$
is the sample mean vector, ${\bf z}$ are unit vector(s) and $e^{\kappa}$ a
real number denoting their norm, where $\kappa$ varied from $3$ up to $8$
increasing with a step size equal to $1$ and the angle between the outliers
${\bf z}$ and the first “clean” eigenvector spanned from $0^{\circ}$ up to
$90^{\circ}$. In all cases, we subtracted the spatial median or the column-
wise median222The results are pretty similar for either type of median and we
here show the results of he column-wise median. and scaled them by the mean
absolute deviation.

At each case, we computed the first Cauchy-PCA eigenvector and the first PP-
PCA eigenvector. The performance metric is the angle (in degrees) between the
first robust (based on Cauchy or PP-PCA) eigenvector and the first ”clean”
eigenvector computed using the classical PCA. All experiments were repeated
$100$ times and the results were averaged.

####  4.1.2 Comparative results

Tables 1-3 present the performance of the first Cauchy-PCA eigenvector and of
the first PP-PCA eigenvector for a variety of norms of the outlier, with
different angles ($\phi$) between the outlier and the leading true
eigenvector, for the $n<p$ case.

The case of $n<p$ was selected as statistical inference in this case is more
challenging than the $p<n$ case333In this paper we focus on high-dimensional
simulations and real-date examples ($p>n)$ but in results not presented in the
paper we found that Cauchy PCA is also very competitive and performs strongly
in low dimensional settings ($p<n$).. Additionally, this case is also
ordinarily met in the field of bioinformatics were the -omics data count tens
of thousands of variables (genes, single nucleotide polymorphisms, etc.) but
only tens or at most hundreds of observations.

As observed in Tables 1-3, the average angular difference between the Cauchy
and the PP PCA ranges from $20^{\circ}$ up to more than $50^{\circ}$, which is
evidently quite substantial, providing evidence that Cauchy PCA has performed
in a superior manner to the projection pursuit method of Croux et al. (2007,
2013). In particular, the tables demonstrate that Cauchy PCA is less error
prone than its competitor but, as is seen in Table 3, the error decreases for
both methods with increasing sample size. Further, the mean angular difference
between the two methods increases as the angle $\phi$ increases. For instance,
in Table 1, when $k=8$ and $\phi=0^{\circ}$ the difference between the two
methods is $20^{\circ}$, whereas when $\phi=90^{\circ}$ the difference
increases to $48^{\circ}$. Further, the error is not highly affected by the
angle $\phi$, or the norm of the outliers. It can be seen that in Table 2 and
Table 3 in the special case of $\phi=90^{\circ}$, the error increases for the
Cauchy PCA by $2^{\circ}-3^{\circ}$, thus corroborating the result of
Corollary 3.1. However, this effect, as in Table 1, is rather small, though
noticeable.

Table 1: Mean angular difference between the robust eigenvectors computed in the contaminated data and the sample eigenvector computed in the clean data when $n=100$ and $p=500$. The norm of the outliers is $e^{k}$ and their angle with the true clean eigenvector is denoted by $\phi$. Angle | Method | k=-Inf | k=3 | k=4 | k=5 | k=6 | k=7 | k=8
---|---|---|---|---|---|---|---|---
$\phi=0^{\circ}$ | Cauchy | 31.17 | 29.79 | 29.54 | 28.83 | 28.86 | 29.24 | 28.78
| PP | 82.45 | 49.91 | 48.84 | 48.22 | 49.08 | 49.61 | 48.14
$\phi=30^{\circ}$ | Cauchy | 31.44 | 29.24 | 29.13 | 28.60 | 28.89 | 29.34 | 29.65
| PP | 82.45 | 65.28 | 65.34 | 63.42 | 62.96 | 66.63 | 65.43
$\phi=60^{\circ}$ | Cauchy | 31.49 | 29.86 | 29.07 | 29.04 | 29.55 | 29.70 | 29.09
| PP | 82.11 | 81.11 | 82.55 | 82.63 | 82.12 | 82.49 | 82.03
$\phi=90^{\circ}$ | Cauchy | 32.32 | 31.67 | 33.00 | 33.13 | 32.86 | 33.19 | 33.06
| PP | 82.38 | 82.06 | 81.69 | 82.12 | 81.73 | 81.74 | 81.88
Table 2: Mean angular difference between the robust eigenvectors computed in the contaminated data and the sample eigenvector computed in the clean data when $n=100$ and $p=1000$. The norm of the outliers is $e^{k}$ and their angle with the true clean eigenvector is denoted by $\phi$. Angle | Method | k=-Inf | k=3 | k=4 | k=5 | k=6 | k=7 | k=8
---|---|---|---|---|---|---|---|---
$\phi=0^{\circ}$ | Cauchy | 36.53 | 33.12 | 33.60 | 33.69 | 32.62 | 32.51 | 33.16
| PP | 83.06 | 80.36 | 80.17 | 81.87 | 80.50 | 80.76 | 80.16
$\phi=30^{\circ}$ | Cauchy | 36.55 | 34.72 | 33.91 | 33.09 | 33.11 | 33.16 | 32.79
| PP | 83.07 | 82.36 | 82.76 | 82.65 | 83.07 | 82.93 | 83.12
$\phi=60^{\circ}$ | Cauchy | 36.42 | 34.46 | 33.96 | 33.61 | 34.41 | 33.07 | 33.47
| PP | 83.78 | 82.86 | 82.71 | 84.05 | 83.46 | 82.71 | 82.78
$\phi=90^{\circ}$ | Cauchy | 36.50 | 36.12 | 36.81 | 37.18 | 39.34 | 39.11 | 38.51
| PP | 83.63 | 83.73 | 83.69 | 83.65 | 84.03 | 83.66 | 83.00
Table 3: Mean angular difference between the robust eigenvectors computed in the contaminated data and the sample eigenvector computed in the clean data when $n=500$ and $p=1000$. The norm of the outliers is $e^{k}$ and their angle with the true clean eigenvector is denoted by $\phi$. Angle | Method | k=-Inf | k=3 | k=4 | k=5 | k=6 | k=7 | k=8
---|---|---|---|---|---|---|---|---
$\phi=0^{\circ}$ | Cauchy | 19.95 | 18.60 | 18.46 | 18.35 | 18.24 | 18.20 | 17.93
| PP | 68.76 | 26.08 | 24.93 | 24.91 | 24.83 | 24.73 | 24.72
$\phi=30^{\circ}$ | Cauchy | 19.43 | 18.30 | 18.39 | 18.22 | 18.16 | 18.01 | 18.13
| PP | 68.98 | 39.72 | 38.88 | 38.44 | 38.20 | 38.15 | 38.14
$\phi=60^{\circ}$ | Cauchy | 19.76 | 18.60 | 18.12 | 18.20 | 18.40 | 18.19 | 18.01
| PP | 69.10 | 64.10 | 63.12 | 62.89 | 62.91 | 62.82 | 62.77
$\phi=90^{\circ}$ | Cauchy | 19.49 | 19.84 | 20.16 | 21.87 | 22.41 | 22.87 | 22.84
| PP | 68.99 | 68.62 | 68.59 | 68.70 | 68.45 | 68.73 | 68.43

###  4.2 High dimensional real datasets

Two real gene expression datasets, GSE13159 and GSE31161444From a biological
standpoint, the data have already been uniformly pre-processed, curated and
automatically annotated., downloaded from the Biodataome platform (Lakiotaki
et al.,, 2018), were used in the experiments. The dimensions of the datasets
were equal to $2,096\times 54,630$ and $1035\times 54,675$, respectively. We
randomly selected $5,000$ variables and computed the outliers using the high
dimensional Minimum Covariance Determinant (MCD) of Ro et al., (2015). In
accordance with the simulations studies, we removed the $2\%$ of the most
extreme outliers detected by MCD and computed the first classical PC
(benchmark eigenvector), the first Cauchy-PCA eigenvector and the first PP-PCA
eigenvector of the ”clean” data. We then added those outliers and increased
their norm by $e^{k}$, where $k=(0,3,4,\ldots,8)$ and computed computed the
first Cauchy-PCA eigenvector and the first PP-PCA eigenvector. In all cases,
we subtracted the spatial median or the column-wise median and scaled them by
the mean absolute deviation. The performance metric is the angle (in degrees)
between the first robust (based on Cauchy or PP-PCA) eigenvector and the first
true “clean” eigenvectors and the time required by each method. This procedure
was repeated $200$ times and the average results are graphically displayed in
Figures 1(a)-(d).

Broadly speaking the effect of the PP PCA does not seem to have been affected
substantially by the centering method, i.e. subtraction of the spatial or the
column-wise median. On the contrary, the Cauchy PCA is affected by the type of
median employed to this end. Centering with the spatial median yields high
error levels for all norms of the outliers, for both datasets, whereas
centering with the column-wise median produces much lower error levels. On
average, the difference in the error between Cauchy PCA and PP PCA is about
$30^{\circ}$ for the GSE31159 dataset (Figure 1(a)) and about $14^{\circ}$ for
the GSE3161 dataset (Figure 1(b)). However, the error of the Cauchy PCA
increases and the stabilizes in the GSE31159 dataset whereas the error of the
PP PCA is stable regardless of the norm of the outliers. A different
conclusion is extracted in the GSE31161 where the error of either method
decreases as the norm of the outliers increases, until it reaches a plateau.

With regards to computational efficiency, the PP PCA is not affected by either
centering method, whereas Cauchy PCA seems to be affected in the GSE31159
dataset but not in the GSE31161 dataset as seen in Figures 1(c) and 1(d).
Cauchy PCA centered with the column-wise median is, on average, 5 times faster
than PP PCA.

GSE31159 | GSE31161
---|---
|
(a) | (b)
|
(c) | (d)
Figure 1: The first row presents the angle between the first Cauchy PC of the
”contaminated” data and the 1st leading eigenvector of the ”clean” data and
the angle between the first Projection Pursuit PC of the ”contaminated” data
and the 1st leading eigenvector of the ”clean” data for increasing norms of
the outliers. The second row contains the time in seconds.

##  5 Conclusion

The starting point for this paper is the observation that classical PCA can be
formulated purely in terms of operations on a Gaussian likelihood. Although
this observation is not new, the specifics of this formulation of classical
PCA do not appear to be as widely known as might be expected. The novel idea
underlying this paper is to formulate a version of PCA in which a Cauchy
likelihood is used instead of a Gaussian likelihood, leading to what we call
Cauchy PCA. Study of the resulting influence functions shows that Cauchy PCA
has very good robustness properties. Moreover, we have provided an
implementation of Cauchy PCA which runs quickly and reliably. Numerous
simulation and real-data examples, mainly in high-dimensional settings, show
that Cauchy PCA typically out-performs alternative robust versions of PCA
whose implementation is in the public domain.

## Appendix

###  A1 Proof of Proposition 2.1

###### Proof.

The perturbed distribution $(1-\epsilon)F({\bf x})+\epsilon\Delta_{\bf z}({\bf
x})$ has perturbed mean value

$\boldsymbol{\mu}_{\epsilon}=\boldsymbol{\mu}+\epsilon({\bf
z}-\boldsymbol{\mu})$

and perturbed covariance matrix

$\boldsymbol{\Sigma}_{\epsilon}=\boldsymbol{\Sigma}+\epsilon(({\bf
z}-\boldsymbol{\mu})({\bf
z}-\boldsymbol{\mu})^{T}-\boldsymbol{\Sigma})+\epsilon^{2}({\bf
z}-\boldsymbol{\mu})({\bf z}-\boldsymbol{\mu})^{T}$

Denoting by $\lambda_{\epsilon}$ the leading eigenvalue of
$\boldsymbol{\Sigma}_{\epsilon}$ and by ${\bf u}_{\epsilon}$ the corresponding
eigenvector, it holds that

$\boldsymbol{\Sigma}_{\epsilon}{\bf u}_{\epsilon}=\lambda_{\epsilon}{\bf
u}_{\epsilon}\ \ \text{and}\ \ {\bf u}_{\epsilon}^{T}{\bf u}_{\epsilon}=1\ .$
(15)

Next, we expand the perturbed eigenvector and eigenvalue around the
unperturbed ones as follows:

${\bf u}_{\epsilon}={\bf u}_{0}+\epsilon{\bf u}_{1}+O(\epsilon^{2})$

and

$\lambda_{\epsilon}=\lambda_{0}+\epsilon\lambda_{1}+O(\epsilon^{2})$

with

$\boldsymbol{\Sigma}{\bf u}_{0}=\lambda_{0}\ \ \text{and}\ \ {\bf
u}_{0}^{T}{\bf u}_{0}=1\ .$

Substituting the formulas into (15), and equating the zero-th and first order
we get

$\boldsymbol{\Sigma}{\bf u}_{0}=\lambda_{0}{\bf u}_{0}\ \ \text{and}\ \ {\bf
u}_{0}^{T}{\bf u}_{0}=1\ .$

and

$(({\bf z}-\boldsymbol{\mu})({\bf
z}-\boldsymbol{\mu})^{T}-\boldsymbol{\Sigma}){\bf
u}_{0}+\boldsymbol{\Sigma}{\bf u}_{1}=\lambda_{0}{\bf u}_{1}+\lambda_{1}{\bf
u}_{0}$ (16)

and

${\bf u}_{0}^{T}{\bf u}_{1}=0\ .$

Multiplying (16) from the left with ${\bf u}_{0}^{T}$, we get

$\lambda_{1}={\bf u}_{0}^{T}(({\bf z}-\boldsymbol{\mu})({\bf
z}-\boldsymbol{\mu})^{T}-\boldsymbol{\Sigma}){\bf u}_{0}+{\bf
u}_{0}^{T}\boldsymbol{\Sigma}{\bf u}_{1}=({\bf u}_{0}^{T}({\bf
z}-\boldsymbol{\mu}))^{2}-\lambda_{0}$

For ${\bf u}_{1}$, we rearrange (16) to

$(\boldsymbol{\Sigma}-\lambda_{0}{\bf I}){\bf u}_{1}=\lambda_{1}{\bf
u}_{0}-(({\bf z}-\boldsymbol{\mu})({\bf
z}-\boldsymbol{\mu})^{T}-\boldsymbol{\Sigma}){\bf u}_{0}$

and then multiply from the left with the pseudo-inverse of
$\boldsymbol{\Sigma}-\lambda_{0}{\bf I}$ to obtain

$(\boldsymbol{\Sigma}-\lambda_{0}{\bf
I})^{+}(\boldsymbol{\Sigma}-\lambda_{0}{\bf I}){\bf
u}_{1}=\lambda_{1}(\boldsymbol{\Sigma}-\lambda_{0}{\bf I})^{+}{\bf
u}_{0}-(\boldsymbol{\Sigma}-\lambda_{0}{\bf I})^{+}(({\bf
z}-\boldsymbol{\mu})({\bf z}-\boldsymbol{\mu})^{T}-\boldsymbol{\Sigma}){\bf
u}_{0}$

Using the properties (Mardia et al., (1979)):
$(\boldsymbol{\Sigma}-\lambda_{0}{\bf
I})^{+}(\boldsymbol{\Sigma}-\lambda_{0}{\bf I})={\bf I}-{\bf u}_{0}{\bf
u}_{0}^{T}$ and $(\boldsymbol{\Sigma}-\lambda_{0}{\bf I})^{+}{\bf u}_{0}={\bf
0}$, we obtain

$\displaystyle{\bf u}_{1}-{\bf u}_{0}{\bf u}_{0}^{T}{\bf
u}_{1}=(\boldsymbol{\Sigma}-\lambda_{0}{\bf I})^{+}({\bf
z}-\boldsymbol{\mu})({\bf z}-\boldsymbol{\mu})^{T}{\bf
u}_{0}-(\boldsymbol{\Sigma}-\lambda_{0}{\bf I})^{+}\lambda_{0}{\bf u}_{0}$
$\displaystyle\Rightarrow$ $\displaystyle{\bf u}_{1}=(({\bf
z}-\boldsymbol{\mu})^{T}{\bf u}_{0})(\boldsymbol{\Sigma}-\lambda_{0}{\bf
I})^{+}({\bf z}-\boldsymbol{\mu})$

and the proof is completed. ∎

###  A2 Proof of Proposition 3.1

Let us first make the symbolism more explicit and denote $l_{F}({\bf
u}|\boldsymbol{\theta})$ the Cauchy log-likelihood function with respect to
the distribution $F$ and $\hat{{\bf u}}_{F}$ the respective leading Cauchy
principal direction. Then, our goal is to calculate the limit of

$\frac{1}{\epsilon}(\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}-\hat{{\bf u}}_{F})$

as $\epsilon\to 0$ where $\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}$ is the leading
Cauchy principal direction for the distribution $F_{\epsilon,{\bf
z}}=(1-\epsilon)F+\epsilon\Delta_{\bf z}$. The optimality condition for the
leading Cauchy principal direction reads

${\bf P}_{\hat{{\bf u}}_{F_{\epsilon,{\bf
z}}}}\left.\frac{\partial}{\partial{\bf u}}l_{F_{\epsilon,{\bf z}}}\big{(}{\bf
u}|\boldsymbol{\theta}_{F_{\epsilon,{\bf z}}}({\bf u})\big{)}\right|_{{\bf
u}=\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}}=0$ (17)

and

${\bf P}_{\hat{{\bf u}}_{F}}\left.\frac{\partial}{\partial{\bf
u}}l_{F}\big{(}{\bf u}|\boldsymbol{\theta}_{F}({\bf u})\big{)}\right|_{{\bf
u}=\hat{{\bf u}}_{F}}=0$

Moreover, $\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}$ is a unit vector which can be
represented as

$\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}=\cos(\rho)\hat{{\bf
u}}_{F}+\sin(\rho){\bf h}$

where ${\bf h}$ is a unit vector perpendicular to $\hat{{\bf u}}_{F}$ and
$\rho$ is a (small) real number. Under these assumptions, $\hat{{\bf
u}}_{F_{\epsilon,{\bf z}}}$ is a unit vector since

$||\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}||_{2}^{2}=\cos^{2}(\rho)||\hat{{\bf
u}}_{F}||_{2}^{2}+\sin^{2}(\rho)||{\bf h}||_{2}^{2}=1$

Obviously, $\rho$ depends on $\epsilon$ and ${\bf z}$ (i.e.,
$\rho=\rho(\epsilon,{\bf z})$) and $\lim_{\epsilon\to 0}\rho=0$ but we choose
to avoid denoting their explicit relationship because it is not required in
our proof. Moreover, a Taylor expansion for the representation leads to

$\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}=\hat{{\bf u}}_{F}+\rho{\bf
h}+O(\rho^{2})$

thus we obtain that

${\bf P}_{\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}}={\bf P}_{\hat{{\bf
u}}_{F}}-\rho(\hat{{\bf u}}_{F}{\bf h}^{T}+{\bf h}\hat{{\bf
u}}_{F}^{T})+O(\rho^{2})$

Next, we compute the partial derivative using the chain rule

$\frac{\partial}{\partial{\bf u}}l_{F}\big{(}{\bf
u}|\boldsymbol{\theta}_{F}({\bf
u})\big{)}=\int_{\mathbb{R}^{p}}\left[\frac{\partial}{\partial c}g(c({\bf
u}),\boldsymbol{\theta}_{F}({\bf u}))\frac{\partial}{\partial{\bf u}}c({\bf
u})+\frac{\partial}{\partial\boldsymbol{\theta}}g(c({\bf
u}),\boldsymbol{\theta}_{F}({\bf u}))\frac{\partial}{\partial{\bf
u}}\boldsymbol{\theta}_{F}({\bf u})\right]dF({\bf x})$

Therefore,

$\displaystyle\left.\frac{\partial}{\partial{\bf u}}l_{F}\big{(}{\bf
u}|\boldsymbol{\theta}_{F}({\bf u})\big{)}\right|_{{\bf u}=\hat{{\bf u}}_{F}}$
$\displaystyle=\int_{\mathbb{R}^{p}}\left[\bar{g}_{c}({\bf x};\hat{{\bf
u}}_{F}){\bf x}+\bar{g}_{\boldsymbol{\theta}}({\bf x};\hat{{\bf
u}}_{F})\frac{\partial}{\partial{\bf u}}\boldsymbol{\theta}_{F}({\bf
u})\Big{|}_{{\bf u}=\hat{{\bf u}}_{F}}\right]dF({\bf x})$
$\displaystyle=\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F}){\bf
x}dF({\bf x})+\int_{\mathbb{R}^{p}}\bar{g}_{\boldsymbol{\theta}}({\bf
x};\hat{{\bf u}}_{F})dF({\bf x})\frac{\partial}{\partial{\bf
u}}\boldsymbol{\theta}_{F}({\bf u})\Big{|}_{{\bf u}=\hat{{\bf u}}_{F}}$
$\displaystyle=\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F}){\bf
x}dF({\bf x})$

The second summand equals to zero because $\hat{{\bf u}}_{F}$ maximizes the
Cauchy log-likelihood function thus it holds that
$\int_{\mathbb{R}^{p}}\bar{g}_{\boldsymbol{\theta}}({\bf x};\hat{{\bf
u}}_{F})dF({\bf x})={\bf 0}$. Similarly,

$\displaystyle\left.\frac{\partial}{\partial{\bf u}}l_{F_{\epsilon,{\bf
z}}}\big{(}{\bf u}|\boldsymbol{\theta}_{F_{\epsilon,{\bf z}}}({\bf
u})\big{)}\right|_{{\bf u}=\hat{{\bf u}}_{F_{\epsilon,{\bf
z}}}}=\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F_{\epsilon,{\bf
z}}}){\bf x}dF_{\epsilon,{\bf z}}({\bf x})$ $\displaystyle=$
$\displaystyle(1-\epsilon)\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf
u}}_{F_{\epsilon,{\bf z}}}){\bf x}dF({\bf x})+\epsilon\bar{g}_{c}({\bf
z};\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}){\bf z}$

Next, we further Taylor expand $\bar{g}_{c}({\bf x};{\bf u}_{F_{\epsilon,{\bf
z}}})$ using $\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}=\hat{{\bf u}}_{F}+\rho{\bf
h}+O(\rho^{2})$

$\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F_{\epsilon,{\bf z}}})=\bar{g}_{c}({\bf
x};\hat{{\bf u}}_{F})+\rho{\bf h}\frac{\partial}{\partial{\bf
u}}\bar{g}_{c}({\bf x};{{\bf u}})\Big{|}_{{\bf u}=\hat{{\bf
u}}_{F}}+O(\rho^{2})$

Using again the chain rule, we obtain that

$\frac{\partial}{\partial{\bf u}}\bar{g}_{c}({\bf x};{{\bf
u}})=\bar{g}_{cc}({\bf x};{{\bf u}}){\bf
x}+\bar{g}_{c\boldsymbol{\theta}}({\bf x};{{\bf
u}})\frac{\partial}{\partial{\bf u}}\boldsymbol{\theta}_{F}({\bf u})$

The computation of the partial derivative $\frac{\partial}{\partial{\bf
u}}\boldsymbol{\theta}_{F}({\bf u})$ follows. Formula
$\boldsymbol{\theta}_{F}({\bf
u})=\operatorname*{\arg\max}_{\boldsymbol{\theta}}l_{F}({\bf x}^{T}{\bf
u}|\boldsymbol{\theta})$ implies that

$\frac{\partial}{\partial\boldsymbol{\theta}}l_{F}(c({\bf
u})|\boldsymbol{\theta})\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta}_{F}({\bf
u})}=0\ .$

Differentiating with respect to ${\bf u}$ and using the implicit function
theorem, we get

$\displaystyle\frac{\partial}{\partial{\bf u}}\boldsymbol{\theta}_{F}({\bf
u})$ $\displaystyle=-\frac{\partial}{\partial{\bf
u}}\frac{\partial}{\partial\boldsymbol{\theta}}l_{F}(c({\bf
u})|\boldsymbol{\theta})\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta}_{F}({\bf
u})}\left[\frac{\partial^{2}}{\partial\boldsymbol{\theta}^{2}}l_{F}(c({\bf
u})|\boldsymbol{\theta})\Big{|}_{\boldsymbol{\theta}=\boldsymbol{\theta}_{F}({\bf
u})}\right]^{-1}$ $\displaystyle=-\int_{\mathbb{R}^{p}}{\bf
x}\bar{g}_{c\boldsymbol{\theta}}({\bf x};{\bf u})dF({\bf
x})\left[\int_{\mathbb{R}^{p}}\bar{g}_{\boldsymbol{\theta}\boldsymbol{\theta}}({\bf
x};{\bf u})dF({\bf x})\right]^{-1}$

Thus,

$\displaystyle\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F_{\epsilon,{\bf
z}}})=\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F})$ $\displaystyle+$
$\displaystyle\rho{\bf h}\left[\bar{g}_{cc}({\bf x};\hat{{\bf u}}_{F}){\bf
x}+\int_{\mathbb{R}^{p}}{\bf x}\bar{g}_{c\boldsymbol{\theta}}({\bf
x};\hat{{\bf u}}_{F})dF({\bf
x})\left[\int_{\mathbb{R}^{p}}\bar{g}_{\boldsymbol{\theta}\boldsymbol{\theta}}({\bf
x};\hat{{\bf u}}_{F})dF({\bf
x})\right]^{-1}\bar{g}_{c\boldsymbol{\theta}}({\bf x};{\hat{{\bf
u}}_{F}})\right]+O(\rho^{2})$

Overall, (17) becomes

$\displaystyle\left[{\bf P}_{\hat{{\bf u}}_{F}}-\rho(\hat{{\bf u}}_{F}{\bf
h}^{T}+{\bf h}\hat{{\bf
u}}_{F}^{T})+O(\rho^{2})\right]\cdot\left[(1-\epsilon)\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf
x};\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}){\bf x}dF({\bf
x})+\epsilon\bar{g}_{c}({\bf z};\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}){\bf
z}\right]=0$ $\displaystyle\Rightarrow$ $\displaystyle{\bf P}_{\hat{{\bf
u}}_{F}}\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf
u}}_{F_{\epsilon,{\bf z}}}){\bf x}dF({\bf x})-\rho(\hat{{\bf u}}_{F}{\bf
h}^{T}+{\bf h}\hat{{\bf u}}_{F}^{T})\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf
x};\hat{{\bf u}}_{F}){\bf x}dF({\bf x})+O(\rho^{2})$
$\displaystyle=\epsilon{\bf P}_{\hat{{\bf
u}}_{F}}\left[\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F}){\bf
x}dF({\bf x})-\bar{g}_{c}({\bf z};\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}){\bf
z}\right]+O(\epsilon\rho)$ $\displaystyle\Rightarrow$ $\displaystyle\rho{\bf
h}\left[\int_{\mathbb{R}^{p}}\bar{g}_{cc}({\bf x};\hat{{\bf u}}_{F}){\bf
x}^{T}{\bf x}dF({\bf x})+\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf
u}}_{F})\hat{{\bf u}}_{F}^{T}{\bf x}dF({\bf x})\right.$
$\displaystyle+\left.\int_{\mathbb{R}^{p}}{\bf
x}\bar{g}_{c\boldsymbol{\theta}}({\bf x};\hat{{\bf u}}_{F})dF({\bf
x})\left[\int_{\mathbb{R}^{p}}\bar{g}_{\boldsymbol{\theta}\boldsymbol{\theta}}({\bf
x};\hat{{\bf u}}_{F})dF({\bf
x})\right]^{-1}\int_{\mathbb{R}^{p}}\bar{g}_{c\boldsymbol{\theta}}({\bf
x};{\hat{{\bf u}}_{F}}){\bf x}dF({\bf x})\right]+O(\rho^{2})$
$\displaystyle=\epsilon{\bf P}_{\hat{{\bf
u}}_{F}}\left[\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F}){\bf
x}dF({\bf x})-\bar{g}_{c}({\bf z};\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}){\bf
z}\right]+O(\epsilon\rho)$

where we use the facts that

${\bf P}_{\hat{{\bf u}}_{F}}{\bf h}={\bf h}$

and

${\bf h}^{T}\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F}){\bf
x}dF({\bf x})={\bf P}_{\hat{{\bf u}}_{F}}\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf
x};\hat{{\bf u}}_{F}){\bf x}dF({\bf x})={\bf P}_{\hat{{\bf
u}}_{F}}\left.\frac{\partial}{\partial{\bf u}}l_{F}\big{(}{\bf
u}|\boldsymbol{\theta}_{F}({\bf u})\big{)}\right|_{{\bf u}=\hat{{\bf
u}}_{F}}=0$

Thus, the influence function is

$IF_{\hat{{\bf u}}_{F}}({\bf z},F)=\lim_{\epsilon\to 0}\frac{\rho{\bf
h}}{\epsilon}={\bf A}^{-1}{\bf b}$

where

$\displaystyle{\bf A}$ $\displaystyle={\bf
I}_{d}\left[\int_{\mathbb{R}^{p}}\bar{g}_{cc}({\bf x};\hat{{\bf u}}_{F}){\bf
x}^{T}{\bf x}dF({\bf x})+\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf
u}}_{F})\hat{{\bf u}}_{F}^{T}{\bf x}dF({\bf x})\right]$
$\displaystyle+\int_{\mathbb{R}^{p}}{\bf x}\bar{g}_{c\boldsymbol{\theta}}({\bf
x};\hat{{\bf u}}_{F})dF({\bf
x})\left[\int_{\mathbb{R}^{p}}\bar{g}_{\boldsymbol{\theta}\boldsymbol{\theta}}({\bf
x};\hat{{\bf u}}_{F})dF({\bf
x})\right]^{-1}\int_{\mathbb{R}^{p}}\bar{g}_{c\boldsymbol{\theta}}({\bf
x};{\hat{{\bf u}}_{F}}){\bf x}dF({\bf x})$

and

${\bf b}={\bf P}_{\hat{{\bf
u}}_{F}}\left[\int_{\mathbb{R}^{p}}\bar{g}_{c}({\bf x};\hat{{\bf u}}_{F}){\bf
x}dF({\bf x})-\bar{g}_{c}({\bf z};\hat{{\bf u}}_{F_{\epsilon,{\bf z}}}){\bf
z}\right]$

## References

  * Bolton and Krzanowski, (1999) Bolton, R. J. and Krzanowski, W. J. (1999).  A Characterization of Principal Components for Projection Pursuit.  American Statistican, 53(2):108–109.
  * Campbell, (1980) Campbell, N. A. (1980).  Robust Procedures in Multivariate Analysis I: Robust Covariance Estimation.  Applied Statistics, 29(3):231–237.
  * Copas, (1975) Copas, J. B. (1975).  On the Unimodality of the Likelihood for the Cauchy Distribution.  Biometrika, 62(3):701–704.
  * Croux et al., (2013) Croux, C., Filzmoser, P., and Fritz, H. (2013).  Robust sparse principal component analysis.  Technometrics, 55(2):202–214.
  * Croux et al., (2007) Croux, C., Filzmoser, P., and Oliveira, M. R. (2007).  Algorithms for projection–pursuit robust principal component analysis.  Chemometrics and Intelligent Laboratory Systems, 87(2):218–225.
  * Filzmoser et al., (2018) Filzmoser, P., Fritz, H., and Kalcher, K. (2018).  pcaPP: Robust PCA by Projection Pursuit.  R package version 1.9-73.
  * Huber, (1964) Huber, P. J. (1964).  Robust Estimation of a Location Parameter.  The Annals of Mathematical Statistics, 35(1):73–101.
  * Hubert and Verboven, (2003) Hubert, M. and Verboven, S. (2003).  A Robust PCR Method for High-Dimensional Regressors.  Journal of Chemometrics, 17:438–452.
  * Lakiotaki et al., (2018) Lakiotaki, K., Vorniotakis, N., Tsagris, M., Georgakopoulos, G., and Tsamardinos, I. (2018).  Biodataome: a collection of uniformly preprocessed and automatically annotated datasets for data-driven biology.  Database, 2018.
  * Li and Chen, (1985) Li, G. Y. and Chen, Z. L. (1985).  Projection-Pursuit Approach to Robust Dispersion Matrices and Principal Components: Primary Theory and Monte Carlo.  Journal of the American Statistical Association, 80:759–766.
  * Mardia et al., (1979) Mardia, K. F., Kent, J. T., and Bibby, J. M. (1979).  Multivariate Analysis.  Academic Press.
  * Maronna, (1976) Maronna, R. A. (1976).  Robust M-Estimators of Multivariate Location and Scatter.  The Annals of Statistics, 4(1):51–67.
  * Ro et al., (2015) Ro, K., Zou, C., Wang, Z., and Yin, G. (2015).  Outlier detection for high-dimensional data.  Biometrika, 102(3):589–599.
  * Sykulski, (2017) Sykulski, M. (2017).  rpca: RobustPCA: Decompose a Matrix into Low-Rank and Sparse Components.  R package version 0.2.3.
  * Todorov, (2016) Todorov, V. (2016).  rrcovHD: Robust Multivariate Methods for High Dimensional Data.  R package version 0.2-5.
  * Vakili, (2018) Vakili, K. (2018).  FastHCS: Robust Algorithm for Principal Component Analysis.  R package version 0.0.6.
  * Xie et al., (1993) Xie, Y. L., Wang, J. H., Liang, L. X., Sun, X. H. S., and Yu, R. Q. (1993).  Robust Principal Component Analysis By Projection-Pursuit.  Journal of Chemometrics, 7(6):527–541.