# Fast robust location and scatter estimation: a depth-based method

Maoyu Zhang Email<EMAIL_ADDRESS>Center for Applied Statistics,
Institute of Statistics and Big Data, Renmin University of China, Beijing,
China Yan Song Email<EMAIL_ADDRESS>Center for Applied Statistics,
Institute of Statistics and Big Data, Renmin University of China, Beijing,
China Wenlin Dai Corresponding author; email<EMAIL_ADDRESS>Center
for Applied Statistics, Institute of Statistics and Big Data, Renmin
University of China, Beijing, China

###### Abstract

The minimum covariance determinant (MCD) estimator is ubiquitous in
multivariate analysis, the critical step of which is to select a subset of a
given size with the lowest sample covariance determinant. The concentration
step (C-step) is a common tool for subset-seeking; however, it becomes
computationally demanding for high-dimensional data. To alleviate the
challenge, we propose a depth-based algorithm, termed as FDB, which replaces
the optimal subset with the trimmed region induced by statistical depth. We
show that the depth-based region is consistent with the MCD-based subset under
a specific class of depth notions, for instance, the projection depth. With
the two suggested depths, the FDB estimator is not only computationally more
efficient but also reaches the same level of robustness as the MCD estimator.
Extensive simulation studies are conducted to assess the empirical performance
of our estimators. We also validate the computational efficiency and
robustness of our estimators under several typical tasks such as principal
component analysis, linear discriminant analysis, image denoise and outlier
detection on real-life datasets. A R package FDB and potential extensions are
available in the Supplementary Materials.

Keywords: Computationally efficient; High-dimensional data; Outliers;
Robustness; Statistical depth.

##  1 Introduction

The Minimum Covariance Determinant (MCD) estimator (Rousseeuw, 1984) is among
the first affine equivariant and highly robust estimators of multivariate
location and scatter. Specifically, for a collection of multivariate data, MCD
seeks a subset of samples that leads to a sample covariance matrix with the
minimum determinant out of all the candidate sets of a specific size. The
location and scatter estimators are then defined as the average and a scaled
covariance matrix of these samples, respectively. Butler et al. (1993) and
Cator and Lopuhaä (2012) established the consistency and asymptotic normality
of the MCD estimator.

MCD has been applied in various fields such as quality control, medicine,
finance, image analysis, and chemistry (Hubert et al., 2008, 2018). Estimating
the covariance matrix is the cornerstone of many multivariate statistical
methods, so MCD has also been used to develop robust and computationally
efficient multivariate techniques, such as principal component analysis (Croux
and Haesbroeck, 2000; Hubert et al., 2005b), factor analysis (Pison et al.,
2003), classification (Hubert and Van Driessen, 2004), clustering (Hardin and
Rocke, 2004), multivariate regression (Rousseeuw et al., 2004b), and others
(Hubert et al., 2008). To cater to its broad applications, extensive effort
has been made to improve the computational efficiency of the approximation
algorithm. For example, Rousseeuw and Driessen (1999) propose the first
computationally efficient algorithm, termed FASTMCD; (Hubert et al., 2012)
suggest an improved version of FASTMCD, termed DetMCD; De Ketelaere et al.
(2020) accelerates DetMCD by refinement of the calculation steps and parallel
computation. Furthermore, Boudt et al. (2020) generalizes the MCD to high-
dimensional cases as the minimum regularized covariance determinant (MRCD).
Other variants include the orthogonalized Gnanadesikan-Kettenring (Maronna and
Zamar, 2002), the minimum (regularized) weighted covariance determinant
(Roelant et al., 2009; Kalina and Tichavskỳ, 2021), and kernel MRCD for non-
elliptical data (Schreurs et al., 2021).

Practically, the MCD-type algorithms are limited by two factors. First, the
computational complexity of the concentration step (C-step), which is critical
for such algorithms, is $O(np^{2}+p^{3})$, and this severely limits the
scalability of the algorithm for massive high-dimensional data. Second, the
approximation to the true MCD subset becomes less accurate due to the curse of
dimensionality. We note that the asymptotic trimmed region induced by a class
of statistical depth shares the same form with the asymptotic MCD subset when
the data are elliptically symmetric distributed (Butler et al., 1993; Zuo and
Serfling, 2000b). This motivates us to investigate the possibility of finding
the MCD subset directly by utilizing statistical depth. Statistical depth was
first considered for ranking multivariate data from the center outward
(Mahalanobis, 1936; Tukey, 1975; Oja, 1983; Liu, 1990; Zuo and Serfling,
2000a; Vardi and Zhang, 2000). Usually, a statistical depth is an increasing
function of the centrality of observations, taking values in $[0,1]$.

Motivated by the connection mentioned above, we propose a fast depth-based
algorithm, denoted as FDB, which approximates the MCD subset with a depth-
induced trimmed region. Specifically, we investigate FDB based on two
representative depth notions, projection depth and $L_{2}$ depth, and denote
the estimators as, $\texttt{FDB}_{\rm pro}$ and $\texttt{FDB}_{\rm L_{2}}$,
respectively. Four main advantages of the proposed algorithm are worth
mentioning: 1) Asymptotically, $\texttt{FDB}_{\rm pro}$ leads to a trimmed
region equivalent to the MCD subset for elliptically symmetric distributions.
2) Both $\texttt{FDB}_{\rm pro}$ and $\texttt{FDB}_{\rm L_{2}}$ achieve the
same level of robustness as the MCD estimator. 3) Empirically, FDB reveals
comparable or even better performance than the MCD estimator regarding
estimation accuracy. 4) Furthermore, the computational efficiency is
dramatically improved by using FDB, especially for high-dimensional cases.

The rest of the paper is organized as follows. Section 2 reviews the MCD
estimator and some related theoretical properties. Section 3 introduces the
idea of FDB estimators and demonstrates the theoretical equivalence between
the MCD subsets and the depth-trimmed regions. Section 4 investigates the
invariance, robustness, and computational complexity of the proposed FDB
estimators. In Section 5, we conduct extensive simulation studies to assess
the performance of the FDB algorithm and compare it with the existing ones
regarding estimation accuracy and computational efficiency. Section 6 applies
the proposed methods to several real applications through typical multivariate
analysis tasks, including principal component analysis, linear discriminant
analysis, image denoise, and outlier detection. We end the paper with some
discussion in Section 7. Proofs of the theoretical results and additional
simulation results are provided in the Supplementary Material.

##  2 MCD Estimators

In this section, we review the theoretical property of the MCD estimator as
well as three widely utilized approximation algorithms. Let
$\bm{x}\in\mathbb{R}^{p}$ be a random variable from an elliptically symmetric
distribution, denoted as ${\rm ES}(f;\bm{\mu},\bm{\Sigma})$, whose density is
of the form

$g(x)=c|\bm{\Sigma}|^{-1/2}f\left((\bm{x}-\bm{\mu})^{T}\bm{\Sigma}^{-1}(\bm{x}-\bm{\mu})\right),$

where $\bm{\Sigma}$ is a symmetric positive definite matrix, and the function
$f:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is assumed to be non-increasing so
that $g(x)$ is unimodal.

Considering random samples $\bm{x}_{1},\dots,\bm{x}_{n}$ independently
generated from the above distribution, MCD aims to solve the following
optimization problem

$\hat{H}_{\alpha_{n},{\rm
MCD}}=\underset{H\in\mathcal{H}_{h}}{\operatorname{argmin}}\left(\operatorname{det}\mathbf{\hat{\Sigma}}(\bm{x}_{H})\right),$
(1)

where $H$ is an index set of $h$ observations (with
$\lfloor(n+p+1)/2\rfloor\leqslant h\leqslant n$, where $\lfloor a\rfloor$
means the largest integer smaller than or equal to $a$), $\alpha_{n}=h/n$ and
$\mathcal{H}_{h}$ is the collection of all such sets. Observations with
corresponding indices $\hat{H}_{\alpha_{n},{\rm MCD}}$ constitute the final
MCD subset, termed $\hat{E}_{\alpha_{n},{\rm MCD}}$,. Define
$\Delta(A,B)=A\cup B-A\cap B$ as the difference between sets $A$ and $B$. The
convergence property of ${\hat{E}}_{\alpha_{n}{\rm MCD}}$ is revisited in
Lemma 1, with the proof provided in Section S1.1 of the Supplementary
Material.

###### Lemma 1

Assume that the random samples $\bm{x}_{i}\sim
ES(f;\bm{\mu},\bm{\Sigma}),i=1,\dots,n$. Then for $\alpha>0$, we have

$\mathbb{P}\left(\Delta\left(\hat{E}_{\alpha_{n},\text{\rm
MCD}},E_{\alpha}\right)\right)\rightarrow 0$

for any sequence $\alpha_{n}\to\alpha$ as $n\to\infty$, where
$E_{\alpha}=\left\\{\bm{x}\in\mathbb{R}^{p}\mid(\bm{x}-\bm{\mu})^{T}\bm{\Sigma}^{-1}(\bm{x}-\bm{\mu})\leq
r_{\alpha}^{2}\right\\}$ with $\mathbb{P}(E_{\alpha})=\alpha$.

Given an $n\times p$ data matrix
$\bm{X}=(\bm{x}_{1},\dots,\bm{x}_{n})^{\prime}$, with its estimated center
$\hat{\bm{\mu}}$ and scatter matrix $\hat{\bm{\Sigma}}$, we denote with
$\mathcal{D}\left(\bm{x}_{i};\hat{\bm{\mu}},\hat{\bm{\Sigma}}\right)=\sqrt{\left(\bm{x}_{i}-\hat{\bm{\mu}}\right)^{T}\hat{\bm{\Sigma}}^{-1}\left(\bm{x}_{i}-\hat{\bm{\mu}}\right)}$
the Mahalanobis distance of $\bm{x}_{i}$. The C-step described in Algorithm 1
is crucial for MCD-type algorithms.

Algorithm 1 The C-step

Input: initial subset $H_{old}$ or the estimates
($\hat{\bm{\mu}}_{\text{old}},\hat{\bm{\Sigma}}_{\text{old}}$), subset size
$h$.
Output: $H_{new}$ or
($\hat{\bm{\mu}}_{\text{new}},\hat{\bm{\Sigma}}_{\text{new}}$)

1:Compute the distances $d_{i,{\rm
old}}=\mathcal{D}\left(\bm{x}_{i};\hat{\bm{\mu}}_{\text{old}},\hat{\bm{\Sigma}}_{\text{old}}\right)$
for $i=1,\ldots,n$.

2:Sort these distances, yielding a permutation $\pi$ for which
$d_{\pi(1),\text{old}}\leqslant d_{\pi(2),\text{old}}\leqslant\ldots\leqslant
d_{\pi(n),\text{old}}$, and set $H_{new}=\\{\pi(1),\pi(2),\ldots,\pi(h)\\}$.

3:Update $\hat{\bm{\mu}}$ and $\hat{\bm{\Sigma}}$ as
$\hat{\bm{\mu}}_{\text{new}}=\frac{1}{h}\sum_{i\in H_{new}}\bm{x}_{i}\quad{\rm
and}\quad\hat{\bm{\Sigma}}_{\text{new}}=\frac{1}{h-1}\sum_{i\in
H_{new}}\left(\bm{x}_{i}-\hat{\bm{\mu}}_{\text{new}}\right)\left(\bm{x}_{i}-\hat{\bm{\mu}}_{\text{new}}\right)^{T}.$

Rousseeuw and Driessen (1999) proposed the first computationally feasible
algorithm, termed FASTMCD, for approximating the MCD subset. Specifically,
they randomly constructed several initial subsets and applied two C-steps for
each subset, yielding the ten subsets with the lowest determinant. Then, they
took C-step iteratively for these ten subsets until the determinant sequence
converged and eventually chose the MCD subset as the one leading to the
smallest determinant. Given this, the computational efficiency of FASTMCD is
thus roughly proportional to the number of the initial subsets. Hubert et al.
(2012) proposed an alternative algorithm DetMCD, which replaces random initial
subsets (of which there could be many) in FASTMCD with six well-designed
deterministic estimators of $(\bm{\mu},\bm{\Sigma})$, and also involves the
C-step in a similar way. Denote the estimates
$\left(\hat{\bm{\mu}},\hat{\bm{\Sigma}}\right)$ as the location and scatter
matrix estimates of the $h$-subset for which the determinant of the sample
covariance matrix is as small as possible. Further, an additional reweighted
step is employed in both algorithms to improve the efficiency of the
estimators. To be more specific, the estimators are renewed as trimmed
estimates for location and scatter,

$\hat{\bm{\mu}}_{\rm
re}=\frac{1}{\sum_{i=1}^{n}W_{i}}\sum_{i=1}^{n}W_{i}\bm{x}_{i}\quad{\rm
and}\quad\hat{\bm{\Sigma}}_{\rm
re}=\frac{1}{\sum_{i=1}^{n}W_{i}-1}\sum_{i=1}^{n}W_{i}\left(\bm{x}_{i}-\hat{\bm{\mu}}_{re}\right)\left(\bm{x}_{i}-\hat{\bm{\mu}}_{re}\right)^{T},$
(2)

where $W_{i}=1$ when
$\mathcal{D}\left(\bm{x}_{i};\hat{\bm{\mu}},c_{0}\hat{\bm{\Sigma}}\right)\leq\sqrt{\chi_{p,0.975}^{2}}$
and 0 otherwise, $\chi_{p,\alpha}^{2}$ is the $\alpha$-quantile of the
$\chi_{p}^{2}$ distribution and
$c_{0}=\operatorname{med}_{i}\mathcal{D}^{2}\left(\bm{x}_{i},\hat{\bm{\mu}},\hat{\bm{\Sigma}}\right)/{\chi_{p,0.5}^{2}}$.

##  3 A Depth-based Alternative

The idea behind the premier step of MCD-type algorithms is to construct
outlier-free subsets as the initial values for the C-step. This motivates us
to approach such a purpose using statistical depth, a popular tool for robust
multivariate data analysis. More importantly, we find the equivalence between
the eventual MCD subset and the depth-trimmed region, which avoids the
implementation of the iterative C-steps and hence improves the computational
efficiency dramatically.

In general, a statistical depth notion is a function
$D:(\bm{x},P)\mapsto[0,1]$, for $\bm{x}\in\mathbb{R}^{p}$ and $P$ from some
class $\mathcal{P}$ of $p$-variate probability distributions, that provides a
center-outward order for a collection of data. Taking into account the
robustness as well as the computational efficiency, to be discussed later, we
specifically consider the following two depth notions for the proposed method.

Projection depth (Zuo and Serfling, 2000a):

$D_{\operatorname{Proj}}(\bm{x};P)=\left(1+\sup_{\|\bm{u}^{\prime}\|=1}\frac{\left|\bm{u}^{\prime}\bm{x}-\operatorname{med}\left(\bm{u}^{\prime}\bm{y}\right)\right|}{\operatorname{MAD}\left(\bm{u}^{\prime}\bm{y}\right)}\right)^{-1},$

where $P$ is the distribution of $\bm{y}$, $\operatorname{med}(V)$ denotes the
median of a univariate random variable $V$, and $\operatorname{MAD}(V)=$
$\operatorname{med}(|V-\operatorname{med}(V)|)$ its median absolute deviation
from the median. Practically, one may choose a finite number of random
directions to approximate the projection depth values.

$\mathbb{L}_{2}$ depth (Zuo and Serfling, 2000a):

$D_{\mathbb{L}_{\mathrm{2}}}(\bm{x};P)=\left(1+E\left[\|\bm{y}-\bm{x}\|_{2}\right]\right)^{-1},$

where $\|\cdot\|_{2}$ is the $L^{2}$ norm.

Consider random samples $\bm{x}_{1},\dots,\bm{x}_{n}$ independently generated
from the elliptically symmetric (ES) distribution. For a given depth function
$D(\cdot;{\rm ES})$ and for $\alpha>0$, we call

$E_{\alpha,\text{depth}}=\left\\{\bm{x}\in\mathbb{R}^{p}\mid D(\bm{x};{\rm
ES})\geq D_{\alpha}\right\\},$

and the sample version is

$\hat{E}_{\alpha_{n},\mathrm{depth}}=\left\\{x\in\mathbb{R}^{p}\mid\hat{D}_{n}(\bm{x};{\rm
ES})\geq D_{\alpha_{n}}\right\\},$

the corresponding $\alpha$-trimmed region with
$\mathbb{P}(E_{\alpha,\text{depth}})=\alpha$.

######  Lemma 2 (Zuo and Serfling (2000b))

Assume that the random samples $\bm{x}_{i}\sim{\rm
ES}(f;\bm{\mu},\bm{\Sigma})$, $i=1,\dots,n$. Then for the projection depth,
the depth trimmed region (subset) $E_{\alpha_{n},\text{\rm depth}}$ satisfies

$\mathbb{P}\left(\Delta\left(\hat{E}_{\alpha_{n},\text{\rm
depth}},E_{\alpha}\right)\right)\rightarrow 0\text{ as }n\rightarrow\infty,$

for any sequence $\alpha_{n}\to\alpha>0$ as $n\to\infty$.

Other depth notions could also lead to the same conclusion except for the
projection depth. We omit those notions since they are either not robust
enough or computationally demanding. Also note that the result does not
necessarily hold for the $L_{2}$ depth, although $\texttt{FDB}_{\rm L_{2}}$
indeed provides satisfactory results in the simulation. Combining Lemmas 1 and
2, it is straightforward that the two subsets are asymptotic equivalent. We
state this result formally in the following theorem.

###### Theorem 1

Assume that the random samples $\bm{x}_{i}\sim{\rm
ES}(f;\bm{\mu},\bm{\Sigma})$, $i=1,\dots,n$. Under the conditions of Lemmas 1
and 2, we have

$\mathbb{P}(\Delta(E_{\alpha_{n},\text{\rm depth}},E_{\alpha_{n},\text{\rm
MCD}}))\to 0,\text{ as }n\rightarrow\infty.$

Proof of Theorem 1 is provided in Section S1.2 of the Supplementary Material.
Motivated by the above result, we propose to approximate the eventual MCD
subset with the depth-trimmed region and avoid the iterative implementation of
the C-step. We provide two toy examples in Figure 1 to illustrate such a
coincidence. Specifically, we generate data from bivariate normal
distributions with unit variance and correlation coefficients 0 and 0.5 for
Figure 1 and 1, respectively. We consider $n=4000$ and $h=3000$. In both
cases, the two subsets match quite well, such that the proportions of the
common elements are no less than 97%, which well supports the high
effectiveness of the proposed method.

Figure 1: The subsets induced by the projection depth and DetMCD. The black
points represent the samples that are outside both subsets; blue points
indicate samples in DetMCD subset only; red points indicate samples in depth-
based subset only; purple points are their intersection. Plots (a) and (b) are
for independent and correlated data, respectively.

For data of low dimensions, both MCD subsets and depth-based trimmed regions
can be computed efficiently, and the result matches well, as shown in Figure
1. However, for high-dimensional data, the MCD algorithms are severely
challenged by the cubically increased computational complexity, and hence the
approximation will be less efficient. To alleviate this challenge, we consider
replacing the MCD subsets with the trimmed regions induced by some
computationally efficient depth notions. By doing so, we may not only reduce
the computational time significantly but also attain comparable or better
robust estimation for both the location and scatter matrix, especially in
high-dimensional cases. In what follows, we introduce the FDB algorithm.

Algorithm 2 FDB

Input: $\bm{x_{i}},i=1,\dots,n$, subset size $h$, selected depth notion.
Output: $\hat{\bm{\mu}}_{\text{FDB}}$, $\hat{\bm{\Sigma}}_{\text{FDB}}$.

1:Calculate the depth value for each observation $\bm{x}_{i}$, denoted as
$\text{GD}(i)$.

2:Sort these values, yielding a permutation $\pi$ of $1,\dots,n$, for which
$\text{GD}(\pi(1))\geq,\dots,\geq\text{GD}(\pi(n))$, and a set
$H_{sub}=\\{\pi(1),\dots,\pi(h)\\}$.

3:Get the location and the scatter matrix estimates as  $\hat{\bm{\mu}}_{\rm
raw}=\frac{1}{h}\sum_{i\in H_{sub}}\bm{x}_{i}\quad{\rm
and}\quad\hat{\bm{\Sigma}}_{\rm raw}=\frac{c_{1}}{h}\sum_{i\in
H_{sub}}\left(\bm{x}_{i}-\bm{\mu}_{\rm
raw}\right)\left(\bm{x}_{i}-\bm{\mu}_{\rm raw}\right)^{T},$ where
$c_{1}=\underset{i}{\operatorname{med}}\mathcal{D}^{2}\left(\bm{x}_{i},\hat{\bm{\mu}}_{\rm
raw},\hat{\bm{\Sigma}}_{0}\right)/{\chi_{p,0.5}^{2}}$ with
$\hat{\bm{\Sigma}}_{0}=\frac{1}{h}\sum_{i\in
H_{sub}}\left(\bm{x}_{i}-\bm{\mu}_{\rm
raw}\right)\left(\bm{x}_{i}-\bm{\mu}_{\rm raw}\right)^{T}$.

4:Apply the reweighted step (2) to the raw estimates, yielding the final FDB
estimates, $\hat{\bm{\mu}}_{\mathrm{FDB}}$ and
$\hat{\bm{\Sigma}}_{\mathrm{FDB}}$.

Algorithm 2 considers the case of $h>p$, which is the condition to guarantee
the invertibility of estimated matrix (Rousseeuw and Van Zomeren, 1990). All
algorithms for the original MCD require that the dimension $p$ be lower than
$h$ to obtain an invertible covariance matrix. It is recommended that $n>5p$
in practice (Rousseeuw and Driessen, 1999).

According to Lemmas 1 and 2, for data from an elliptically symmetric
distribution, MCD and FDB algorithms both approximate the optimal subset,
though from different perspectives. MCD approaches the solution by combining
well selected (or random) initial subset and the iterative implementation of
the C-step, which could be computationally demanding for high dimensional
data. In contrast, FDB relies on ordering the data from the center outward,
and hence its computational complexity is mainly determined by the cost of
assigning depth values to each sample.

The idea of incorporating depth (outlyingness) to construct MCD estimators has
been considered in the literature. For example, the Stahel–Donoho outlyingness
(Donoho, 1982), equivalent to the projection depth, is applied to determine an
$h$-subset consisting of the $h$ points with the lowest outlyingness, and the
corresponding sample mean and covariance matrix are used as one initial value
for the C-step (Hubert et al., 2005a; Schreurs et al., 2021). Debruyne and
Hubert (2009) studied the influence function and asymptotic relative
efficiency of the estimators obtained directly based on such a subset (without
the reweighted step). For the first time, we establish the equivalence of the
two subsets, which indicates that the depth-based subset is a reasonable
approximation to the optimal subset rather than just one option of the initial
value for the C-step.

##  4 Properties of FDB

This section focuses on the properties of the FDB estimators. Specifically, We
discuss three types of properties, that are of main interest for such methods
(Maronna and Zamar, 2002; Hubert et al., 2012), invariance, robustness and
computational complexity. We show that the proposed estimators are quite
satisfactory in these aspects.

###  4.1 Invariant Properties

Affine equivariance makes the analysis independent of any affine
transformation of the data. For any nonsingular $p\times p$ matrix $\bm{A}$
and $p\times 1$ vector $\bm{v}$, the estimators $\hat{\bm{\mu}}$ and
$\hat{\bm{\Sigma}}$ are affine equivariant if they satisfy

$\hat{\bm{\mu}}\left(\bm{X}\bm{A}+\mathbf{1}_{n}\bm{v}^{T}\right)=\hat{\bm{\mu}}(\bm{X})\bm{A}+\bm{v}{\rm\leavevmode\nobreak\
and\leavevmode\nobreak\
}\hat{\bm{\Sigma}}\left(\bm{X}\bm{A}+\mathbf{1}_{n}\bm{v}^{T}\right)=\bm{A}^{T}\hat{\bm{\Sigma}}(\bm{X})\bm{A},$

where $\mathbf{1}_{n}=(1,1,\ldots,1)^{T}$. The projection depth has been shown
affine equivariant (Zuo and Serfling, 2000a; Zuo, 2006), that is the depth
value does not vary through affine transformation for any sample, and hence
the indexes of samples forming the trimmed region remain the same.
Consequently, $\texttt{FDB}_{\rm pro}$ is obviously affine equivariant. For
$\texttt{FDB}_{\rm L_{2}}$, a similar property holds for rigid transformation
(Mosler and Mozharovskyi, 2022), which is a bit more restrictive than the
affine transformation. For high dimensional situations, the affine
equivariance may be less important under nonstandard data contamination such
as componentwise outliers (Alqallaf et al., 2009).

Permutation invariance provides an effective way to guaranteeing the
robustness of analysis to the perturbation of observations. An estimator
$T(\cdot)$ is said to be permutation invariant if $T(\bm{P}\bm{X})=T(\bm{X})$
for any permutation matrix $\bm{P}$. Permutation invariance holds for both
projection depth and $L_{2}$ depth since they do not involve any random
subsets, and hence the depth values remain the same through permutation and so
will the $h$-subset.

###  4.2 Robustness

Robustness is the property of main interest when outlier contamination of the
data is suspected. As aforementioned, the MCD estimator is highly robust that
it achieves the highest possible asymptotic breakdown point, about 1/2, with
$h\approx n/2$. The robustness of FDB is determined by the property of the
employed depth notion. According to Zuo (2006) and Lopuhaa and Rousseeuw
(1991), the breakdown points of the trimmed regions induced by the projection
depth and $L_{2}$ depth are both 1/2, with $\alpha\approx 0.5$. That is,
$\texttt{FDB}_{\rm pro}$ and $\texttt{FDB}_{\rm L_{2}}$ both have a breakdown
point as high as that of the MCD estimator. Another indicator is the influence
function, which captures the local robustness of estimators. Zuo (2006) and
Niinimaa and Oja (1995) showed the influence functions of depth regions
induced by projection depth and $L_{2}$ depth are both bounded. Therefore,
$\texttt{FDB}_{\rm pro}$ and $\texttt{FDB}_{\rm L_{2}}$ are highly robust
locally as well as globally.

The robustness discussed above is from the theoretical aspect. Practically,
both FDB and MCD approximate the theoretically optimal subsets, and their
empirical performances do not necessarily match the theoretical results under
all the scenarios. This means that the subset selected by the two methods
under finite samples may still contain outliers. To show this point, we
provide an example in Figure 2, where we generate $n=4000$ samples from a
40-dimensional normal distribution with standard normal marginal distributions
and a correlation coefficient of 0.5. We consider two levels of contamination,
10% and 40%, and two types of outliers, point and cluster (see for details in
Section 5.1), respectively. For the first column with 10% outliers, we set
$h=\lfloor 75\%n\rfloor$; for the second column with 40% outliers, $h=\lfloor
50\%n\rfloor$. $\texttt{FDB}_{\rm pro}$ performs perfectly for the first three
cases and fails for the last scenario; in contrast, DetMCD is only
satisfactory for the first case but fails for the rest. This indicates that
for high-dimensional data, $\texttt{FDB}_{\rm pro}$ provides more reliable
approximations to the optimal subset.

(a) 10% cluster

(b) 40% cluster

(c) 10% point

(d) 40% point

Figure 2: The subsets induced by $\texttt{FDB}_{\rm pro}$ and DetMCD. Dots
denote normal sample and crosses denote outliers. Blue ones form the
intersection of two subsets; orange (purple) ones indicate the samples unique
to the $\texttt{FDB}_{\rm pro}$ (DetMCD) subset; blue ones are samples dropped
by both subsets.

###  4.3 Computational complexity

(a) $n=1000$ (b) $p=200$

Figure 3: The average computation time (seconds) for different settings over
20 replicates. (a): $\log(t)$ versus $\log(p)$ with $n=1000$; (b): $t$ versus
$n$ with $p=200$.

The computational complexity for finding the $h$-subset by MCDs is
$O(\psi(np^{2}+p^{3}))$. Specifically, for each C-step, it requires computing
the covariance matrix and Mahalanobis distances, with complexities $O(np^{2})$
and $O(p^{3}+np^{2})$ respectively, and $\psi$ depends on the number of
initial estimates and the times of C-step iteration. For FASTMCD, the number
of initial estimates defaults to 500; for DetMCD, it defaults to 6; and for
RT-DetMCD, it is further reduced to 2. However, these efforts only reduced the
value of $\phi$ but the order term still remains the same.

In contrast, the computational complexity of $\texttt{FDB}_{\rm pro}$ for
finding the subset is $O(knp)$ with $k$ as the number of projection
directions. According to our numerical experiments, The performance of
$\texttt{FDB}_{\rm pro}$ is quite stable when the number of random directions
is set around $1000$; see Figure S2 of the Supplement. Hence,
$\texttt{FDB}_{\rm pro}$ leads to a significant improvement over the MCD
estimators. We remark that it is possible to further reduce the number of
projection directions according to some elaborate generative algorithms
(Dyckerhoff et al., 2021). However, these algorithms may instead lengthen the
total computational time due to the tedious procedure for searching “better”
directions, and hence we stick to selecting the directions randomly. For the
case of ultra-high dimensional data, we suggest an adaptive rule,
$k=\max(1000,10p)$. As for $\texttt{FDB}_{\rm L_{2}}$, the computational
complexity is $O(n^{2}p)$, which scales linearly with the dimension of data.

To show the improvement, we provide some numerical results for computation
time in Figure 3. Notably, the speed of DetMCD is also influenced by the way
of constructing the initial estimates. Specifically, $Q_{n}$ (Rousseeuw and
Croux, 1993) is applied to construct initial estimates for DetMCD, which is
computationally demanding. To improve speed, Hubert et al. (2012) suggested
substituting $Q_{n}$ with the $\tau$-scale of Yohai and Zamar (1988) when
$n>1000$. We follow this suggestion by using the $Q_{n}$ estimator for DetMCD
in Figure 3(a), and the $\tau$-estimator in Figure 3(b), respectively. For
$\texttt{FDB}_{\rm pro}$, we let $k=1000$. All experiments are run using
R-package ddalpha for DetMCD on an Intel(R) Xeon(R) with 3.10GHz and 192 GB
memory processor.

FDBs show significant improvement over DetMCD under all the settings.
Specifically, in Figure 3(a), the line of DetMCD is steeper than those of the
other two methods, which matches well with different orders of dimension $p$
in their theoretical computational complexities aforementioned. In Figure
3(b), both DetMCD and $\texttt{FDB}_{\rm pro}$ reveal linear trends with the
increasing sample size, while $\texttt{FDB}_{\rm L_{2}}$ shows a quadratic
trend, though its computation time is the least when $n<4000$.

##  5 Simulations

We conduct extensive simulations with data from symmetric distributions to
assess the performance of our proposed algorithms, $\texttt{FDB}_{\rm pro}$
and $\texttt{FDB}_{L_{2}}$, and make a comparison with DetMCD (Hubert et al.,
2012). Besides, we also provide some exploration for the scenarios of
asymmetric distributions in Section S3.1 of the Supplement. To evaluate the
estimation results, we use the following five measures (the smaller the
better).

  * •

An error measure of the location estimator, given by
$e_{\mu}=\|\hat{\bm{\mu}}-\bm{\mu_{0}}\|$, where $\bm{\mu_{0}}$ denotes the
true sample mean.

  * •

An error measure of the scatter estimator, defined as the logarithm of the
condition number of $\hat{\mathbf{\Sigma}}\mathbf{\Sigma^{-1}}$,
$e_{\Sigma}=\log_{10}(\operatorname{cond}(\hat{\mathbf{\Sigma}}\mathbf{\Sigma}^{-1})).$

  * •

The mean squared error (MSE) of $\mathbf{\Sigma}$, ${\rm
MSE}=\frac{1}{Sp^{2}}\sum_{s=1}^{S}||\hat{\mathbf{\Sigma}}-\mathbf{\Sigma}||_{F}^{2}.$

  * •

The Kullback Leibler (KL) divergence between $\hat{\mathbf{\Sigma}}$ and
${\mathbf{\Sigma}}$, ${\rm
KL}\left(\hat{\mathbf{\Sigma}},\mathbf{\Sigma}\right)=\operatorname{trace}\left(\hat{\mathbf{\Sigma}}\mathbf{\Sigma}^{-1}\right)-\log\left(\operatorname{det}\left(\hat{\mathbf{\Sigma}}\mathbf{\Sigma}^{-1}\right)\right)-p,$,
which is identical the KL divergence between the two Gaussian distributions
with the same mean.

  * •

The computation time $t$ (in seconds) of the whole procedure, including the
optimal subset pursuit and the reweighted step.

###  5.1 Estimation performance

In this subsection, we generate the bulk of non-outlying samples as
$\bm{x}_{i}=\bm{G}\bm{y}_{i}$, where $\bm{y}_{i}$ are from
$N_{p}(\mathbf{0},\mathbf{I})$ and $\bm{G}$ is a $p\times p$ matrix with unit
diagonal elements and off-diagonal elements equal to 0.75. The number of
outliers is $m=\lfloor n\varepsilon\rfloor$ and $\varepsilon$ denotes the
level of contamination. Four contamination types are considered: point,
random, cluster, and radial outliers. Point outliers are obtained by
generating $\bm{y}_{i}\sim
N_{p}\left(r\bm{a}\sqrt{p},0.01^{2}\mathbf{I}\right)$, where $\bm{a}$ is a
unit vector generated orthogonal to $\bm{a_{0}}=(1,1,\ldots,1)^{T}$. Random
outliers are obtained by generating $\bm{y}_{i}\sim
N_{p}\left(\bm{\mu}_{ir},\mathbf{I}\right)$, where
$\bm{\mu}_{ir}=rp^{1/4}\bm{\nu}/\|\bm{\nu}\|$ with $\bm{\nu}$ a random vector
from $N_{p}(\bm{0},\mathbf{I})$. Cluster outliers are obtained by generating
$\bm{y}_{i}\sim N_{p}\left(rp^{-1/4}\bm{a_{0}},\mathbf{I}\right)$, where $r$
is a constant. Radial outliers are obtained by generating $\bm{y}_{i}\sim
N_{p}(\mathbf{0},5\mathbf{I})$. Except for the random outliers, the other
three types have been considered by Hubert et al. (2012).

Different contamination levels are considered, namely $\varepsilon=0\%,10\%$,
and $40\%$. Let $h$ be $\lfloor 0.75n\rfloor$ when $\varepsilon=0\%$ or
$10\%$, and $\lfloor 0.5n\rfloor$ when $\varepsilon=40\%$ for each method
under investigation. The number of directions for projection depth is set as
$k=\max(1000,10p)$ as suggested in Section 4.3. For each method, we compute
the reweighted location vectors $\hat{\bm{\mu}}_{\bm{X}}$ and the reweighted
scatter matrices $\hat{\bm{\Sigma}}_{\bm{X}}$. The corresponding estimators
for the data set $\bm{Y}$ are obtained as
$\hat{\bm{\mu}}_{\bm{Y}}=\bm{G}^{-1}\hat{\bm{\mu}}_{\bm{X}}$ and
$\hat{\bm{\Sigma}}_{\bm{Y}}=\bm{G}^{-1}\hat{\mathbf{\Sigma}}_{\bm{X}}\bm{G}^{-1}$,
which are compared to the true values using the aforementioned measures. In
this part, the true covariance matrix of $\bm{Y}$ is
${\mathbf{\Sigma}}=\mathbf{I}_{p}$.

(a) point

(b) random

(c) cluster

(d) point

(e) random

(f) cluster

Figure 4: The average $e_{\Sigma}$ for levels of abnormality ($r$), with
$n=400$ and $p=40$, the first row is for $\varepsilon=0.1$ and $\alpha=0.75$,
and the second row is for $\varepsilon=0.4$ and $\alpha=0.5$.

(a) point

(b) cluster

(c) point

(d) cluster

(e) point

(f) cluster

(g) point

(h) cluster

Figure 5: The average $e_{\Sigma}$ for different $p$, with $n=2000$, the first
row is for $\varepsilon=0.1$ and $\alpha=0.75$; the second row is for
$\varepsilon=0.4$ and $\alpha=0.5$; the first two columns are for $r=2$, and
the last two columns are for $r=20$, respectively. The cases of random
contamination is omitted since all methods perform similarly.

We first provide a full picture for the performance of each method by
conducting simulation studies under a broad range of settings. To be more
specific, we generated $1000$ data sets consisting of different types of
contamination under a broad range of $r$ and $p$, and computed the average
$e_{\Sigma}$ with the three methods. The results are illustrated in Figures 4
and 5. Other measures, $e_{\mu}$, MSE, and KL divergence, all reveal similar
patterns as $e_{\Sigma}$ and hence are omitted. As shown in Figure 4,
$\texttt{FDB}_{\rm pro}$ and $\texttt{FDB}_{\rm L_{2}}$ outperform DetMCD
across all values of $r$ for both random and cluster outliers under either low
or high contamination levels. For the case of point contamination,
$\texttt{FDB}_{\rm pro}$ still holds the upper hand, with $\texttt{FDB}_{\rm
L_{2}}$ slightly worse than DetMCD, especially for small $r$ values.

Figure 5 shows that $\texttt{FDB}s$ is among the best estimators for all the
settings. However, in general, DetMCD’s performance deviates more seriously
with the increase of dimension $p$, suggesting that it is less suitable for
high-dimensional data. Overall, $\texttt{FDB}_{\rm pro}$ provides the best
performance among the three options, that both $\texttt{FDB}_{\rm L_{2}}$ and
DetMCD produce large $e_{\Sigma}$ under point contamination since their
induced subsets may contain outliers, which is consistent with the results in
Figure 2. Between $\texttt{FDB}_{\rm L_{2}}$ and DetMCD, the later is better
for point contamination while the former is better (or even the best) for
cluster contamination.

Next, we provide more detailed numerical outputs for typical simulation
settings under study. Following Hubert et al. (2012), we consider three
options, A: $n=200$ and $p=5$, B: $n=400$ and $p=40$, and C: $n=2000$ and
$p=200$, representing low, moderate and high dimensions, respectively. Other
settings remain the same as those in Figure 4 except that $r$ is fixed at $5$
for point, random, and cluster outliers. We report the average measures over
1000 runs in Tables 1, 2 and 3, corresponding to 0% (clean data), 10% and 40%
contamination, respectively.

Table 1: Simulation results for clean data. | $e_{\mu}$ | $e_{\Sigma}$ | MSE | KL | t
---|---|---|---|---|---
A | DetMCD | 0.158 (0.050) | 0.238 (0.051) | 0.008 (0.003) | 0.119 (0.046) | 0.022
$\texttt{FDB}_{\rm L_{2}}$ | 0.157 (0.051) | 0.245 (0.054) | 0.008 (0.003) | 0.122 (0.049) | 0.003
$\texttt{FDB}_{\rm pro}$ | 0.157 (0.050) | 0.232 (0.049) | 0.007 (0.003) | 0.112 (0.044) | 0.011
B | DetMCD | 0.344 (0.040) | 0.615 (0.029) | 0.003 (0.000) | 2.771 (0.179) | 0.388
$\texttt{FDB}_{\rm L_{2}}$ | 0.329 (0.038) | 0.571 (0.025) | 0.003 (0.000) | 2.404 (0.135) | 0.014
$\texttt{FDB}_{\rm pro}$ | 0.331 (0.038) | 0.573 (0.025) | 0.003 (0.000) | 2.435 (0.137) | 0.039
C | DetMCD | 0.355 (0.016) | 0.656 (0.011) | 6e-4 (0.000) | 14.202 (0.249) | 3.876
$\texttt{FDB}_{\rm L_{2}}$ | 0.331 (0.016) | 0.596 (0.008) | 5e-4 (0.000) | 11.796 (0.166) | 1.043
$\texttt{FDB}_{\rm pro}$ | 0.334 (0.016) | 0.598 (0.008) | 5e-4 (0.000) | 11.876 (0.177) | 1.109
Table 2: Simulation results for data with $\varepsilon=10\%$ and $r=5$. | Point | Random | Cluster | Radial
---|---|---|---|---
| Det | $\texttt{FDB}_{\rm L_{2}}$ | $\texttt{FDB}_{\rm pro}$ | Det | $\texttt{FDB}_{\rm L_{2}}$ | $\texttt{FDB}_{\rm pro}$ | Det | $\texttt{FDB}_{\rm L_{2}}$ | $\texttt{FDB}_{\rm pro}$ | Det | $\texttt{FDB}_{\rm L_{2}}$ | $\texttt{FDB}_{\rm pro}$
A | $e_{\mu}$ | 0.166 | 1.188 | 0.165 | 0.164 | 0.166 | 0.164 | 0.163 | 0.163 | 0.163 | 0.165 | 0.171 | 0.163
(0.054) | (0.071) | (0.054) | (0.054) | (0.056) | (0.054) | (0.053) | (0.053) | (0.053) | (0.051) | (0.053) | (0.051)
$e_{\Sigma}$ | 0.236 | 1.333 | 0.236 | 0.228 | 0.236 | 0.230 | 0.228 | 0.227 | 0.229 | 0.233 | 0.251 | 0.233
(0.052) | (0.070) | (0.052) | (0.045) | (0.047) | (0.046) | (0.045) | (0.045) | (0.045) | (0.049) | (0.054) | (0.048)
MSE | 0.008 | 5.526 | 0.008 | 0.008 | 0.008 | 0.008 | 0.007 | 0.007 | 0.007 | 0.008 | 0.009 | 0.008
(0.003) | (0.268) | (0.003) | (0.003) | (0.003) | (0.003) | (0.03) | (0.003) | (0.003) | (0.003) | (0.003) | (0.003)
KL | 0.107 | 9.375 | 0.112 | 0.101 | 0.112 | 0.105 | 0.101 | 0.101 | 0.102 | 0.103 | 0.117 | 0.104
(0.041) | (0.314) | (0.043) | (0.038) | (0.042) | (0.039) | (0.038) | (0.038) | (0.038) | (0.038) | (0.046) | (0.037)
B | $e_{\mu}$ | 2.403 | 3.456 | 0.339 | 0.347 | 0.339 | 0.340 | 0.349 | 0.340 | 0.341 | 0.347 | 0.337 | 0.338
(1.597) | (0.67) | (0.037) | (0.040) | (0.038) | (0.039) | (0.037) | (0.036) | (0.038) | (0.039) | (0.037) | (0.037)
$e_{\Sigma}$ | 1.787 | 2.420 | 0.591 | 0.622 | 0.595 | 0.596 | 0.616 | 0.588 | 0.592 | 0.616 | 0.593 | 0.594
(0.911) | (0.028) | (0.026) | (0.031) | (0.026) | (0.027) | (0.026) | (0.025) | (0.024) | (0.027) | (0.026) | (0.026)
MSE | 4.049 | 5.912 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003
(3.145) | (0.147) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000)
KL | 64.051 | 95.728 | 2.591 | 2.789 | 2.591 | 2.597 | 2.781 | 2.558 | 2.591 | 2.801 | 2.574 | 2.587
(47.593) | (1.304) | (0.131) | (0.172) | (0.141) | (0.138) | (0.168) | (0.144) | (0.140) | (0.162) | (0.147) | (0.142)
C | $e_{\mu}$ | 7.824 | 8.056 | 0.348 | 0.355 | 0.340 | 0.341 | 0.359 | 0.344 | 0.346 | 0.362 | 0.345 | 0.347
(2.507) | (0.082) | (0.018) | (0.015) | (0.015) | (0.015) | (0.017) | (0.016) | (0.017) | (0.018) | (0.018) | (0.018)
$e_{\Sigma}$ | 2.965 | 3.151 | 0.618 | 0.654 | 0.617 | 0.618 | 0.654 | 0.617 | 0.618 | 0.653 | 0.613 | 0.616
(0.773) | (0.012) | (0.009) | (0.010) | (0.009) | (0.009) | (0.010) | (0.008) | (0.008) | (0.012) | (0.010) | (0.009)
MSE | 6.488 | 6.366 | 6e-04 | 6e-04 | 6e-04 | 6e-04 | 6e-04 | 6e-04 | 6e-04 | e-04 | e-04 | e-04
(2.185) | (0.103) | (8e-06) | (7e-06) | (7e-06) | (7e-06) | (8e-06) | (8e-06) | (8e-06) | (8e-06) | (8e-06) | (8e-06)
KL | 495.2 | 513.9 | 12.61 | 14.08 | 12.648 | 12.92 | 14.02 | 12.57 | 12.65 | 14.06 | 12.58 | 12.65
(161.4) | (4.320) | (0.174) | (0.167) | (0.151) | (0.150) | (0.195) | (0.141) | (0.142) | 0.174 | 0.156 | 0.142
Table 3: Simulation results for data with $\varepsilon=40\%$ and $r=5$. | Point | Random | Cluster | Radial
---|---|---|---|---
| Det | $\texttt{FDB}_{\rm L_{2}}$ | $\texttt{FDB}_{\rm pro}$ | Det | $\texttt{FDB}_{\rm L_{2}}$ | $\texttt{FDB}_{\rm pro}$ | Det | $\texttt{FDB}_{\rm L_{2}}$ | $\texttt{FDB}_{\rm pro}$ | Det | $\texttt{FDB}_{\rm L_{2}}$ | $\texttt{FDB}_{\rm pro}$
A | $e_{\mu}$ | 7.226 | 6.946 | 6.335 | 0.193 | 0.280 | 0.191 | 0.604 | 0.179 | 0.240 | 0.216 | 0.233 | 0.215
(0.982) | (0.354) | (0.376) | (0.066) | (0.100) | (0.066) | (0.215) | (0.061) | (0.106) | (0.072) | (0.071) | (0.072)
$e_{\Sigma}$ | 2.989 | 2.981 | 2.563 | 0.279 | 0.623 | 0.278 | 0.421 | 0.263 | 0.284 | 0.299 | 0.447 | 0.291
(0.290) | (0.208) | (0.184) | (0.062) | (0.126) | (0.062) | (0.386) | (0.058) | (0.170) | (0.058) | (0.081) | (0.053)
MSE | 31.34 | 33.45 | 36.43 | 0.011 | 0.474 | 0.011 | 1.037 | 0.009 | 0.071 | 0.041 | 0.120 | 0.045
(3.592) | (2.941) | (2.329) | (0.006) | (0.360) | (0.005) | (2.017) | (0.004) | (0.111) | (0.019) | (0.054) | (0.020)
KL | 29.66 | 29.46 | 29.91 | 0.142 | 2.363 | 0.141 | 1.673 | 0.130 | 0.349 | 0.365 | 0.963 | 0.400
(0.918) | (0.894) | (0.795) | (0.060) | (1.453) | (0.058) | (3.523) | (0.053) | (0.483) | (0.142) | (0.324) | (0.152)
B | $e_{\mu}$ | 22.98 | 25.30 | 25.10 | 0.408 | 0.408 | 0.407 | 4.796 | 0.411 | 0.409 | 0.412 | 0.413 | 0.410
(0.079) | (0.041) | (1.138) | (0.047) | (0.048) | (0.048) | (0.371) | (0.045) | (0.045) | (0.044) | (0.046) | (0.045)
$e_{\Sigma}$ | 4.373 | 6.161 | 6.003 | 0.732 | 0.729 | 0.724 | 2.000 | 0.718 | 0.720 | 0.737 | 0.739 | 0.726
(0.112) | (0.151) | (0.470) | (0.033) | (0.036) | (0.035) | (0.093) | ( 0.031) | (0.032) | (0.037) | (0.036) | (0.033)
MSE | 24.51 | 15.98 | 16.55 | 0.004 | 0.004 | 0.004 | 0.875 | 0.004 | 0.004 | 0.004 | 0.004 | 0.004
(0.568) | (0.366) | (3.196) | (0.000) | (0.000) | (0.000) | (0.074) | (0.000) | (0.000) | (0.000) | (0.000) | (0.000)
KL | 242.1 | 229.1 | 231.1 | 3.868 | 3.815 | 3.796 | 36.654 | 3.731 | 3.747 | 3.877 | 3.890 | 3.799
(2.489) | (2.716) | (5.105) | (0.217) | (0.217) | (0.216) | (2.466) | (0.201) | (0.210) | (0.218) | (0.236) | (0.204)
C | $e_{\mu}$ | 51.43 | 56.56 | 55.64 | 0.420 | 0.409 | 0.413 | 6.345 | 0.418 | 0.417 | 0.421 | 0.410 | 0.415
(0.016) | (0.012) | (3.242) | (0.020) | (0.020) | (0.021) | (0.161) | (0.019) | (0.020) | (0.021) | (0.020) | (0.021)
$e_{\Sigma}$ | 5.061 | 6.996 | 6.756 | 0.769 | 0.745 | 0.753 | 2.317 | 0.757 | 0.758 | 0.769 | 0.747 | 0.760
(0.039) | (0.027) | (0.674) | (0.012) | (0.011) | (0.011) | (0.009) | (0.012) | (0.011) | (0.011) | (0.009) | (0.011)
MSE | 24.57 | 16.01 | 17.33 | 9e-04 | 8e-04 | 4e-04 | 0.155 | 8e-04 | 8e-04 | 9e-04 | 8e-04 | 8e-04
(0.085) | (0.055) | (4.426) | (1e-05) | (1e-05) | (1e-05) | (0.004) | (1e-05) | (1e-05) | (1e-05) | (9e-06) | (1e-05)
KL | 1204 | 1151 | 1140 | 19.20 | 18.09 | 18.45 | 88.81 | 18.63 | 18.66 | 19.14 | 18.12 | 18.70
(2.094) | (2.722) | (31.32) | (0.238) | (0.199) | (0.224) | (0.975) | (0.233) | (0.237) | (0.240) | (0.198) | (0.222)

For clean data (Table 1), the three estimators are comparable for low-
dimensional cases; FDBs achieve slightly smaller values for $e_{\mu}$,
$e_{\Sigma}$ and KL when the data dimension is moderate or high. More
importantly, the running time of DetMCD is reduced in all settings, with the
relative computational efficiency, defined as $t_{\rm MCD}/t_{\rm FDB}$,
ranging between 2 and 10 for $\texttt{FDB}_{\rm pro}$, and between 3 and 27
for $\texttt{FDB}_{\rm L_{2}}$. Such a comparison of computation time holds
similarly when contamination presents, and hence is omitted in the remaining
tables.

When the amount of contamination is 10% (Table 2), the performances of three
methods are all satisfying under settings with random, cluster or radial
contamination, and the comparison is similar to that for clean data. For point
contamination, $\texttt{FDB}_{\rm L_{2}}$ performs the worse for each of the
three settings, that it generates the largest values for the four measures of
estimation accuracy; DetMCD gets problematic for moderate and high dimensional
data (settings B and C), which indicates its deficiency for such cases; in
contrast, $\texttt{FDB}_{\rm pro}$ remains robust under all three settings and
provides values of measures quite close to corresponding ones for the clean
data in Table 1.

When the amount of contamination increases to 40% (Table 3), the three methods
all work well and are comparable under settings with random or radial
contamination. However, none of them is satisfactory when the contamination
presents as point outliers, and this weak performance was also observed for
both DetMCD and FASTMCD under similar settings in Hubert et al. (2012). For
cluster contamination, DetMCD leads to larger estimation errors in each cell,
especially when the dimension of data is moderate or high; however, both
$\texttt{FDB}_{\rm L_{2}}$ and $\texttt{FDB}_{\rm pro}$ instead remains very
stable across all three settings. Additional results for $r=2$ and $20$ are
provided in Tables S1–S4 of the Supplementary Material, from which we may draw
similar conclusions for the comparison of the three methods.

(a) point

(b) random

(c) cluster

(d) radial

Figure 6: The average $e_{\Sigma}$ for different $\alpha$, with $n=400$,
$p=40$, $r=100$, and $\varepsilon=0.4$.

###  5.2 Robustness assessment

In addition, we assess the tolerance of each method to the core-set size
$h=\lfloor\alpha n\rfloor$. Specifically, we generate data from Setting B with
40% contamination, which is high enough for most practical implementations,
and we set $r=100$ for point, random and cluster outliers. Meanwhile, we
consider $\alpha$ ranging from 0.5 to 0.6, which is the highest value that
possibly produces a clean subset. The average $e_{\rm\Sigma}$ from 1000
replications are reported in Figure 6. $\texttt{FDB}_{\rm pro}$ remain robust
for different values of $\alpha$ under all the investigated settings;
$\texttt{FDB}_{\rm L_{2}}$ is satisfactory for random, cluster and radial
contamination, while its estimation error grows substantially with $\alpha$
for point contamination; DetMCD generates table estimation when $\alpha=0.5$
and $0.55$ but becomes deficient when $\alpha$ raises up to 0.6 in each plot
of Figure 6. The other three measures show similar patterns as illustrated in
Figure S1 of the Supplementary Material. In conclusion, $\texttt{FDB}_{\rm
pro}$ reaches the strongest tolerance to the core-set size, when the
proportion of outliers in the data is very high.

To sum it up, FDBs improves the computational efficiency significantly, which
is the main motivation of this work, and $\texttt{FDB}_{\rm L_{2}}$ generally
achieves the highest computational efficiency. Besides, we surprisingly find
that $\texttt{FDB}_{\rm pro}$ shows superiority in terms of both estimation
accuracy and robustness, especially in high-dimensional cases. One may safely
substitute DetMCD with $\texttt{FDB}_{\rm pro}$ for practical implementation.

##  6 Real data examples

In this section, we apply the FDB methods to four real datasets of various
dimensions and sizes. The resultant robust multivariate location and scatter
estimates are evaluated via some typical tasks in multivariate analysis such
as outlier detection, linear discriminant analysis (LDA), principal component
analysis (PCA), and image denoise. The same three methods from Section 5 are
utilized. The computations are performed in R (R Core Team 2021) on a laptop
with a 10-core and 32GB memory processor.

###  6.1 Robust PCA for forged bank notes data

The first dataset is the forged Swiss bank notes data (Milo, 1990), which is
also used in Hubert et al. (2012). The data are of size $n=100$ and dimension
$p=6$, denoted as $\bm{X}\in\mathbb{R}^{n\times p}$. Since this dataset
includes outliers and highly correlated variables (Rousseeuw et al., 2004a;
Willems et al., 2009), we employ the proposed algorithms and DetMCD to get
robust estimation $(\hat{\bm{\mu}},\hat{\bm{\Sigma}})$ first and conduct
robust PCA based on these estimates. The classical PCA obtained by sample
location and scatter is also demonstrated for comparison. We use the first two
principal components $\bm{P}\in\mathbb{R}^{p\times 2}$, which explain over
$80\%$ of the total variance. The projections of data on the $2$-dimensional
PCA subspace,
$\bm{T}=\\{t_{ik}\\}=(\bm{X}-\mathbf{1}_{n}\hat{\bm{\mu}}^{T})\bm{P}$, are
shown in Fig. 7(e)–7(h).

(a) $\texttt{FDB}_{\rm L_{2}}$

(b) $\texttt{FDB}_{\rm pro}$

(c) DetMCD

(d) Sample

(e) $\texttt{FDB}_{\rm L_{2}}$

(f) $\texttt{FDB}_{\rm pro}$

(g) DetMCD

(h) Sample

Figure 7: Diagnostic plots of the forged bank notes using robust PCAs based on
(a) $\texttt{FDB}_{\rm L_{2}}$, (b) $\texttt{FDB}_{\rm pro}$ and (c) DetMCD
and classical PCA based on (d) Sample. Projection of data on 2-dimensional
subspace and diagnostic results obtained by (e) $\texttt{FDB}_{\rm L_{2}}$,
(f) $\texttt{FDB}_{\rm pro}$, and (g) DetMCD and (h) classical PCA. Purple
points in (c) and (g) represent samples that may be misclassified by using
DetMCD.

The score distance (SD) and orthogonal distance (OD) represent the robust
distance of samples in the two-dimensional PCA subspace and the orthogonal
distance of samples to the PCA subspace, respectively. For each sample
$\bm{x}_{i}$, we have $\mathbf{SD}_{i}=\sum_{k=1}^{2}t_{ik}^{2}/\lambda_{k}$,
where $\lambda_{1}$ and $\lambda_{2}$ are the first two eigenvalues of
$\hat{\bm{\Sigma}}$, and $\mathbf{OD}_{i}=\sum_{j=1}^{p}e_{i,j}^{2}$, where
$\bm{X}-\bm{TP}^{T}=\\{e_{i,j}\\}$. Figure 7(a)–7(d) illustrate diagnostic
plots. By two cutoff values in each diagnostic plot, we categorize samples
into four types and assign different colors to them in Fig. 7(e)–7(h). Regular
samples gather in the bottom-left region of diagnostic plots with both score
distances and orthogonal distances relatively small, which form the main body
of the data cloud. Good leverage samples are close to the PCA subspace but far
from the regular samples, e.g., the samples $13$ and $23$ in the bottom-right
region of Fig. 7(a)–7(c), and the green points in Fig. 7(e)–Fig. 7(g).
Orthogonal outliers are far from the PCA subspace but not distinguishable by
only observing their projections. With larger orthogonal distances but smaller
score distances, they locate on the top-left region of diagonal plots, e.g.,
samples $11$, $62$ and $67$ in Fig. 7(a)–7(c), and red points in Fig.
7(e)–Fig. 7(g). For bad leverage samples, both score and orthogonal distances
are large. They lie on the top-right region of diagonal plots and are
represented as orange points in projection plots. In general, the three robust
methods lead to comparable analysis results and all significantly improve the
performance of the classical PCA, which agrees with Theorem 1. However, DetMCD
may identify some regular points as outliers according to its specific cutoff
value for the orthogonal distance; see the purple points in the third column
of Figure 7.

###  6.2 Outlier detection for phoneme data

In the second example, we detect outliers for a phoneme dataset, which comes
from a speech recognition database TIMIT and has been discussed in Hastie et
al. (2009). The data includes $1050$ speech frames, $1000$ of which are “ao”
and $50$ of which are “iy”. Each data frame has been transformed to be a log-
periodogram of length $256$. First, we reduce the dimensions by smoothing
splines. For each sample, we replace the original variables
$\bm{x}\in\mathbb{R}^{256}$ with $50$-dimensional variables
$\tilde{\bm{x}}=\bm{N}^{T}\bm{x}$, where $\bm{N}\in\mathbb{R}^{256\times 50}$
is the basis matrix of natural cubic splines. We use $50$ basis functions with
knots uniformly placed over $1,\ldots,256$. To this end, we are dealing with
data of $n=1050$ and $p=50$.

Table 4: Performance of various methods in outlier detection for Phoneme data. Method | $\texttt{FDB}_{\rm L_{2}}$ | $\texttt{FDB}_{\rm pro}$ | DetMCD | Sample
---|---|---|---|---
Number | $49$ | $49$ | $49$ | $13$
AUC | $0.9895$ | $0.9895$ | $0.9895$ | $0.6115$
Time | $0.5390$ | $0.5225$ | $1.7911$ | $0.2865$

We take “ao” and “iy” to be regular cases and outliers, respectively, and
perform outlier detection. To be more specific, we calculate the robust
Mahalanobis distances based on FDBs and DetMCD, and classical Mahalanobis
distances based on Sample. Then, we treat cases with the top $50$ largest
distances as outliers, and the remaining ones as regular cases. Table 4
records the number of “iy” that are flagged as outliers, the area under the
ROC curve (AUC), and the average computational time (second) of $100$
replicates for these methods. Again, the proposed methods and DetMCD perform
similarly and they all outperform Sample with higher AUCs. In addition, the
proposed methods reveal advantages in computation time.

###  6.3 Denoise for MNIST data

The Modified National Institute of Standards and Technology (MNIST) database
is widely used for training various image processing systems. It contains a
large set of handwritten images representing the digits zero through nine.
Each digit is stored as a gray-scale image with a size $p=28\times 28$.
Training data $\bm{X}$ and testing data $\bm{X}_{t}$ of $n=10000$ are randomly
selected from MNIST. Noises
$\bm{\epsilon}\sim\mathcal{N}(\mathbf{0}_{p},60^{2}\mathbf{I}_{p})$ are
generated and added to $20\%$ of the images in the training set and all images
in the testing set. Our task is to denoise the testing images, which has been
considered in Schreurs et al. (2021).

The detailed process is similar to the robust PCA in Section 6.1. First, we
apply $\texttt{FDB}_{\rm L_{2}}$, $\texttt{FDB}_{\rm pro}$, DetMCD, and Sample
to training images and obtain the estimated location and scatter.
$\texttt{FDB}_{\rm pro}$ with $k=1000$, denoted as $\texttt{FDB}_{\rm
pro1000}$, is also used for better comparison. Here all methods have
comparable computational time. Then, we calculate the eigenvectors of the
robust estimated scatter. Next, we project the testing data to the subspace
spanned by the first $K$ eigenvectors and then transform the projected data,
i.e., scores, back to the original space. We denote the reconstructed data as
$\hat{\bm{X}}_{t}$. Fig. 8 illustrates $\hat{\bm{X}}_{t}$ obtained by various
methods with $K=75$. DetMCD is removed since it returns an error of high
condition numbers and hence is not applicable to this example. We can see that
the denoised images for the proposed methods are more clear than those for the
Sample, which verifies the influence of adding noise on evaluating scatter and
the efficiency of the proposed methods. As in Schreurs et al. (2021), we also
calculate the mean absolute error
$\mathrm{MAE}=\sum_{i=1}^{n}\sum_{j=1}^{p}|x_{t(ij)}-\hat{x}_{t(ij)}|/(np)$
between the original and denoised images. Table 5 shows that the MAE for
Sample is obviously larger than those for proposed methods.

Figure 8: The first and second columns show images before and after adding noise, respectively. The third to the sixth columns show denoised images obtained by $\texttt{FDB}_{\rm L_{2}}$, $\texttt{FDB}_{\rm pro}$, $\texttt{FDB}_{\rm pro1000}$, and Sample with $K=75$, respectively. Table 5: MAEs for proposed robust methods and classical method with various $K$s. $K$ | $15$ | $30$ | $45$ | $60$ | $75$
---|---|---|---|---|---
$\texttt{FDB}_{\rm L_{2}}$ | $23.04525$ | $19.27832$ | $17.52109$ | $16.63028$ | $16.26773$
$\texttt{FDB}_{\rm pro}$ | $23.06536$ | $19.27085$ | $17.51270$ | $16.62653$ | $16.26147$
$\texttt{FDB}_{\rm pro1000}$ | $23.09654$ | $19.30584$ | $17.54031$ | $16.65916$ | $16.29452$
Sample | $24.58077$ | $20.60175$ | $18.92803$ | $18.16014$ | $17.95447$

###  6.4 Outlier detection for Musk data

Musk data is commonly used in high-dimensional classification and outlier
detection problems (Aggarwal and Sathe, 2015; Porwal and Mukund, 2017). The
original data, which can be found in UCI includes $6598$ samples, divide into
musk and non-musk classes. Each sample has $p=166$ features characterizing the
molecule structure. Here we use a preprocessed musk dataset of $n=3062$, which
consists of $2965$ non-musk samples as inliers and $97$ musk sample as
outliers. Our task is to detect the outliers with the proposed methods, DetMCD
and Sample. As shown in Fig. 9, our proposed methods and DetMCD can exactly
pick out outliers, which are represented as red points, whereas Sample would
take many regular samples as outliers. Furthermore, the computational time of
$\texttt{FDB}_{\rm L_{2}}$, $\texttt{FDB}_{\rm pro}$, and DetMCD are $13.547$,
$12.884$, and $53.391$ seconds, respectively. The proposed methods are
computationally more efficient than DetMCD.

(a) $\texttt{FDB}_{\rm L_{2}}$

(b) $\texttt{FDB}_{\rm pro}$

(c) DetMCD

(d) Sample

Figure 9: (a), (b), and (c) are robust distances obtained by
$\texttt{FDB}_{\rm L_{2}}$, $\texttt{FDB}_{\rm pro}$, and DetMCD,
respectively. (d) Mahalanobis distance obtained by Sample.

##  7 Discussion

MCD-type methods suffer from high computational complexity due to the
iteratively used C-step. To tackle the issue, we directly approximate the MCD
subset with the trimmed subset induced by statistical depth. Two depth
notions, the projection depth and the $L_{2}$ depth, are recommended due to
high computational efficiency and robustness. In addition, we establish the
equivalence between the desired MCD subset and the trimmed subset induced by
the projection depth. Bypassing the iteration of the C-step, we manage to
reduce the computational complexity from $O(\psi np^{2}+\psi p^{3})$ to
$O(knp)$ and $O(n^{2}p)$ with $\texttt{FDB}_{\rm pro}$ and $\texttt{FDB}_{\rm
L_{2}}$, respectively. Moreover, the proposed estimators also reach the same
level of robustness as the MCD estimator. We conduct extensive simulation
studies and show that our estimators are comparable with MCD-type estimators
for low-dimensional data and significantly outperform MCD-type estimators for
high-dimensional cases.

The real data examples provide strong evidence that FDB is a valuable
complement to the toolset of robust multivariate analysis, including but not
limited to PCA, LDA, image denoise and outlier detection. The FDB algorithm
may benefit other applications that directly or indirectly rely on robust
covariance matrix estimation, such as robust linear regression (Coakley and
Hettmansperger, 1993), regression with continuous and categorical regressors
(Hubert and Rousseeuw, 1997), MCD-regression (Rousseeuw et al., 2004b),
multivariate least trimmed squares estimation (Agulló et al., 2008), and
robust errors-in-variables regression (Fekri and Ruiz-Gazen, 2004).

The present study primarily focuses on data from elliptical symmetric
distributions with enough samples, which may be violated in practical
scenarios (Schreurs et al., 2021). In Section S3.1 of the Supplementary
Material, we evaluate the effectiveness of our estimators when dealing with
skewed distributions. In addition, we did preliminary work to extend the
proposed algorithm to the scenario of “small $n$, large $p$” and evaluate our
idea with a real-world dataset in Section S3.2 of the Supplementary Material.
Our preliminary exploration shows promising results for both cases. Further
investigation is needed to address more general scenarios. One possible
solution is to adapt depth notions applicable to more general distributions to
obtain the $h$-subset and then apply the kernel trick to map the subset to a
feature space, where outlier detection can be conducted. The computational
time is significantly reduced for high-dimensional scenarios ($n>p$); however,
the estimation accuracy still desires further improvement. Rather than a
shrinkage estimator, it is also of interest to extend the MCD framework by
considering a low-rank and sparse estimator to alleviate the curse of
dimensionality.

## Supplementary materials

  * •

We provide an R-package named FDB and R codes of the FDB algorithm proposed in
this paper.

  * •

The file of supplement involves proofs of theoretical results, additional
simulation results as well as preliminary explorations of potential
extensions.

## Acknowledgement

We are very grateful to three anonymous referees, an associate editor, and the
Editor for their valuable comments that have greatly improved the manuscript.
The first two authors contribute equally to the paper.

## References

  * Aggarwal and Sathe (2015) Aggarwal, C. C. and Sathe, S. (2015), “Theoretical foundations and algorithms for outlier ensembles,” ACM SIGKDD Explorations Newsletter, 17, 24–47.
  * Agulló et al. (2008) Agulló, J., Croux, C., and Van Aelst, S. (2008), “The multivariate least-trimmed squares estimator,” Journal of Multivariate Analysis, 99, 311–338.
  * Alqallaf et al. (2009) Alqallaf, F., Van Aelst, S., Yohai, V. J., and Zamar, R. H. (2009), “Propagation of outliers in multivariate data,” The Annals of Statistics, 311–331.
  * Boudt et al. (2020) Boudt, K., Rousseeuw, P. J., Vanduffel, S., and Verdonck, T. (2020), “The minimum regularized covariance determinant estimator,” Statistics and Computing, 30, 113–128.
  * Butler et al. (1993) Butler, R., Davies, P., and Jhun, M. (1993), “Asymptotics for the minimum covariance determinant estimator,” The Annals of Statistics, 21, 1385–1400.
  * Cator and Lopuhaä (2012) Cator, E. A. and Lopuhaä, H. P. (2012), “Central limit theorem and influence function for the MCD estimators at general multivariate distributions,” Bernoulli, 18, 520–551.
  * Coakley and Hettmansperger (1993) Coakley, C. W. and Hettmansperger, T. P. (1993), “A bounded influence, high breakdown, efficient regression estimator,” Journal of the American Statistical Association, 88, 872–880.
  * Croux and Haesbroeck (2000) Croux, C. and Haesbroeck, G. (2000), “Principal component analysis based on robust estimators of the covariance or correlation matrix: influence functions and efficiencies,” Biometrika, 87, 603–618.
  * De Ketelaere et al. (2020) De Ketelaere, B., Hubert, M., Raymaekers, J., Rousseeuw, P. J., and Vranckx, I. (2020), “Real-time outlier detection for large datasets by RT-DetMCD,” Chemometrics and Intelligent Laboratory Systems, 199, 103957\.
  * Debruyne and Hubert (2009) Debruyne, M. and Hubert, M. (2009), “The influence function of the Stahel–Donoho covariance estimator of smallest outlyingness,” Statistics & probability letters, 79, 275–282.
  * Donoho (1982) Donoho, D. L. (1982), “Breakdown Properties of Multivariate Location Estimators,” Ph. D. Qualifying Paper, Harvard University.
  * Dyckerhoff et al. (2021) Dyckerhoff, R., Mozharovskyi, P., and Nagy, S. (2021), “Approximate computation of projection depths,” Computational Statistics & Data Analysis, 157, 107166.
  * Fekri and Ruiz-Gazen (2004) Fekri, M. and Ruiz-Gazen, A. (2004), “Robust weighted orthogonal regression in the errors-in-variables model,” Journal of Multivariate Analysis, 88, 89–108.
  * Hardin and Rocke (2004) Hardin, J. and Rocke, D. M. (2004), “Outlier detection in the multiple cluster setting using the minimum covariance determinant estimator,” Computational Statistics & Data Analysis, 44, 625–638.
  * Hastie et al. (2009) Hastie, T., Tibshirani, R., and Friedman, J. (2009), The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Springer Science & Business Media.
  * Hubert et al. (2018) Hubert, M., Debruyne, M., and Rousseeuw, P. J. (2018), “Minimum covariance determinant and extensions,” Wiley Interdisciplinary Reviews: Computational Statistics, 10, e1421.
  * Hubert and Rousseeuw (1997) Hubert, M. and Rousseeuw, P. J. (1997), “Robust regression with both continuous and binary regressors,” Journal of Statistical Planning and Inference, 57, 153–163.
  * Hubert et al. (2005a) Hubert, M., Rousseeuw, P. J., and Branden, K. V. (2005a), “ROBPCA: A New Approach to Robust Principal Component Analysis,” Technometrics, 47, 64 – 79.
  * Hubert et al. (2008) Hubert, M., Rousseeuw, P. J., and Van Aelst, S. (2008), “High-breakdown robust multivariate methods,” Statistical Science, 23, 92–119.
  * Hubert et al. (2005b) Hubert, M., Rousseeuw, P. J., and Vanden Branden, K. (2005b), “ROBPCA: a new approach to robust principal component analysis,” Technometrics, 47, 64–79.
  * Hubert et al. (2012) Hubert, M., Rousseeuw, P. J., and Verdonck, T. (2012), “A deterministic algorithm for robust location and scatter,” Journal of Computational and Graphical Statistics, 21, 618–637.
  * Hubert and Van Driessen (2004) Hubert, M. and Van Driessen, K. (2004), “Fast and robust discriminant analysis,” Computational Statistics & Data Analysis, 45, 301–320.
  * Kalina and Tichavskỳ (2021) Kalina, J. and Tichavskỳ, J. (2021), “The minimum weighted covariance determinant estimator for high-dimensional data,” Advances in Data Analysis and Classification, 1–23.
  * Liu (1990) Liu, R. Y. (1990), “On a notion of data depth based on random simplices,” The Annals of Statistics, 18, 405–414.
  * Lopuhaa and Rousseeuw (1991) Lopuhaa, H. P. and Rousseeuw, P. J. (1991), “Breakdown points of affine equivariant estimators of multivariate location and covariance matrices,” The Annals of Statistics, 229–248.
  * Mahalanobis (1936) Mahalanobis, P. C. (1936), “On the generalized distance in statistics,” National Institute of Science of India, 12, 49–55.
  * Maronna and Zamar (2002) Maronna, R. A. and Zamar, R. H. (2002), “Robust estimates of location and dispersion for high-dimensional datasets,” Technometrics, 44, 307–317.
  * Milo (1990) Milo, W. (1990), “Multivariate statistics: a practical approach,” Journal of the Royal Statistical Society: Series A, 153, 112–113.
  * Mosler and Mozharovskyi (2022) Mosler, K. and Mozharovskyi, P. (2022), “Choosing among notions of multivariate depth statistics,” Statistical Science, 37, 348–368.
  * Niinimaa and Oja (1995) Niinimaa, A. and Oja, H. (1995), “On the influence functions of certain bivariate medians,” Journal of the Royal Statistical Society: Series B (Methodological), 57, 565–574.
  * Oja (1983) Oja, H. (1983), “Descriptive statistics for multivariate distributions,” Statistics & Probability Letters, 1, 327–332.
  * Pison et al. (2003) Pison, G., Rousseeuw, P. J., Filzmoser, P., and Croux, C. (2003), “Robust factor analysis,” Journal of Multivariate Analysis, 84, 145–172.
  * Porwal and Mukund (2017) Porwal, U. and Mukund, S. (2017), “Outlier detection by consistent data selection method,” arXiv preprint arXiv:1712.04129.
  * Roelant et al. (2009) Roelant, E., Van Aelst, S., and Willems, G. (2009), “The minimum weighted covariance determinant estimator,” Metrika, 70, 177–204.
  * Rousseeuw et al. (2004a) Rousseeuw, P., Van Aelst, S., Driessen, K., and Agullo, J. (2004a), “Robust multivariate regression,” Technometrics, 46, 293–305.
  * Rousseeuw (1984) Rousseeuw, P. J. (1984), “Least median of squares regression,” Journal of the American Statistical Association, 79, 871–880.
  * Rousseeuw and Croux (1993) Rousseeuw, P. J. and Croux, C. (1993), “Alternatives to the median absolute deviation,” Journal of the American Statistical Association, 88, 1273–1283.
  * Rousseeuw and Driessen (1999) Rousseeuw, P. J. and Driessen, K. V. (1999), “A fast algorithm for the minimum covariance determinant estimator,” Technometrics, 41, 212–223.
  * Rousseeuw et al. (2004b) Rousseeuw, P. J., Van Aelst, S., Van Driessen, K., and Gulló, J. A. (2004b), “Robust multivariate regression,” Technometrics, 46, 293–305.
  * Rousseeuw and Van Zomeren (1990) Rousseeuw, P. J. and Van Zomeren, B. C. (1990), “Unmasking multivariate outliers and leverage points,” Journal of the American Statistical Association, 85, 633–639.
  * Schreurs et al. (2021) Schreurs, J., Vranckx, I., Ketelaere, B. D., Hubert, M., Suykens, J. A. K., and Rousseeuw, P. J. (2021), “Outlier detection in non-elliptical data by kernel MRCD,” Statistics and Computing, 31, 1–18.
  * Tukey (1975) Tukey, J. W. (1975), “Mathematics and the picturing of data,” in Proceedings of the International Congress of Mathematicians, Vancouver, 1975, vol. 2, pp. 523–531.
  * Vardi and Zhang (2000) Vardi, Y. and Zhang, C.-H. (2000), “The multivariate L1-median and associated data depth,” Proceedings of the National Academy of Sciences, 97, 1423–1426.
  * Willems et al. (2009) Willems, G., Joe, H., and Zamar, R. (2009), “Diagnosing multivariate outliers detected by robust estimators,” Journal of Computational and Graphical Statistics, 18, 73–91.
  * Yohai and Zamar (1988) Yohai, V. J. and Zamar, R. H. (1988), “High breakdown-point estimates of regression by means of the minimization of an efficient scale,” Journal of the American Statistical Association, 83, 406–413.
  * Zuo (2006) Zuo, Y. (2006), “Multidimensional trimming based on projection depth,” The Annals of Statistics, 34, 2211–2251.
  * Zuo and Serfling (2000a) Zuo, Y. and Serfling, R. (2000a), “General notions of statistical depth function,” The Annals of Statistics, 28, 461–482.
  * Zuo and Serfling (2000b) — (2000b), “Structural properties and convergence results for contours of sample statistical depth functions,” The Annals of Statistics, 28, 483–499.