# Statistically Optimal Robust Mean and Covariance Estimation for Anisotropic
Gaussians

Arshak Minasyan and Nikita Zhivotovskiy CREST, ENSAE, Institut Polytechnique
de Paris<EMAIL_ADDRESS>of California, Berkeley.
Department of Statistics<EMAIL_ADDRESS>
###### Abstract

Assume that $X_{1},\ldots,X_{N}$ is an $\varepsilon$-contaminated sample of
$N$ independent Gaussian vectors in $\mathbb{R}^{d}$ with mean $\mu$ and
covariance $\Sigma$. In the strong $\varepsilon$-contamination model we assume
that the adversary replaced an $\varepsilon$ fraction of vectors in the
original Gaussian sample by any other vectors. We show that there is an
estimator $\widehat{\mu}$ of the mean satisfying, with probability at least
$1-\delta$, a bound of the form

$\|\widehat{\mu}-\mu\|_{2}\leqslant
c\left(\sqrt{\frac{\operatorname{Tr}(\Sigma)}{N}}+\sqrt{\frac{\|\Sigma\|\log(1/\delta)}{N}}+\varepsilon\sqrt{\|\Sigma\|}\right),$

where $c>0$ is an absolute constant and $\|\Sigma\|$ denotes the operator norm
of $\Sigma$. In the same contaminated Gaussian setup, we construct an
estimator $\widehat{\Sigma}$ of the covariance matrix $\Sigma$ that satisfies,
with probability at least $1-\delta$,

$\left\|\widehat{\Sigma}-\Sigma\right\|\leqslant
c\left(\sqrt{\frac{\|\Sigma\|\operatorname{Tr}(\Sigma)}{N}}+\|\Sigma\|\sqrt{\frac{\log(1/\delta)}{N}}+\varepsilon\|\Sigma\|\right).$

Both results are optimal up to multiplicative constant factors. Despite the
recent significant interest in robust statistics, achieving both dimension-
free bounds in the canonical Gaussian case remained open. In fact, several
previously known results were either dimension-dependent and required $\Sigma$
to be close to identity, or had a sub-optimal dependence on the contamination
level $\varepsilon$.

As a part of the analysis, we derive sharp concentration inequalities for
central order statistics of Gaussian, folded normal, and chi-squared
distributions.

##  1 Robust multivariate mean estimation

Arguably the first rigorously studied question in robust statistics is the
mean (or the location parameter) estimation for contaminated Gaussian
distributions [Hub64]. A natural extension of this question is a problem of
multivariate Gaussian mean estimation when the data is contaminated by a
malicious adversary. When working with uncontaminated data, the celebrated
Borell, Tsirelson-Ibragimov-Sudakov [Bor75, TIS76] Gaussian concentration
inequality implies the sharp non-asymptotic bound for the performance of the
sample mean: If $X_{1}\ldots,X_{N}$ are independent Gaussian random vectors in
$\mathbb{R}^{d}$ with mean $\mu$ and covariance $\Sigma$, then, with
probability at least $1-\delta$,

$\Bigg{\|}\frac{1}{N}\sum\limits_{i=1}^{N}X_{i}-\mu\Bigg{\|}_{2}\leqslant\sqrt{\frac{\operatorname{Tr}(\Sigma)}{N}}+\sqrt{\frac{2\|\Sigma\|\log(1/\delta)}{N}}.$
(1)

Our question is to estimate the mean $\mu$ of a Gaussian random vector when an
$\varepsilon$-fraction of all observations is corrupted by a malicious
adversary, who knows both the clean sample and our estimator. We focus on the
multivariate case and make no assumptions on corrupted observations. This
model is usually called the strong contamination model [DK19, DKP20] (for the
exact definition of our model see [DM22, Definition 1]) and includes many
known contamination models such as Huber’s $\varepsilon$-contamination model
[Hub64]. Such a contaminated sample of Gaussian vectors will be referred to as
the $\varepsilon$-_contaminated sample_. Clearly, the sample mean can be
compromised even if there is a single outlier, so we will be aiming to provide
an analog of (1) for a different estimator. Despite significant recent
progress in robust statistics (we refer to the recent surveys where both the
statistical [LM19a] and algorithmic [DK19] aspects are discussed in detail),
there is still no sharp analog of the Gaussian concentration inequality (1)
when the $\varepsilon$-strong adversarial contamination is allowed. Although
the Gaussian case is historically the starting point in the theory of robust
statistics, our question remains open even from an information-theoretic point
of view, without considering computational aspects. Before stating our first
bound, we need an additional definition. Given a covariance matrix $\Sigma$,
its effective rank is defined as

$\mathbf{r}(\Sigma)=\frac{\operatorname{Tr}(\Sigma)}{\|\Sigma\|},$

where $\operatorname{Tr}(\Sigma)$ is the trace of matrix $\Sigma$. Obviously,
$1\leqslant\mathbf{r}(\Sigma)\leqslant d$ for a $d$ by $d$ covariance matrix
$\Sigma$, but can be much smaller if the distribution of the data is
anisotropic and is defined by several principal directions. We are now ready
to present our first bound.

######  Theorem 1 (Robust mean estimation in the Gaussian case).

There are absolute constants $c_{1},c_{2}>0$ such that the following holds.
Assume that $X_{1},\ldots,X_{N}$ is an $\varepsilon$-contaminated sample of
Gaussian random vectors in $\mathbb{R}^{d}$ with mean $\mu$ and covariance
$\Sigma$. Let $\varepsilon<c_{1}$, then there is an estimator $\widehat{\mu}$
satisfying, with probability at least $1-\delta$,

$\|\widehat{\mu}-\mu\|_{2}\leqslant
c_{2}\sqrt{\|\Sigma\|}\left(\sqrt{\frac{\mathbf{r}(\Sigma)}{N}}+\sqrt{\frac{\log(1/\delta)}{N}}+\varepsilon\right).$

Up to multiplicative constant factors, no estimator of the Gaussian mean can
perform better.

Although the generality of our result is its main strength, we focus for a
moment on the _isotropic case_ , that is $\Sigma=I_{d}$. In this case, if the
adversary corrupts at most $O(\sqrt{dN})$ elements of the sample, we still get
the optimal performance (1) of the sample mean in the setup where the data is
not contaminated. This dependence stands in contrast with the more usual
$\sqrt{\varepsilon}$-dependence on the contamination level that is known to be
achievable under the two moments assumption [LM21, DKP20]. In particular, the
latter results only allow $O(d)$ outliers to maintain the optimal performance.

Another interesting aspect of our analysis is that we do not make any
assumptions on the sample size. In comparison, the existing estimators that
recover the optimal bound in the isotropic case [CGR18, DK19] require
$N\geqslant c(d+\log(1/\delta))$, or $N\geqslant cd\varepsilon^{-2}$ as in
[DL22], where $c>0$ is some absolute constant. We will encounter some weaker
assumptions in Section 6, but only when tuning a single real-valued parameter
of our estimator.

In the context of mean estimation of anisotropic sub-Gaussian distributions,
the sharpest known bound is due to Lugosi and Mendelson [LM21]. These authors
proposed a multivariate version of a trimmed mean estimator that achieves the
following error rate

$\|\widehat{\mu}_{\textrm{LM}}-\mu\|_{2}\leqslant
c\sqrt{\|\Sigma\|}\left(\sqrt{\frac{\mathbf{r}(\Sigma)}{N}}+\sqrt{\frac{\log(1/\delta)}{N}}+\varepsilon\sqrt{\log\left(\frac{1}{\varepsilon}\right)}\right).$
(2)

The same rate has also been provided by Dalalyan and the first author of this
paper [DM22] for a different estimator under an additional assumption that
$\Sigma$ is known, though using a computationally efficient algorithm. One may
think that the presence of an additional $\sqrt{\log(1/\varepsilon)}$ term is
an artifact of the analysis in [LM21] and [DM22]. This is in fact not true,
and the presence of this term is known to be necessary for sub-Gaussian
distributions [CDG19, LM21]. By trimming the observations as in [LM21] one can
lose some of the specific properties of the Gaussian distribution. Thus, our
result asks for both a different estimator and a different analysis. The
presence of an additional $\sqrt{\log(1/\varepsilon)}$ term is also
interesting from the computational perspective. In particular, the authors of
[DKS17] argue that for any polynomial time Statistical Query algorithm the
factor $\sqrt{\log(1/\varepsilon)}$ is unavoidable in the error bound (see
also [HL19]). We note that in the rich literature on robust statistics, a
number of polynomial running time estimators were proposed with
$\varepsilon\sqrt{\log(1/\varepsilon)}$ dependence on the contamination level
[LRV16, DKK+17, DK19, DM22, BMD22].

It is worth mentioning that there are several results showing that the linear
dependence on $\varepsilon$ can be achieved in the isotropic case, where
$\Sigma$ is identity (or close to it). We refer to the analysis of the Tukey
median, the direction-dependent median as well as the Stahel-Donoho median of
means in respectively [CGR18, DK19, DL22]. Apart from the fact that the
isotropic assumption is quite restrictive when working with real data, the
presence of VC-type or sphere covering arguments is inherent to the analysis
of existing dimension-dependent estimators. Unfortunately, these arguments
cannot help us to prove the dimension-free bound of Theorem 1. It is well
understood that the bound on the Gaussian complexity of ellipsoids, which
corresponds to the $\sqrt{\operatorname{Tr}(\Sigma)}$ term in our result,
follows neither from the ellipsoid covering, nor from the Dudley integral, nor
from VC-type arguments. This feature is especially pronounced in the Gaussian
covariance estimation problem discussed in Section 2. We refer to Chapter 2.5
in the monograph of Talagrand [Tal14] for a thorough analysis of these
questions.

The starting point of our analysis is the folklore property of the sample
median of the Gaussian distribution in the one-dimensional case. We denote the
sample median by $\operatorname{Med}(\cdot)$ in what follows. If
$X_{1},\ldots,X_{N}$ is an $\varepsilon$-contaminated sample of independent
standard Gaussians with mean $\mu$ and variance $\sigma^{2}$, then, with
probability at least $1-\delta$,

$\left|\operatorname{Med}(X_{1},\ldots,X_{N})-\mu\right|\leqslant
c\sigma\left(\sqrt{\frac{\log(1/\delta)}{N}}+\varepsilon\right),$

whenever $N\geqslant c\log(1/\delta)$ and $\varepsilon$ is smaller than some
absolute constant. When going to higher dimensions, instead of working with
Tukey’s median, whose sharp analysis is only known in the isotropic case
[CGR18], or Stahel-Donoho-type estimators as in [DL22], we base our solution
on what we call the _smoothed median estimator_.

Let $x_{1},\ldots,x_{N}$ be a set of vectors in $\mathbb{R}^{d}$, and let
$\xi=(\xi_{1},\ldots,\xi_{N})$ be a zero mean random vector in
$\mathbb{R}^{N}$ whose covariance matrix $H$ is given by
$H_{i,j}=\beta^{-1}\langle x_{i},x_{j}\rangle$ for $i,j=1,\ldots,N$. That is,
$H$ is proportional to the Gram matrix of the original data. Here $\beta>0$ is
any integer (chosen by the statistician) satisfying
$\mathbf{r}(\Sigma)/10\leqslant\beta\leqslant 10\mathbf{r}(\Sigma)$. For any
direction $v\in S^{d-1}$, we are interested in the following quantity we call
the _smoothed median_ :

$\operatorname{SmoothMed}_{v}(x_{1},\ldots,x_{N})={\mathbf{E}}_{\xi}\operatorname{Med}(\langle
x_{1},v\rangle+\xi_{1},\ldots,\langle x_{N},v\rangle+\xi_{N}).$ (3)

Observe that $\operatorname{SmoothMed}_{v}(x_{1},\ldots,x_{N})$ is a function
of $x_{1},\ldots,x_{N}$ and $v$. The estimator of Theorem 1 has a simple form.
Given an $\varepsilon$-contaminated sample $X_{1},\ldots,X_{N}$, we set

$\widehat{\mu}=\mathop{\mathrm{arg}\,\mathrm{min}}_{\nu\in\mathbb{R}^{d}}\sup\limits_{v\in
S^{d-1}}\left|\operatorname{SmoothMed}_{v}(X_{1},\ldots,X_{N})-\langle\nu,v\rangle\right|.$
(4)

######  Remark 1.

In Section 6, we discuss how one may choose an integer $\beta$ satisfying
$\mathbf{r}(\Sigma)/10\leqslant\beta\leqslant 10\mathbf{r}(\Sigma)$ based only
on the contaminated sample. Importantly, we avoid a sample-splitting approach
when tuning this parameter.

From the practical perspective, our estimator has complexity exponential in
dimension. This is a typical limitation for all existing estimators that have
a linear dependence on the contamination level $\varepsilon$ in the Gaussian
case111Recall that existing estimators of this kind lead to dimension-
dependent bounds.. However, in the case when the dimension $d$ is small
enough, one can replace the computation over the sphere $S^{d-1}$ by an
appropriate $\varepsilon$-net and approximate the smoothing integration
uniformly over all the elements of this net using a Monte Carlo sampling
technique.

The appearance of the smoothed median follows from the proof technique that
guarantees a dimension-free nature of our bound. Our approach uses the so-
called _PAC-Bayesian lemma_ , whose applications were pioneered by O. Catoni
and co-authors [AC11, Cat16, CG17] in the context of mean/covariance
estimation/linear regression in the heavy-tailed setup. Our application
further develops these techniques but in the context of adversarial
contamination. An additional discussion appears in Section 4.

#### Notation.

Throughout the text $c,c_{1},c_{2},\ldots$ denote absolute constants that may
change from line to line. For two positive semi-definite matrices $A$ and $B$
we write $A\preceq B$, if $B-A$ is positive semi-definite. The symbol
$\|\cdot\|$ denotes the operator norm of a matrix or the Euclidean norm of a
vector depending on the context. Let $\mathbb{S}_{+}^{d}$ denote the set of
$d\times d$ positive semi-definite matrices. The symbol $I_{d}$ denotes the
identity $d\times d$ matrix. We denote the indicator of the event $A$ by
$\bm{1}\left(A\right)$. For any integer $N$, $[N]$ is the shortened notation
of the set $\\{1,\ldots,N\\}$. For a random variable $Y$ and $\alpha\in[1,2]$,
its $\psi_{\alpha}$ Orlicz norm is defined as follows

$\|Y\|_{\psi_{\alpha}}=\inf\\{c>0:{\mathbf{E}}\exp(|Y|^{\alpha}/c^{\alpha})\leqslant
2\\}.$

Using the standard convention, we say that $\|\cdot\|_{\psi_{2}}$ is the sub-
Gaussian norm and $\|\cdot\|_{\psi_{1}}$ is the sub-exponential norm. Let
$\mathcal{KL}(\rho,\gamma)=\int\log\left(\frac{d\rho}{d\gamma}\right)d\rho$
denote the Kullback-Leibler divergence between the two measures $\rho$ and
$\gamma$ such that $\rho\ll\gamma$. The notation $\rho\ll\gamma$ means that
the measure $\rho$ is absolutely continuous with respect to the measure
$\gamma$. We define the order statistics. Given a set of real numbers
$x_{1},\ldots,x_{n}$, let $x_{(1)},\ldots,x_{(n)}$ denote their non-decreasing
rearrangement. That is, $x_{(1)}\leqslant x_{(2)}\leqslant\ldots\leqslant
x_{(n)}$. For $\alpha\in[0,1]$, assuming that $\alpha n$ is an integer, we set

$\operatorname{Quant}_{\alpha}(x_{1},\ldots,x_{n})=x_{(\alpha n)}.$

In particular, assuming for simplicity that $n$ is odd, the sample median is
given by

$\operatorname{Med}(x_{1},\ldots,x_{n})=x_{((n+1)/2)}.$

#### Related literature.

Robust statistics is a well-developed topic with several explanatory texts
published recently. In our context, the most relevant are the surveys of
Lugosi and Mendelson [LM19a], and of Diakonikolas and Kane [DK19]. Some
classical references on robust statistics include the textbooks [HRRS80,
Hub81, RL05]. When discussing covariance estimation, we refer to the survey
[KMR+19] where the focus is on heavy-tailed distributions. We also mention
several recent papers on covariance estimation [CGR18, MW22, AZ22, OR22],
where the focus is on adversarial contamination.

Instead of working with the Euclidean norm as in Theorem 1, some authors focus
on the Mahalanobis norm. That is, one aims to construct an estimator
$\widehat{\mu}$ such that
$(\widehat{\mu}-\mu)^{\mathsf{T}}\Sigma^{-1}(\widehat{\mu}-\mu)$ is small with
high probability. It appears that the bounds with respect to this norm are
necessarily dimension-dependent, and a simple VC-type/sphere covering argument
is sufficient to obtain the optimal rates of convergence [DL22]. Similar
observations are also valid for the covariance estimation problem. We focus on
the operator norm, where the analysis allows for dimension-free bounds. This
is not the case for the (weighted) Frobenius norm commonly analyzed in the
literature.

##  2 Covariance estimation

We now move to a more challenging problem of covariance estimation. For
simplicity, we assume that our uncontaminated distribution is zero mean. We
first need to present a sharp analog of inequality (1) in the case where no
contamination is allowed. Such a result has been shown only recently by
Koltchinskii and Lounici [KL17]. Their analysis is based on the generic
chaining for quadratic processes. This non-trivial approach is motivated by
the difficulty of replacing $d$ with the effective rank $\mathbf{r}(\Sigma)$.
Let us formulate their result. Assume that $Y_{1},\ldots,Y_{N}$ are
independent zero mean Gaussian vectors in $\mathbb{R}^{d}$ with covariance
$\Sigma$. There are absolute constants $c_{1},c_{2}>0$ such that, with
probability at least $1-\delta$,

$\left\|\frac{1}{N}\sum\limits_{i=1}^{N}Y_{i}Y_{i}^{\mathsf{T}}-\Sigma\right\|\leqslant
c_{1}\|\Sigma\|\left(\sqrt{\frac{\mathbf{r}(\Sigma)}{N}}+\sqrt{\frac{\log(1/\delta)}{N}}\right),$
(5)

provided that $N\geqslant c_{2}(\mathbf{r}(\Sigma)+\log(1/\delta))$. When
adversarial contamination is allowed, the sharpest known dimension-free result
is implied by the bound of Abdalla and the second author of this paper [AZ22].
Their work suggests a trimmed-mean-based estimator that achieves the rate

$\left\|\widehat{\Sigma}_{\textrm{AZh}}-\Sigma\right\|\leqslant
c_{1}\|\Sigma\|\left(\sqrt{\frac{\mathbf{r}(\Sigma)}{N}}+\sqrt{\frac{\log(1/\delta)}{N}}+\varepsilon\log\left(\frac{1}{\varepsilon}\right)\right),$
(6)

whenever $N\geqslant c_{2}(\mathbf{r}(\Sigma)+\log(1/\delta))$. The above
bound is valid for any sub-Gaussian distribution and cannot be improved in
general. However, similarly to the case of mean estimation, we expect a better
dependence on the contamination parameter $\varepsilon$ in the Gaussian case.
We also note that simpler versions of the bound (6) based on the median-of-
means estimators [MZ20] only lead to a $\sqrt{\varepsilon}$-dependence on the
contamination level.

######  Remark 2.

As a side note, when aiming for a sharp leading constant in (5) in the
uncontaminated setup, an almost optimal performance follows from the recent
result of Han [Han22] combined with the second order concentration inequality
derived in [KZ20]. A similar bound in the sub-Gaussian case with explicit
constants can be found in [Zhi21].

In our analysis, we first make some additional assumptions. Our estimator
depends on some additional parameters that could be pre-estimated based only
on the observed $\varepsilon$-contaminated sample. A careful analysis of pre-
estimation procedures is deferred to Section 6. For the rest of this section,
we assume that we have an access to the following quantities.

  1. 1.

There is an integer $\beta$ and a real number $\omega$ satisfying respectively

$\frac{1}{10}\mathbf{r}(\Sigma)\leqslant\beta\leqslant
10\mathbf{r}(\Sigma),\quad\text{and}\quad\frac{1}{10}\|\Sigma\|\leqslant\omega\leqslant
10\|\Sigma\|.$

  2. 2.

Let $H$ be a known positive semi-definite matrix. Assume that we know a real
number $\alpha=\alpha(H)$ satisfying

$|\alpha-\operatorname{Tr}(\Sigma H)|\leqslant c\operatorname{Tr}(\Sigma
H)\left(\sqrt{\frac{\mathbf{r}(\Sigma)+\log(1/\delta)}{N}}+\varepsilon\right),$
(7)

for some absolute constant $c>0$.

  3. 3.

We have an access to a positive semi-definite matrix $G$ satisfying

$\frac{1}{10}\Sigma\preceq
G,\quad\textrm{and}\quad\operatorname{Tr}(G)\leqslant
10\operatorname{Tr}(\Sigma).$ (8)

Except for the matrix $G$, we only need to tune real valued parameters. This
can be usually done under the minimal assumption $N\geqslant
c(\mathbf{r}(\Sigma)+\log(1/\delta))$ for some absolute constant $c>0$.
Observe that the assumption $\Sigma\preceq 10G$ does not imply that $\Sigma$
is close to $G$ in the operator norm. At the same time, this assumption
requires some control over the smallest singular value of $\Sigma$. We show,
in particular, that whenever $N\geqslant c(d+\log(1/\delta))$, we can always
efficiently construct such a matrix $G$ based on contaminated data, while
still maintaining the dimension-free nature of our upper bound. Moreover, it
appears that our guarantees are uniform with respect to the choice of the
matrix $G$. One can rerun our estimator on the same data multiple times with
any admissible $G$ satisfying (8) without affecting the performance of our
estimator. We discuss this formally in Section 6.

######  Theorem 2 (Robust covariance estimation in the Gaussian case).

There are absolute constants $c_{1},c_{2}>0$ such that the following holds.
Assume that $X_{1},\ldots,X_{N}$ is an $\varepsilon$-contaminated sample of
zero mean Gaussian vectors in $\mathbb{R}^{d}$ with covariance $\Sigma$. Let
$\varepsilon<c_{1}$, then there is an estimator
$\widehat{\Sigma}=\widehat{\Sigma}_{\alpha,\beta,\omega,G,\varepsilon}(X_{1},\ldots,X_{N})$
satisfying, with probability at least $1-\delta$,

$\left\|\widehat{\Sigma}-\Sigma\right\|\leqslant
c_{2}\|\Sigma\|\left(\sqrt{\frac{\mathbf{r}(\Sigma)}{N}}+\sqrt{\frac{\log(1/\delta)}{N}}+\varepsilon\right).$

Up to multiplicative constant factors, no estimator of the covariance matrix
performs better.

We are now ready to define our estimator.

We first construct the following distribution. For any $v\in S^{d-1}$, let
$\rho_{v}$ be a distribution in $\mathbb{R}^{d}$ whose density $f_{v}$ is
given by

$f_{v}(x)=\frac{1}{p(2\pi\beta^{-1})^{d/2}}\exp\left(-\frac{\beta\|x-v\|^{2}}{2}\right)\bm{1}\left(\|G^{1/2}(x-v)\|\leqslant
100\sqrt{\omega}\right).$

Here $p>0$ is a normalization factor. Assume that $\theta$ is a random vector
distributed according to $\rho_{v}$. Let
$H={\mathbf{E}}_{\rho_{v}}(\theta-v)(\theta-v)^{\mathsf{T}}$ be a covariance
matrix of $\rho_{v}$ (by the symmetry of $\rho_{v}$ around $v$, the matrix $H$
does not depend on $v$). For some specifically chosen absolute constant $c>0$
and $\alpha=\alpha(H)$, define the set

$\mathcal{H}=\Bigg{\\{}\Gamma\in\mathbb{S}_{+}^{d}:|\operatorname{Tr}(\Gamma
H)-\alpha|\leqslant c\operatorname{Tr}(\Gamma
H)\Bigg{(}\sqrt{\frac{\mathbf{r}(G)+\log({1}/{\delta})}{N}}+\varepsilon\Bigg{)};\Gamma\preceq
10G;\|\Gamma\|\leqslant 10\omega\Bigg{\\}}.$ (9)

For an $\varepsilon$-contaminated sample $X_{1},\dots,X_{N}$, our estimator is
defined as follows:

$\widehat{\Sigma}=\mathop{\mathrm{arg}\,\mathrm{min}}_{\Gamma\in\mathcal{H}}\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\left|\operatorname{Med}\left(|\langle
X_{1},\theta\rangle|,\ldots,|\langle
X_{N},\theta\rangle|\right)-\Phi^{-1}(3/4)\sqrt{\theta^{\mathsf{T}}\Gamma\theta}\right|.$
(10)

This estimator is a more complex version of our smoothed median estimator.
First, instead of working with the Gaussian smoothing measure, we restrict
this distribution to an elliptic set
$\\{x\in\mathbb{R}^{d}:\|G^{1/2}(x-v)\|\leqslant 100\sqrt{\omega}\\}$. Second,
we need to restrict the eigenvalues of the output matrix and introduce the set
$\mathcal{H}$. Our estimator is related to minimizing the so-called _median
absolute deviation_ (see [DG92] for related definitions). The proof of Theorem
2 exploits the fact that quantiles of $|\langle X,\theta\rangle|$ are tightly
connected with corresponding variances. This is reflected in the term
$\Phi^{-1}(3/4)$ appearing in the definition of our estimator.

##  3 Concentration inequalities for sample quantiles

In this section, we obtain sub-Gaussian and sub-exponential concentration
inequalities for quantiles of i.i.d. observations sampled according to several
regular distributions. The analysis of sample quantiles is a standard question
in statistics. The early work of Kolmogorov [Kol31] focused on proving a
central limit theorem for the sample median of some symmetric distributions.
Subsequently, the focus was on explicit expansions for this limit law [Bur97].
Another line of research focused on studying the explicit formulas for the
distribution of order statistics. We refer to the monograph of David and
Nagaraja [DN04] on order statistics and to the monograph of De Haan and
Ferreira [DHF07] on the extreme value theory, a topic covering the properties
of smallest and largest values in samples. Finally, many authors focused on
the analysis of Bahadur’s representation of sample quantiles (see e.g.,
[Bah66, Kie67], [Sha03, Theorem 5.11]). Unfortunately, neither the exact
expressions for the distribution of sample quantiles nor various asymptotic
expansions lead to the exact concentration inequalities we are interested in.

Less is known about concentration inequalities for sample quantiles. Some
explicit non-asymptotic bounds appear in the monograph of Shao [Sha03, Section
5.3], where a reduction to a concentration of Bernoulli random variables is
made. Several related concentration inequalities appear in the work of
Boucheron and Thomas [BT12], though their bounds are not sharp enough for our
purposes. In particular, these authors provide a sub-exponential concentration
inequality for the sample median of the Gaussian distribution, while our
results will lead to sub-Gaussian concentration inequalities. More recent
results on sample quantiles are inspired by the problems in robust statistics.
In fact, the analysis in [CGR18, ABM19] provides some sharp bounds for sample
quantiles, though the existing bounds do imply the sub-Gaussian concentration
only for small enough deviations from the mean. Another line of results is due
to Bobkov and Ledoux [BL19]. Their results provide sharp concentration
inequalities for log-concave distributions (recall that the Gaussian
distribution is log-concave), but only lead to sub-exponential tails due to
their generality222For the special case of the uniform distribution on the
real line Bobkov and Ledoux [BL19] provide a sub-Gaussian concentration
inequality for all order statistics of the uniform distribution in $[0,1]$..

Our approach is quite simple, though, to the best of our knowledge, it is not
used explicitly in the literature. When proving concentration inequalities for
sample quantiles, we consider two regimes. For small deviations, we use the
regularity of the density function and follow the reduction to a concentration
of Bernoulli random variables as in [Sha03, CGR18, ABM19, DK19, Xia19], while
for large deviations we use the sub-Gaussian/sub-exponential tails of our
distribution. This leads to desired sharp concentration inequalities. We
discuss some straightforward extensions of our analysis in Section 7.

Before providing our first concentration inequality, recall that the
cumulative distribution function of a standard Gaussian is denoted by
$\Phi(\cdot)$. Denote its inverse by $\Phi^{-1}(\cdot)$.

######  Lemma 1 (Concentration for Gaussian quantiles).

There are absolute constants $c_{1},c_{2}>0$ such that the following holds.
Let $\varepsilon\in[0,1/4]$. Assume without loss of generality that
$(1/2\pm\varepsilon)N$ are integers. Let $Y_{1},\ldots,Y_{N}$ be a sample of
independent standard Gaussian random variables. Then, for any $t\geqslant 0$,

$\Pr(|Y_{((1/2\pm\varepsilon)N)}-\Phi^{-1}(1/2\pm\varepsilon)|\geqslant
t)\leqslant 2\exp(-c_{1}Nt^{2}).$

Equivalently,

$\left\|Y_{((1/2\pm\varepsilon)N)}-\Phi^{-1}(1/2\pm\varepsilon)\right\|_{\psi_{2}}\leqslant\frac{c_{2}}{\sqrt{N}}.$

###### Proof.

We only analyze the quantile $Y_{((1/2+\varepsilon)N)}$, since the analysis
for $Y_{((1/2-\varepsilon)N)}$ is the same. We analyze two parts of the tail
separately. First, we show that

$\Pr(Y_{((1/2+\varepsilon)N)}-\Phi^{-1}(1/2+\varepsilon)\geqslant
t)\leqslant\exp(-c_{1}Nt^{2}).$

This analysis is also split into two regimes. For some absolute constant
$C>0$, we first show the above inequality for $0\leqslant t\leqslant C$ and
then proceed with the case $t\geqslant C$. In the first regime, we follow the
standard reduction to binomial tails (see similar computations in [Sha03,
Theorem 5.9] and [Xia19]). Let $Z_{1},\ldots,Z_{N}$ be independent Bernoulli
random variables with expectation $p=1-\Phi(s)$ for some fixed
$s\in\mathbb{R}$. We have

$\Pr(Y_{((1/2+\varepsilon)N)}\geqslant
s)=\Pr\left(\sum\nolimits_{i=1}^{N}Z_{i}>(1/2-\varepsilon)N\right).$

We set $s=\Phi^{-1}(1/2+\varepsilon)+t$ and obtain, using that
${\mathbf{E}}Z_{i}=1-\Phi(\Phi^{-1}(1/2+\varepsilon)+t)$,

$\displaystyle\Pr(Y_{((1/2+\varepsilon)N)}\geqslant\Phi^{-1}(1/2+\varepsilon)+t)$
$\displaystyle=\Pr\left(\frac{1}{N}\sum\limits_{i=1}^{N}Z_{i}-{\mathbf{E}}Z_{i}>\Phi(\Phi^{-1}(1/2+\varepsilon)+t)-\frac{1}{2}-\varepsilon\right).$

Denoting $\varphi_{+}(t)=\Phi(\Phi^{-1}(1/2+\varepsilon)+t)-1/2-\varepsilon$,
we have by Hoeffding’s inequality applied to independent Bernoulli random
variables, whenever $\varphi_{+}(t)\geqslant 0$,

$\Pr(Y_{((1/2+\varepsilon)N)}-\Phi^{-1}(1/2+\varepsilon)\geqslant
t)\leqslant\exp\left(-2N(\varphi_{+}(t))^{2}\right).$ (11)

Let us lower bound the function $\varphi_{+}$. Using the density formula, we
have

$\displaystyle\varphi_{+}(t)$
$\displaystyle=\Phi(\Phi^{-1}(1/2+\varepsilon)+t)-\Phi(\Phi^{-1}(1/2+\varepsilon))$
$\displaystyle\geqslant\frac{t}{\sqrt{2\pi}}\exp\left(-(\Phi^{-1}(1/2+\varepsilon)+t)^{2}/2\right)$
$\displaystyle\geqslant\frac{t}{\sqrt{2\pi}}\exp\left(-(\Phi^{-1}(3/4)+t)^{2}/2\right).$

We combine these computations with the tail for large values of $t$. Since
$\varepsilon\leqslant 1/4$, we need that at least $N/4$ (assume that it is an
integer without loss of generality) of all observations are above
$\Phi^{-1}(1/2+\varepsilon)+t$. This can be controlled as follows

$\displaystyle\Pr(Y_{((1/2+\varepsilon)N)}-\Phi^{-1}(1/2+\varepsilon)\geqslant
t)$
$\displaystyle\leqslant\binom{N}{N/4}\left(\Pr(Y_{1}\geqslant\Phi^{-1}(1/2+\varepsilon)+t)\right)^{N/4}$
$\displaystyle\leqslant 2^{N}\left(\Pr(Y_{1}\geqslant t)\right)^{N/4}$
$\displaystyle\leqslant 2^{N}\exp(-Nt^{2}/8)$
$\displaystyle=\exp(N\log(2)-Nt^{2}/8)$
$\displaystyle\leqslant\exp(-Nt^{2}/16),$

whenever $t\geqslant 4\sqrt{\log(2)}$. The inequality (11) and the lower bound
on $\varphi_{+}(t)$ give us

$\Pr(Y_{((1/2+\varepsilon)N)}-\Phi^{-1}(1/2+\varepsilon)\geqslant
t)\leqslant\exp\left(\frac{-2Nt^{2}}{2\pi\exp((\Phi^{-1}(3/4)+4\sqrt{\log(2)})^{{2}})}\right),$

whenever $0\leqslant t\leqslant 4\sqrt{\log(2)}$. Combining two regimes and
adjusting the absolute constant, we prove an upper tail. Let us prove the
lower tail bound. The proof is similar, though the computations are slightly
different. We want to show for any $t\geqslant 0$,

$\Pr(Y_{((1/2+\varepsilon)N)}-\Phi^{-1}(1/2+\varepsilon)\leqslant-t)\leqslant\exp(-c_{1}Nt^{2}).$

We have for any $s\in\mathbb{R}$ and $Z_{1},\ldots,Z_{N}$ as above

$\Pr(Y_{((1/2+\varepsilon)N)}\leqslant
s)=\Pr\left(\sum\nolimits_{i=1}^{N}Z_{i}\leqslant(1/2-\varepsilon)N\right).$

Define

$\varphi_{-}(t)=\Phi(\Phi^{-1}(1/2+\varepsilon)-t)-1/2-\varepsilon.$

We have similarly

$\Pr(Y_{((1/2+\varepsilon)N)}-\Phi^{-1}(1/2+\varepsilon)\leqslant-t)\leqslant\exp\left(-2N(\varphi_{-}(t))^{2}\right).$

Now we lower bound the quantity $|\varphi_{-}(t)|$ as follows

$\displaystyle|\varphi_{-}(t)|$
$\displaystyle=\Phi(\Phi^{-1}(1/2+\varepsilon))-\Phi(\Phi^{-1}(1/2+\varepsilon)-t)$
$\displaystyle\geqslant\frac{t}{\sqrt{2\pi}}\exp\left(-\max\\{(\Phi^{-1}(1/2+\varepsilon))^{2},(\Phi^{-1}(1/2+\varepsilon)-t)^{2}\\}/2\right)$
$\displaystyle\geqslant\frac{t}{\sqrt{2\pi}}\exp\left(-\max\\{(\Phi^{-1}(3/4))^{2},t^{2}\\}/2\right).$

As above, we need to get the tail for large values of $t$. We need that at
least $N/{4}$ of all observations are below $\Phi^{-1}(1/2+\varepsilon)-t$. In
what follows, we assume $t\geqslant 2\Phi^{-1}(3/4)$. This can be controlled
as follows

$\displaystyle\Pr(Y_{((1/2+\varepsilon)N)}-\Phi^{-1}(1/2+\varepsilon)\leqslant-t)$
$\displaystyle=\Pr(-Y_{((1/2+\varepsilon)N)}\geqslant
t-\Phi^{-1}(1/2+\varepsilon))$
$\displaystyle\leqslant\binom{N}{N/{4}}\left(\Pr(-Y_{1}\geqslant
t-\Phi^{-1}(1/2+\varepsilon))\right)^{N/{4}}$ $\displaystyle\leqslant
2^{N}\left(\Pr(Y_{1}\geqslant t-\Phi^{-1}(3/4))\right)^{N/{4}}$
$\displaystyle\leqslant 2^{N}\left(\Pr(Y_{1}\geqslant t/2)\right)^{N/{4}}$
$\displaystyle\leqslant 2^{N}\exp(-Nt^{2}/{32})$
$\displaystyle\leqslant\exp(-Nt^{2}/{64}),$

whenever $t\geqslant\max\\{2\Phi^{-1}(3/4),{8}\sqrt{\log(2)}\\}$. The proof of
the lower tail follows. The union bound concludes the proof. Finally, our
bound on the $\psi_{2}$-norm follows from [Ver16, Proposition 2.5.2]. ∎

Our second result presents a similar concentration bound for the empirical
quantiles of i.i.d. observations drawn from $\chi^{2}_{1}$ distribution. This
distribution coincides with the distribution of the squared standard Gaussian
random variable. Denote the cumulative distribution function by
$F_{\chi^{2}_{1}}(\cdot)$ and its inverse by $F_{\chi^{2}_{1}}^{-1}(\cdot)$.
The key difference is that we only show a sub-exponential tail in this case.
We remark that when considering the $\chi^{2}_{k}$ distribution with
$k\geqslant 2$ degrees of freedom, the desired concentration inequality
follows from log-concavity and [BL19, Lemma 6.5].

######  Lemma 2 (Quantiles of the $\chi^{2}_{1}$ distribution).

There is an absolute constant $c_{1}>0$ such that the following holds. Assume
without loss of generality that $(1/2\pm\varepsilon)N$ are integers. Let
$Y_{1},\ldots,Y_{N}$ be a sample of independent $\chi^{2}_{1}$ random
variables and $\varepsilon\in[0,1/4]$. Then,

$\left\|Y_{((1/2\pm\varepsilon)N)}-F_{\chi^{2}_{1}}^{-1}(1/2\pm\varepsilon)\right\|_{\psi_{1}}\leqslant\frac{c_{1}}{\sqrt{N}}.$

###### Proof.

First, the density $g$ of $\chi_{1}^{2}$ is given by

$g(x)=\frac{\exp(-x/2)}{\sqrt{2\pi x}},\quad x>0.$ (12)

It is also easy to show that
$F^{-1}_{\chi_{1}^{2}}(1/2)=(\Phi^{-1}(3/4))^{2}$. Using the same notation
again, we denote
$\varphi_{+}(t)=F_{\chi_{1}^{2}}(F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)+t)-1/2-\varepsilon$.
Repeating the lines of the proof of Lemma 1, whenever $\varphi_{+}(t)\geqslant
0$, we have

$\Pr(Y_{((1/2+\varepsilon)N)}-F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)\geqslant
t)\leqslant\exp\left(-2N(\varphi_{+}(t))^{2}\right).$ (13)

Using (12), we have

$\displaystyle\varphi_{+}(t)$
$\displaystyle=F_{\chi_{1}^{2}}(F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)+t)-F_{\chi_{1}^{2}}(F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon))$
$\displaystyle\geqslant\frac{t\exp(-F_{\chi_{1}^{2}}(F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)+t))}{\sqrt{2\pi
F_{\chi_{1}^{2}}(F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)+t)}}.$

This gives us a sub-Gaussian tail for as long as $\varepsilon\leqslant 1/4$
and $t\geqslant 0$ is bounded by some absolute constant. By the concentration
of the $\chi^{2}$ distribution [LM00, Lemma 1] we have for $t\geqslant 0$,

$\Pr(Y_{1}\geqslant 1+\sqrt{2t}+t)\leqslant\exp(-t),\ \textrm{and thus},\
\Pr(Y_{1}\geqslant t)\leqslant\exp(-t/2),\ \textrm{whenever}\ t\geqslant
4+2\sqrt{3}.$

Since $\varepsilon\leqslant 1/4$, we need that at least $N/4$ of all
observations are above $F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)+t$. Therefore,
whenever $t\geqslant 4+2\sqrt{3}$, we have

$\displaystyle\Pr(Y_{((1/2+\varepsilon)N)}-F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)\geqslant
t)$
$\displaystyle\leqslant\binom{N}{N/4}\left(\Pr(Y_{1}\geqslant\Phi^{-1}(1/2+\varepsilon)+t)\right)^{N/4}$
$\displaystyle\leqslant\binom{N}{N/4}\left(\Pr(Y_{1}\geqslant t)\right)^{N/4}$
$\displaystyle\leqslant 2^{N}\exp(-Nt/8)$ $\displaystyle=\exp(N\log(2)-Nt/8)$
$\displaystyle\leqslant\exp(-Nt/16),$

where the last inequality requires additionally $t\geqslant 16\log(2)$.
Combining the above bounds and adjusting the absolute constant $c_{1}\geqslant
0$ we show that

$\Pr(Y_{((1/2+\varepsilon)N)}-F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)\geqslant
t)\leqslant\exp(-c_{1}N\min\\{t,t^{2}\\}).$

We continue with the bound on the lower tail. We want to show for any
$t\geqslant 0$,

$\Pr(Y_{((1/2+\varepsilon)N)}-F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)\leqslant-t)\leqslant\exp(-c_{2}Nt^{2}),$

where $c_{2}$ is an absolute constant. For $0\leqslant
t<F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)$ we define

$\varphi_{-}(t)=F^{-1}_{\chi_{1}^{2}}(F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)-t)-1/2-\varepsilon.$

Using the same argument we show

$\Pr\left(Y_{((1/2+\varepsilon)N)}-F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)\leqslant-t\right)\leqslant\exp\left(-2N(\varphi_{-}(t))^{2}\right).$

We have

$|\varphi_{-}(t)|=F_{\chi_{1}^{2}}(F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon))-F_{\chi_{1}^{2}}(F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)-t)\geqslant\frac{t\exp\left(-F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)\right)}{\sqrt{2\pi
F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)}}.$

Since $\varepsilon\leqslant 1/4$ we conclude the proof in the regime
$t\leqslant F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)$. Observe that due to the
non-negativity of $Y_{((1/2+\varepsilon)N)}$ we can extend this bound to all
$t>F^{-1}_{\chi_{1}^{2}}(1/2+\varepsilon)$. Lemma 6 in Section 4 concludes the
proof. The analysis of the $(1/2-\varepsilon)$-th quantile repeats the same
lines. ∎

Our final result proves a similar bound for the (standard) half-normal
distribution. Namely, we want to prove the concentration inequality for
quantiles of the absolute values of standard Gaussian random variables. Denote
the cumulative distribution function of this distribution by
$\Phi_{\operatorname{H}}(\cdot)$ and its inverse by
$\Phi^{-1}_{\operatorname{H}}(\cdot)$.

######  Lemma 3 (Quantiles of the half-normal distribution).

There is an absolute constant $c_{1}$ such that the following holds. Assume
without loss of generality that $(1/2\pm\varepsilon)N$ are integers. Let
$Y_{1},\ldots,Y_{N}$ be a sample of independent half-normal random variables
and $\varepsilon\in[0,1/4]$. Then,

$\left\|Y_{((1/2\pm\varepsilon)N)}-\Phi^{-1}_{\operatorname{H}}(1/2\pm\varepsilon)\right\|_{\psi_{2}}\leqslant\frac{c_{1}}{\sqrt{N}}.$

The proof of this result repeats the same computations used in the proofs of
Lemma 1 and Lemma 2. We omit the details.

##  4 Auxiliary results

The following section contains several technical results used throughout the
paper. We start with a bound usually referred to as the _PAC-Bayesian lemma_ ,
which is a direct consequence of the Donsker–Varadhan’s variational formula
for the relative entropy [DV75].

######  Lemma 4.

Assume that $X$ is a random variable defined on some measurable space
$\mathcal{X}$. Assume also that $\Theta$ (called the parameter space) is a
subset of $\mathbb{R}^{d}$. Let $\gamma$ be a distribution (called prior) on
$\Theta$ and let $\rho$ be any distribution (called posterior) on $\Theta$
such that $\rho\ll\gamma$. Let $f:\mathcal{X}\times\Theta\to\mathbb{R}$ be
such that ${\mathbf{E}}_{X}\exp(f(X,\theta))$ is finite $\gamma$-almost
surely. Then, we have

$\displaystyle\Pr_{X}\Big{(}\textrm{for all}\
\rho\ll\gamma~{}:~{}{\mathbf{E}}_{\rho}f(X,\theta)\leqslant{\mathbf{E}}_{\rho}\log({\mathbf{E}}_{X}\exp(f(X,\theta)))+\mathcal{KL}(\rho,\gamma)+t\Big{)}\geqslant
1-e^{-t}.$

One of the key arguments, used in several recent papers on mean and covariance
estimation of heavy-tailed distributions [CG17, Giu18, AZ22, OR22], is a
skillful application of this lemma allowing to bypass the sphere-covering and
VC-type arguments. Lemma 4 will play the same key role in our analysis.
However, previous applications of this lemma were based on sums of (truncated)
random variables in place of $f(X,\theta)$ (in this case $X$ is essentially a
vector of independent random variables $X_{1},\ldots,X_{N}$), while we are
exploiting the interplay between Lemma 4 and sample quantiles of particular
univariate distributions. The proof of Lemma 4 and some of its applications
can be found in [CG17, Zhi21].

###  4.1 Analysis of the posterior distribution

Another technical aspect of our analysis is the introduction of truncated
posterior distributions in the context of robust estimation. For a given
positive semi-definite matrix $G$ and $r\geqslant 0$, we truncate the
multivariate Gaussian distribution with mean $v\in S^{d-1}$ and covariance
$\beta^{-1}I_{d}$ as follows. Define the density function

$f_{v}(x)=\frac{1}{p(2\pi\beta^{-1})^{d/2}}\exp\left(-\frac{\beta\|x-v\|^{2}}{2}\right)\bm{1}\left(\|G^{1/2}(x-v)\|\leqslant
r\right),$ (14)

where $p>0$ is a normalization factor. We proceed with the following result.

######  Lemma 5 (Properties of the truncated posterior).

Let $r,\beta>0$ and let $G$ denote a positive semi-definite in the definition
(14). Let $\Sigma$ be a covariance matrix of a zero mean random vector $X$ in
$\mathbb{R}^{d}$ satisfying

$\frac{1}{10}\Sigma\preceq
G,\quad\textrm{and}\quad\operatorname{Tr}(G)\leqslant
10\operatorname{Tr}(\Sigma).$

Let $\gamma$ be a Gaussian measure in $\mathbb{R}^{d}$ with mean zero and
covariance $\beta^{-1}I_{d}$. If, additionally

$r\geqslant\sqrt{20\beta^{-1}\operatorname{Tr}(\Sigma)},$

then we have

$\mathcal{KL}(\rho_{v},\gamma)\leqslant\log(2)+\beta/2.$

Furthermore, let $\theta$ be distributed according to $\rho_{v}$. Then,
${\mathbf{E}}_{\rho_{v}}\theta=v$, and almost surely with respect to the
realization of $\theta$, we have

$\theta^{\mathsf{T}}\Sigma\theta\leqslant 2\|\Sigma\|+20r^{2}.$ (15)

###### Proof.

We use that for $\theta$ distributed according to $\rho_{v}$ it holds that
${\mathbf{E}}_{\rho_{v}}\theta=v$. This follows from the symmetry of the
density around $v$. Let $g$ denote the density of a Gaussian random vector
with mean zero and covariance $\beta^{-1}I_{d}$. To control
$\mathcal{KL}(\rho_{v},\mu)$ we write

$\displaystyle\int\log\left(\frac{f_{v}(x)}{g(x)}\right)f_{v}(x)dx$
$\displaystyle={\mathbf{E}}_{\rho_{v}}\log\left(\frac{1}{p}\exp\left(\frac{-\beta\|\theta-v\|^{2}+\beta\|\theta\|^{2}}{2}\right)\right)$
$\displaystyle=\log\left(\frac{1}{p}\right)+{\mathbf{E}}_{\rho_{v}}\left(\frac{-\beta\|v\|^{2}+2\beta\langle\theta,v\rangle}{2}\right)$
$\displaystyle=\log\left(\frac{1}{p}\right)+\frac{\beta}{2}.$

To prove the desired inequality we observe that

$p=\Pr(\|G^{1/2}W\|\leqslant r),$

where $W$ is a zero mean Gaussian random vector with covariance
$\beta^{-1}I_{d}$. Since $\operatorname{Tr}(G)\leqslant
10\operatorname{Tr}(\Sigma)$, a simple computation shows that

$\Pr(\|G^{1/2}W\|\geqslant
r)\leqslant{\mathbf{E}}W^{\mathsf{T}}GW/r^{2}=\beta^{-1}\operatorname{Tr}(G)/r^{2}\leqslant
10\beta^{-1}\operatorname{Tr}(\Sigma)/r^{2}\leqslant 1/2,$

as long as $r\geqslant\sqrt{20\beta^{-1}\operatorname{Tr}(\Sigma)}$. Thus,
under this assumption $p\geqslant 1/2$, and hence $\log(1/p)\leqslant\log 2$.
This proves the first inequality. Using the second property of the matrix $G$,
we have

$\theta^{\mathsf{T}}\Sigma\theta\leqslant 2v^{\mathsf{T}}\Sigma
v+2(\theta-v)^{\mathsf{T}}\Sigma(\theta-v)\leqslant 2v^{\mathsf{T}}\Sigma
v+20(\theta-v)^{\mathsf{T}}G(\theta-v)\leqslant 2\|\Sigma\|+20r^{2}.$

The claim follows. ∎

Our next result convertes a mixed sub-Gaussian/sub-exponential tail bound into
a bound on the $\|\cdot\|_{\psi_{1}}$-norm. We present this standard
computation for the sake of completeness.

######  Lemma 6.

Assume that a random variable $X$ satisfies for all $t\geqslant 0$,

$\Pr(|X|\geqslant t)\leqslant 2\exp(-K\min\\{t^{2},t\\}),$

where $K>1$ is some constant. Then, there is an absolute constant $c>0$ such
that

$\|X\|_{\psi_{1}}\leqslant\frac{c}{\sqrt{K}}.$

###### Proof.

We can simply compute the moments of $X$. For fixed $p\geqslant 1$, we have

$\displaystyle{\mathbf{E}}|X|^{p}$
$\displaystyle=\int\limits_{0}^{\infty}\Pr(|X|^{p}\geqslant
t)~{}dt=\int\limits_{0}^{\infty}\Pr(|X|\geqslant t)pt^{p-1}~{}dt\leqslant
2\int\limits_{0}^{\infty}\exp(-K\min\\{t^{2},t\\})pt^{p-1}dt$
$\displaystyle\leqslant
2\int\limits_{0}^{\infty}\exp(-Kt^{2})pt^{p-1}dt+2\int\limits_{0}^{\infty}\exp(-Kt)pt^{p-1}dt=\frac{1}{K^{p/2}}p\Gamma(p/2)+\frac{1}{K^{p}}2p\Gamma(p)$
$\displaystyle\leqslant\frac{3p(p/2)^{p/2}}{K^{p/2}}+\frac{2p^{p}}{K^{p}}\leqslant\frac{3p^{p}}{K^{p/2}}+\frac{2p^{p}}{K^{p}}\leqslant\frac{5p^{p}}{K^{p/2}},$

where $\Gamma(\cdot)$ stands for the gamma function, and we used
$\Gamma(x)\leqslant 3x^{x}$ for all $x\geqslant 1/2$ together with
$p\Gamma(p)=\Gamma(p+1)\leqslant p^{p}$. Finally, [Ver16, Proposition 2.7.1,
(b)] implies the desired bound. ∎

##  5 Proofs of main results

We begin with the proof of our first main result that yields that the
estimator defined in (4) achieves an optimal error bound for the robust mean
estimation problem. We discuss the optimality of our results at the end of
this section.

###### Proof of Theorem 1..

First, by the definition of our estimator, we have

$\widehat{\mu}=\mathop{\mathrm{arg}\,\mathrm{min}}_{\nu\in\mathbb{R}^{d}}\sup\limits_{v\in
S^{d-1}}|{\mathbf{E}}_{\rho_{v}}\operatorname{Med}(\langle
X_{1},\theta\rangle,\ldots,\langle X_{N},\theta\rangle)-\langle\nu,v\rangle|,$

where $\rho_{v}$ is a multivariate Gaussian distribution in $\mathbb{R}^{d}$
with mean $v$ and covariance $\beta^{-1}I_{d}$. By the triangle inequality and
the definition of our estimator, we have

$\displaystyle\|\widehat{\mu}-\mu\|$ $\displaystyle=\sup_{v\in
S^{d-1}}\langle\widehat{\mu}-\mu,v\rangle$
$\displaystyle\leqslant\sup\limits_{v\in
S^{d-1}}|{\mathbf{E}}_{\rho_{v}}\operatorname{Med}(\langle
X_{1},\theta\rangle,\ldots,\langle
X_{N},\theta\rangle)-\langle\widehat{\mu},v\rangle|$
$\displaystyle\quad+\sup\limits_{v\in
S^{d-1}}|{\mathbf{E}}_{\rho_{v}}\operatorname{Med}(\langle
X_{1},\theta\rangle,\ldots,\langle X_{N},\theta\rangle)-\langle\mu,v\rangle|$
$\displaystyle\leqslant 2\sup\limits_{v\in
S^{d-1}}|{\mathbf{E}}_{\rho_{v}}\operatorname{Med}(\langle
X_{1},\theta\rangle,\ldots,\langle X_{N},\theta\rangle)-\langle\mu,v\rangle|$
$\displaystyle=2\sup\limits_{v\in
S^{d-1}}|{\mathbf{E}}_{\rho_{v}}\operatorname{Med}(\langle
X_{1}-\mu,\theta\rangle,\ldots,\langle X_{N}-\mu,\theta\rangle)|.$

We need to bound the last quantity. From now on we can assume without loss of
generality that $\mu=0$. Assume $Y_{1},\ldots,Y_{N}$ is an uncontaminated
sample of zero mean independent Gaussians with covariance $\Sigma$. That is,
at most $\varepsilon N$ elements among $X_{1},\ldots,X_{N}$ are different from
their $Y_{1},\ldots,Y_{N}$ counterparts. Observe that the sample median of
projections of the contaminated sample on any direction cannot be too far away
from $1/2\pm\varepsilon$ quantiles of the corresponding projections for the
uncontaminated sample. Formally, assuming that both the sample median and
$1/2\pm\varepsilon$ sample quantiles are unique, we have $\rho_{v}$-almost
surely

$\displaystyle\operatorname{Quant}_{\frac{1}{2}-\varepsilon}(\langle
Y_{1},\theta\rangle,\ldots,\langle Y_{N},\theta\rangle)$
$\displaystyle\leqslant\operatorname{Med}(\langle
X_{1},\theta\rangle,\ldots,\langle X_{N},\theta\rangle)$
$\displaystyle\leqslant\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(\langle
Y_{1},\theta\rangle,\ldots,\langle Y_{N},\theta\rangle),$

and thus, taking the expectation with respect to $\rho_{v}$, we have

$\displaystyle{\mathbf{E}}_{\rho_{v}}\operatorname{Quant}_{\frac{1}{2}-\varepsilon}(\langle
Y_{1},\theta\rangle,\ldots,\langle Y_{N},\theta\rangle)$
$\displaystyle\leqslant{\mathbf{E}}_{\rho_{v}}\operatorname{Med}(\langle
X_{1},\theta\rangle,\ldots,\langle X_{N},\theta\rangle)$
$\displaystyle\leqslant{\mathbf{E}}_{\rho_{v}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(\langle
Y_{1},\theta\rangle,\ldots,\langle Y_{N},\theta\rangle).$

Therefore, we have for any $v\in S^{d-1}$,

$\displaystyle\left|{\mathbf{E}}_{\rho_{v}}\operatorname{Med}(\langle
X_{1},\theta\rangle,\ldots,\langle X_{N},\theta\rangle)\right|$
$\displaystyle\leqslant\left|{\mathbf{E}}_{\rho_{v}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(\langle
Y_{1},\theta\rangle,\ldots,\langle Y_{N},\theta\rangle)\right|$
$\displaystyle\qquad+\left|{\mathbf{E}}_{\rho_{v}}\operatorname{Quant}_{\frac{1}{2}-\varepsilon}(\langle
Y_{1},\theta\rangle,\ldots,\langle Y_{N},\theta\rangle)\right|.$

Both terms will be analyzed similarly. We only analyze the first one. Observe
that due to the spherical symmetry, we have that $S_{N}=\\{\langle
Y_{1},\theta\rangle/\sqrt{\theta^{\mathsf{T}}\Sigma\theta},\ldots,\langle
Y_{N},\theta\rangle/\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\\}$ consists of
independent standard Gaussian random variables (in our case, $\theta\neq 0$
almost surely). We have

$\displaystyle\left|{\mathbf{E}}_{\rho_{v}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(\langle
Y_{1},\theta\rangle,\ldots,\langle Y_{N},\theta\rangle)\right|$
$\displaystyle\leqslant\left|{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left(\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})-{\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})\right)\right|$
$\displaystyle\qquad+\left|{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\cdot({\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})-\Phi^{-1}(1/2+\varepsilon))\right|$
$\displaystyle\qquad+{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\cdot\Phi^{-1}(1/2+\varepsilon)$
$\displaystyle=(I)+(II)+(III).$

To upper bound $(I)$ we want to apply Lemma 4. Fix $\lambda>0$ and let
$\gamma$ be a multivariate Gaussian distribution in $\mathbb{R}^{d}$ with zero
mean and covariance $\beta^{-1}I_{d}$. The standard formula implies
$\mathcal{KL}(\rho_{v},\gamma)=\beta/2$. Thus, by Lemma 4 we have
simultaneously for all $v\in S^{d-1}$, with probability at least $1-\delta$,

$\displaystyle\lambda{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left(\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})-{\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})\right)$
$\displaystyle\quad\leqslant{\mathbf{E}}_{\rho_{v}}\log{\mathbf{E}}\exp\left(\lambda\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left(\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})-{\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})\right)\right)+\beta/2+\log(1/\delta).$

Since centering multiplies the $\psi_{2}$-norm by at most an absolute constant
factor (see e.g., [Ver16, Lemma 2.6.8]), we have by Lemma 1, for some absolute
constant $c_{1}>0$,

$\left\|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})-{\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})\right\|_{\psi_{2}}\leqslant\frac{c_{1}}{\sqrt{N}}.$

Thus, by [Ver16, Proposition 2.5.2, (v)] (conditioned on $\theta$, we take
$\lambda\sqrt{\theta^{\mathsf{T}}\Sigma\theta}$ instead of $\lambda$ in that
result), we have for some absolute constant $c_{2}>0$,

$\displaystyle{\mathbf{E}}_{\rho_{v}}\log{\mathbf{E}}\exp\left(\lambda\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left(\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})-{\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})\right)\right)$
$\displaystyle\leqslant\frac{{\mathbf{E}}_{\rho_{v}}c_{2}\lambda^{2}\theta^{\mathsf{T}}\Sigma\theta}{N}$
$\displaystyle=\frac{c_{2}\lambda^{2}(v^{\mathsf{T}}\Sigma
v+\beta^{-1}\operatorname{Tr}(\Sigma))}{N}$
$\displaystyle\leqslant\frac{11c_{2}\lambda^{2}\|\Sigma\|}{N},$

where the last lines are based on a direct computation and our choice of
$\beta$ (we have $\beta^{-1}\leqslant
10\|\Sigma\|/\operatorname{Tr}(\Sigma)$). Optimizing the bound on
$\lambda{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left(\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})-{\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})\right)$
with respect to $\lambda>0$ and since $\beta\leqslant
10\operatorname{Tr}(\Sigma)/\|\Sigma\|$, we obtain that uniformly over
$S^{d-1}$,

${\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left(\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})-{\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})\right)\leqslant
c_{3}\sqrt{\frac{\operatorname{Tr}(\Sigma)+\|\Sigma\|\log(1/\delta)}{N}},$

where $c_{3}>0$ is an absolute constant. Repeating the proof for $\lambda<0$
and using the union bound, we extend this bound to the upper bound $(I)$.

We now focus on bounding $(II)$. Observe that for any scalar $C$, we have
$\|C\|_{\psi_{2}}=|C|/\sqrt{\log 2}$. Using this observation, together with
Jensen’s inequality and Lemma 1 we have for some $c_{4}>0$,

$\displaystyle\left|{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\cdot({\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}+\varepsilon}(S_{N})-\Phi^{-1}(1/2+\varepsilon))\right|$
$\displaystyle\leqslant\frac{c_{4}\sqrt{\log
2}\cdot{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}}{\sqrt{N}}$
$\displaystyle\leqslant\frac{c_{4}\sqrt{\log
2}\cdot\sqrt{11\|\Sigma\|}}{\sqrt{N}}.$

Finally, we bound $(III)$. First, we notice that the function
$\Phi^{-1}(\cdot)$ is locally Lipschitz on a closed interval
$[1/2-\varepsilon,1/2+\varepsilon]$ for $\varepsilon\in[0,1/4]$. We compute
and bound the local Lipschitz constant of $\Phi^{-1}$ as follows

$\frac{d}{dx}\Phi^{-1}(1/2+x)=\sqrt{2\pi}\exp\Big{(}(\Phi^{-1}(1/2+x))^{2}/{2}\Big{)}\leqslant\sqrt{2\pi}\exp\Big{(}(\Phi^{-1}(3/4))^{2}/{2}\Big{)}\leqslant
4,$

where we used the fact that the function $\exp(\Phi^{-1}(\cdot)^{2})$ is
increasing. Since the standard Gaussian distribution is symmetric, we have
$\Phi^{-1}(1/2)=0$. Hence, the bound for $\Phi^{-1}(1/2+\varepsilon)$ reads as

$\Phi^{-1}(1/2+\varepsilon)=\Phi^{-1}(1/2+\varepsilon)-\Phi^{-1}(1/2)\leqslant
4\varepsilon.$

Therefore, we have

${\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\cdot\Phi^{-1}(1/2+\varepsilon)\leqslant
4\varepsilon\sqrt{11\|\Sigma\|}.$

Combining the above bounds concludes the proof. ∎

We are now ready to prove our second main result.

###### Proof of Theorem 2.

Recall that $H$ is a covariance matrix of $\rho_{v}$ and does not depend on
direction $v\in S^{d-1}$. Observe that
${\mathbf{E}}_{\rho_{v}}\theta^{\mathsf{T}}\Sigma\theta=v^{\mathsf{T}}\Sigma
v+\operatorname{Tr}(\Sigma H)$. Moreover, since our choice of parameters
implies $p\geqslant 1/2$ in (14), we have

$H\preceq 2\beta^{-1}I_{d},\quad\textrm{and}\quad\|H\|\leqslant 2\beta^{-1}.$

We also observe that $\mathbf{r}(G)=\operatorname{Tr}(G)/\|G\|\leqslant
100\mathbf{r}(\Sigma)$. Using the triangle inequality, as well as the
definition of our estimator combined with the definition of the set
$\mathcal{H}$ from (9), we have

$\displaystyle\|\widehat{\Sigma}-\Sigma\|$ $\displaystyle=\sup\limits_{v\in
S^{d-1}}\left|v^{\mathsf{T}}{\Sigma}v-v^{\mathsf{T}}\widehat{\Sigma}v\right|$
$\displaystyle\leqslant\sup\limits_{v\in
S^{d-1}}\left|{\mathbf{E}}_{\rho_{v}}\theta^{\mathsf{T}}({\Sigma}-\widehat{\Sigma})\theta\right|+|(\operatorname{Tr}(\widehat{\Sigma}H)-\alpha(H))-(\operatorname{Tr}({\Sigma}H)-\alpha(H))|$
$\displaystyle\leqslant\sup\limits_{v\in
S^{d-1}}\left|{\mathbf{E}}_{\rho_{v}}\theta^{\mathsf{T}}({\Sigma}-\widehat{\Sigma})\theta\right|+c\operatorname{Tr}((\Sigma+\widehat{\Sigma})H)\left(\sqrt{\frac{\max\\{\mathbf{r}(G),\mathbf{r}(\Sigma)\\}+\log(1/\delta)}{N}}+\varepsilon\right)$
$\displaystyle\leqslant\sup\limits_{v\in
S^{d-1}}\left|{\mathbf{E}}_{\rho_{v}}\theta^{\mathsf{T}}({\Sigma}-\widehat{\Sigma})\theta\right|+2\beta^{-1}c\operatorname{Tr}(\Sigma+10G)\left(\sqrt{\frac{100\mathbf{r}(\Sigma)+\log(1/\delta)}{N}}+\varepsilon\right)$
$\displaystyle\leqslant\sup\limits_{v\in
S^{d-1}}\left|{\mathbf{E}}_{\rho_{v}}\theta^{\mathsf{T}}({\Sigma}-\widehat{\Sigma})\theta\right|+202\beta^{-1}c\operatorname{Tr}(\Sigma)\left(\sqrt{\frac{100\mathbf{r}(\Sigma)+\log(1/\delta)}{N}}+\varepsilon\right).$

Since $\beta^{-1}\leqslant 10\mathbf{r}(\Sigma)$, the last term in the last
inequality is not larger than the rate of convergence in the statement of
Theorem 2. We now can focus only on bounding the first term in the last line
of the inequalities from the previous display. We first need some auxiliary
computations. Using the definition of the set $\mathcal{H}$, we have
$\rho_{v}$-almost surely

$\displaystyle\sqrt{\theta^{\mathsf{T}}{\Sigma}\theta}+\sqrt{\theta^{\mathsf{T}}\widehat{\Sigma}\theta}$
$\displaystyle\leqslant\sqrt{2v^{\mathsf{T}}{\Sigma}v+2(\theta-v)^{\mathsf{T}}{\Sigma}(\theta-v)}+\sqrt{2v^{\mathsf{T}}\widehat{\Sigma}v+2(\theta-v)^{\mathsf{T}}\widehat{\Sigma}(\theta-v)}$
$\displaystyle\leqslant\sqrt{2\|\Sigma\|+20(\theta-v)^{\mathsf{T}}{G}(\theta-v)}+\sqrt{20\omega+20(\theta-v)^{\mathsf{T}}G(\theta-v)}$
$\displaystyle\leqslant\sqrt{2\|\Sigma\|+20r^{2}}+\sqrt{20\omega+20r^{2}}$
$\displaystyle\leqslant c_{1}\sqrt{\|\Sigma\|},$

where $c_{1}>0$ is some absolute constant. This implies the following lines

$\displaystyle\sup\limits_{v\in
S^{d-1}}\left|{\mathbf{E}}_{\rho_{v}}\theta^{\mathsf{T}}({\Sigma}-\widehat{\Sigma})\theta\right|$
$\displaystyle\leqslant\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\left|\sqrt{\theta^{\mathsf{T}}{\Sigma}\theta}-\sqrt{\theta^{\mathsf{T}}\widehat{\Sigma}\theta}\right|\left|\sqrt{\theta^{\mathsf{T}}{\Sigma}\theta}+\sqrt{\theta^{\mathsf{T}}\widehat{\Sigma}\theta}\right|$
$\displaystyle\leqslant c_{1}\sqrt{\|\Sigma\|}\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\left|\sqrt{\theta^{\mathsf{T}}{\Sigma}\theta}-\sqrt{\theta^{\mathsf{T}}\widehat{\Sigma}\theta}\right|$
$\displaystyle\leqslant c_{1}\sqrt{\|\Sigma\|}\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\left|\operatorname{Med}\left(|\langle
X_{1},\theta\rangle|,\ldots,|\langle
X_{N},\theta\rangle|\right)/(\Phi^{-1}(3/4))-\sqrt{\theta^{\mathsf{T}}\widehat{\Sigma}\theta}\right|$
$\displaystyle\qquad+c_{1}\sqrt{\|\Sigma\|}\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\left|\operatorname{Med}\left(|\langle
X_{1},\theta\rangle|,\ldots,|\langle
X_{N},\theta\rangle|\right)/(\Phi^{-1}(3/4))-\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\right|$
$\displaystyle\leqslant 2c_{1}\sqrt{\|\Sigma\|}\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\left|\operatorname{Med}\left(|\langle
X_{1},\theta\rangle|,\ldots,|\langle
X_{N},\theta\rangle|\right)/(\Phi^{-1}(3/4))-\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\right|,$

where in the last line we used the definition of $\widehat{\Sigma}$ and that
$\Sigma\in\mathcal{H}$. We focus on upper bounding the last expression. Let
$Y_{1},\ldots,Y_{N}$ be the uncontaminated version of our
$\varepsilon$-contaminated sample. Using the same argument as in the proof of
Theorem 1, we have

$\displaystyle\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\left|\operatorname{Med}\left(|\langle
X_{1},\theta\rangle|,\ldots,|\langle
X_{N},\theta\rangle|\right)-\Phi^{-1}(3/4)\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\right|$
$\displaystyle\leqslant\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\left|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}\left(|\langle
Y_{1},\theta\rangle|,\ldots,|\langle
Y_{N},\theta\rangle|\right)-\Phi^{-1}(3/4)\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\right|$
$\displaystyle\qquad+\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\left|\operatorname{Quant}_{\frac{1}{2}-\varepsilon}\left(|\langle
Y_{1},\theta\rangle|,\ldots,|\langle
Y_{N},\theta\rangle|\right)-\Phi^{-1}(3/4)\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\right|.$

We only analyze the first term. Observe that due to the spherical symmetry, we
have that $S_{N}=\\{|\langle
Y_{1},\theta\rangle|/\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\ ,\ldots,\
|\langle Y_{N},\theta\rangle|/\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\\}$
consists of independent half-normal random variables (in our case, $\theta\neq
0$ almost surely). By the triangle inequality, we have

$\displaystyle\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}\left(S_{N}\right)-\Phi^{-1}(3/4)\right|$
$\displaystyle\leqslant\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left(\left|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}\left(S_{N}\right)-\Phi^{-1}_{\operatorname{H}}(1/2+\varepsilon)\right|-{\mathbf{E}}\left|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}\left(S_{N}\right)-\Phi^{-1}_{\operatorname{H}}(1/2+\varepsilon)\right|\right)$
$\displaystyle\qquad+\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}{\mathbf{E}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}\left(S_{N}\right)-\Phi^{-1}_{\operatorname{H}}(1/2+\varepsilon)\right|$
$\displaystyle\qquad+\sup\limits_{v\in
S^{d-1}}{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left|\Phi^{-1}(3/4)-\Phi^{-1}_{\operatorname{H}}(1/2+\varepsilon)\right|$
$\displaystyle=(I)+(II)+(III).$

We want to apply Lemma 4 to control $(I)$. Fix $\lambda>0$ and let $\gamma$ be
a multivariate Gaussian distribution in $\mathbb{R}^{d}$ with zero mean and
covariance $\beta^{-1}I_{d}$. Lemma 5 implies that for our choice of
parameters $\mathcal{KL}(\rho_{v},\gamma)\leqslant\log(2)+\beta/2$. Denote

$Q(S_{N})=\left|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}\left(S_{N}\right)-\Phi^{-1}_{\operatorname{H}}(1/2+\varepsilon)\right|-{\mathbf{E}}\left|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}\left(S_{N}\right)-\Phi^{-1}_{\operatorname{H}}(1/2+\varepsilon)\right|.$

Observe that conditioned on $\theta$, the random variable $Q(S_{N})$ is a
centered version of the random variable
$|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}\left(S_{N}\right)-\Phi^{-1}_{\operatorname{H}}(1/2+\varepsilon)|$
whose $\|\cdot\|_{\psi_{2}}$ is controlled by Lemma 3. Since centering
multiplies the $\psi_{2}$-norm by at most an absolute constant factor, we have
(conditioned on $\theta$) that
$\left\|Q(S_{N})\right\|_{\psi_{2}}\leqslant\frac{c_{2}}{\sqrt{N}}$ for some
absolute constant $c_{2}>0$. By Lemma 4 we have, simultaneously for all $v\in
S^{d-1}$, with probability at least $1-\delta$,

$\displaystyle\lambda{\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\
Q(S_{N})\leqslant{\mathbf{E}}_{\rho_{v}}\log{\mathbf{E}}\exp\left(\lambda\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\
Q(S_{N})\right)+\beta/2+\log(2/\delta).$

Thus, by [Ver16, Proposition 2.5.2, (v)] (conditioned on $\theta$, we take
$\lambda\sqrt{\theta^{\mathsf{T}}\Sigma\theta}$ instead of $\lambda$ in that
result), repeating the lines of the proof of Theorem 1, we have for some
absolute constants $c_{3},c_{4}>0$,

${\mathbf{E}}_{\rho_{v}}\log{\mathbf{E}}\exp\Big{(}\lambda\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\
Q(S_{N})\Big{)}\leqslant{\mathbf{E}}_{\rho_{v}}\frac{c_{3}\lambda^{2}\theta^{\mathsf{T}}\Sigma\theta}{N}\leqslant\frac{c_{4}\lambda^{2}\|\Sigma\|}{N}.$

Combining the bounds and optimizing with respect to $\lambda$, we have
simultaneously for all $v\in S^{d-1}$, with probability at least $1-\delta$,

$(I)\leqslant
c_{5}{\|\Sigma\|^{1/2}}\Bigg{(}\sqrt{\frac{\mathbf{r}(\Sigma)}{N}}+\sqrt{\frac{\log(1/\delta)}{N}}\Bigg{)},$

where $c_{5}>0$ is some absolute constant. We now bound the term $(II)$.
Similarly to the proof of Theorem 1, we use Lemma 3 to get, for some absolute
constant $c_{6}>0$, the following bound

${\mathbf{E}}_{\rho_{v}}{\mathbf{E}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left|\operatorname{Quant}_{\frac{1}{2}+\varepsilon}\left(S_{N}\right)-\Phi^{-1}_{\operatorname{H}}(1/2+\varepsilon)\right|\leqslant
c_{6}\sqrt{\frac{\|\Sigma\|}{N}}.$

To bound $(III)$ we first observe that
$\Phi^{-1}_{\operatorname{H}}(1/2)=\Phi^{-1}(3/4)$. Now we show that the
difference
$\Phi_{\operatorname{H}}^{-1}(1/2+\varepsilon)-\Phi^{-1}_{\operatorname{H}}(1/2)$
is bounded by $\varepsilon$ (up to multiplicative constant) for
$\varepsilon\in[0,1/4]$. Similarly to the arguments used in the proof of
Theorem 1 for the quantile function of standard Gaussian distribution, we
compute and bound the derivative of $\Phi_{\operatorname{H}}^{-1}(1/2+x)$ when
$x\in[0,1/4]$ as follows

$\frac{d}{dx}\Phi_{\operatorname{H}}^{-1}(1/2+x)=\sqrt{\frac{\pi}{2}}\exp\Big{(}(\Phi^{-1}_{\operatorname{H}}(1/2+x))^{2}/{2}\Big{)}\leqslant\sqrt{\frac{\pi}{2}}\exp\Big{(}(\Phi^{-1}_{\operatorname{H}}(3/4))^{2}/{2}\Big{)}\leqslant
3.$

Therefore, we have for some $c_{7}>0$,

${\mathbf{E}}_{\rho_{v}}\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\cdot\big{|}\Phi_{\operatorname{H}}^{-1}(1/2+\varepsilon)-\Phi_{\operatorname{H}}^{-1}(1/2)\big{|}\leqslant
c_{7}\varepsilon\sqrt{\|\Sigma\|}.$

Combining the obtained bounds, we complete the proof. ∎

#### Statistical optimality of our estimators.

We shortly discuss the claimed optimality of our bounds. The optimality
results follow immediately from existing lower bounds. The bounds in [CGR18,
Theorem 2.2 and Theorem 3.2] show that Theorem 1 and Theorem 2 both have the
optimal dependence on the contamination level with correct dimension-free
parametric rate. For covariance estimation, the optimality of the remaining
terms is described in detail in [AZ22, Section 5]. Matching lower bounds for
the mean estimation problem are shown in [LM19b].

##  6 Tuning the unknown parameters

Our focus is now on tuning a few parameters used in our estimators. For the
sake of simplicity, we assume that either $\varepsilon$ is known exactly or a
known upper bound $\varepsilon_{0}$ is available, such that
$\varepsilon\leqslant\varepsilon_{0}<1/2$. This is a standard assumption in
the literature [LM21]. Observe that at least in mean estimation the value of
$\varepsilon$ is only used to tune the parameter $\beta$. The most standard
approach to estimating other parameters is the sample-splitting idea. One
splits the sample into several independent blocks of equal sizes. For each
block, we can bound the number of contaminated points. This will allow us to
state our result for any $\varepsilon\in[0,c]$, where $c$ is some small enough
absolute constant. An interesting aspect of our analysis is that we can tune
different parameters on the same sample. We will now discuss this in more
details.

#### Handling the dependencies.

It is clear that in the strong contamination setup, the adversary can make the
aforementioned blocks dependent. That is, the outliers in any sub-sample may
depend on the entire sample. Some authors assume implicitly that the splitting
of the sample results in independent subsamples. For example, the analysis of
the trimmed-mean estimator in [LM21, Theorem 1] uses this independence, which
holds, for example, in Huber’s contamination model, but is not true in the
general strong contamination model. Taking care of the sample splitting step
in the strong contamination model requires some additional stability-type
analysis. We refer to [DKLP22, Section 6] as an example of this approach.

We now show that our approach allows one to tune the parameters on the same
sample. Thus, our result is valid in the strong contamination model without
additional assumptions. For clarity, we only focus on the mean estimation
problem. Assume we are given an $\varepsilon$-contaminated sample of size $N$.
We denote it by $S_{N}$. Given $S_{N}$, we first find an integer
$\beta=\beta(S_{N})$ satisfying, with probability at least $1-\delta/2$,

$\mathbf{r}(\Sigma)/10\leqslant\beta(S_{N})\leqslant 10\mathbf{r}(\Sigma).$
(16)

We then compute our estimator defined in (4) on the same sample $S_{N}$ with
$\beta=\beta(S_{N})$. Denote the event where (16) holds by $E$. We show that
due to the nature of Lemma 4, this dependence does not lead to additional
technical issues. First, observe that since $\beta$ is an integer, we can use
the union bound over at most $10\mathbf{r}(\Sigma)$ prior Gaussian
distributions $\gamma$ to handle potential dependence of $\beta$ on $S_{N}$.
One can verify that this application of the union bound does not change the
bound of Theorem 1. Importantly, the result of Lemma 4 is uniform with respect
to the posterior distribution $\rho_{v}$ and allows $\beta$ to depend on the
sample as long as $\mathcal{KL}(\rho_{v},\gamma)=\beta(S_{N})/2\leqslant
5\mathbf{r}(\Sigma)$, which holds on the event $E$. Finally, one can easily
verify that, on the same event $E$, the desired upper bound on the term

${\mathbf{E}}_{\rho_{v}}\log{\mathbf{E}}\exp\left(\lambda\sqrt{\theta^{\mathsf{T}}\Sigma\theta}\left(\operatorname{Quant}_{\frac{1}{2}\pm\varepsilon}(S_{N})-{\mathbf{E}}\operatorname{Quant}_{\frac{1}{2}\pm\varepsilon}(S_{N})\right)\right),$

appearing in the proof of Theorem 1 is not affected by the fact that
$\beta=\beta(S_{N})$. This argument allows us to use $\beta(S_{N})$ in our
estimator.

Similar ideas can also be applied in the covariance estimation setup. To avoid
unnecessary technicalities, we assume that for covariance estimation we can
indeed split the sample into several blocks and the adversary is not allowed
to create dependencies between these blocks. This covers many standard
contamination models, including Huber’s $\varepsilon$-contamination model.

#### Estimating $\beta$ and $\omega$.

This step follows from existing results. In particular, in the Gaussian case
Proposition 6 in [AZ22] provides an estimator $\omega$ satisfying
$\|\Sigma\|/4\leqslant\omega\leqslant 4\|\Sigma\|$ whenever $N\geqslant
c(\mathbf{r}(\Sigma)+\log(1/\delta))$, where $c>0$ is an absolute constant. We
also need to estimate $\operatorname{Tr}(\Sigma)$. This problem reduces to
mean estimation. The linear dependence on $\varepsilon$ will not play any role
since we only need to know $\operatorname{Tr}(\Sigma)$ up to a multiplicative
constant factor. In particular, one can use any _sub-Gaussian mean estimator_
in $\mathbb{R}$ (see [LM19a] for the exact definition) that is tolerant to
strong contamination and gives a $\sqrt{\varepsilon}$-dependence on the
contamination level to find $\tau$ satisfying
$\operatorname{Tr}(\Sigma)/2\leqslant\tau\leqslant
2\operatorname{Tr}(\Sigma)$, whenever $\varepsilon$ is small enough and
$N\geqslant c\log(1/\delta)$. This allows us to find an integer $\beta$
satisfying (16).

#### Constructing the matrix $G$.

We discuss how to construct a positive semi-definite matrix $G$, satisfying

$\Sigma\preceq 10G,\quad\textrm{and}\quad\operatorname{Tr}(G)\leqslant
10\operatorname{Tr}(\Sigma).$ (17)

The following result allows us to construct such a matrix efficiently whenever
$N\geqslant c(d+\log(1/\delta))$, where $c>0$ is some absolute constant.

######  Proposition 1.

There are absolute constants $c,c_{1}>0$ such that the following holds. Assume
that $X$ is a Gaussian zero mean vector in $\mathbb{R}^{d}$ with covariance
$\Sigma$. Let $X_{1},\ldots,X_{N}$ be an $\varepsilon$-contaminated set of
independent copies of $X$. Fix $\delta\in(0,1)$. Assume that
$\varepsilon\leqslant c$ and $N\geqslant c_{1}(d+\log(1/\delta))$. Then, with
probability at least $1-\delta$, simultaneously for all
$I^{\prime}\subseteq[N]$ such that $|I^{\prime}|=cN$, we have

$\Sigma\preceq\frac{10}{N}\sum\limits_{i\in[N]\setminus
I^{\prime}}X_{i}X_{i}^{\mathsf{T}}.$

Moreover, on the same event, there exists $I\subseteq[N]$ such that $|I|=cN$,
and

$\frac{1}{N}\sum\limits_{i\in[N]\setminus I}\|X_{i}\|^{2}\leqslant
10\operatorname{Tr}(\Sigma).$

This result implies immediately that the matrix
$G=\frac{1}{N}\sum\nolimits_{i\in[N]\setminus I}X_{i}X_{i}^{\mathsf{T}}$
satisfies the desired property (17). In order to find this set, one only needs
to find a set $I$ of size $\varepsilon N$ such that
$\sum\nolimits_{i\in[N]\setminus I}\|X_{i}\|^{2}\leqslant
10N\operatorname{Tr}(\Sigma)$. This can be done simply by removing the
$\varepsilon N$ vectors with the largest norms.

###### Proof.

Without loss of generality, we assume that $cN$ is an integer. Fix any
$I\subset[N]$ of size $2cN$. Let $Y_{1},\ldots,Y_{N}$ denote an uncontaminated
sample. The total number of such subsets is upper bounded by
$\binom{N}{2cn}\leqslant\left(\frac{2e}{c}\right)^{2cN}$. By the bound of
Oliveira [Oli16, Theorem 4.1 with $\mathrm{h}=3$] and the union bound over all
sets $I$ of size $2cN$, we have

$\Sigma\left(1-27\sqrt{\frac{d+4cN\log(\frac{2e}{c})+2\log(2/\delta)}{N-2cN}}\right)\preceq\frac{1}{N-2cN}\sum\limits_{i\in[N]\setminus
I}Y_{i}Y_{i}^{\mathsf{T}}.$

When $c$ is small enough and $N\geqslant c_{1}(d+\log(1/\delta))$ for large
enough $c_{1}>0$, on the same event, we have

$\Sigma\preceq\frac{10}{N}\sum\limits_{i\in[N]\setminus
I}Y_{i}Y_{i}^{\mathsf{T}}.$

Observe that since each term $Y_{i}Y_{i}^{\mathsf{T}}$ is a positive semi-
definite matrix and $\varepsilon\leqslant c$, we have that for any
$I^{\prime}$ of size $cN$, there is a set $I$ of size $2cN$ such that

$\sum\nolimits_{i\in[N]\setminus
I}Y_{i}Y_{i}^{\mathsf{T}}\preceq\sum\nolimits_{i\in[N]\setminus
I^{\prime}}X_{i}X_{i}^{\mathsf{T}}.$

Indeed, to build such a set $I$ we consider the union of the set of
contaminated points with the set $I^{\prime}$ (we can add any additional
elements if the cardinality of this union is less than $2cN$). This implies
that under our assumption for all $I^{\prime}\subset[N]$ of size $cN$, with
probability at least $1-\delta$,

$\Sigma\preceq\frac{10}{N}\sum\nolimits_{i\in[N]\setminus
I^{\prime}}X_{i}X_{i}^{\mathsf{T}}.$

We now consider the second part of the statement. Combining the Gaussian
concentration inequality [BLM13, Example 5.7] and [Ver16, Proposition 2.5.2],
we get that there is an absolute constant $c_{2}>0$ such that

$\bigl{\|}\|X\|-{\mathbf{E}}\|X\|\bigr{\|}_{\psi_{2}}\leqslant
c_{2}\sqrt{\|\Sigma\|}.$

It is now standard to verify that
$\bigl{\|}\|X\|^{2}-{\mathbf{E}}\|X\|^{2}\bigr{\|}_{\psi_{1}}\leqslant
c_{3}\|\Sigma\|$, where $c_{3}>0$ is an absolute constant. By the Bernstein
inequality [Ver16, Theorem 2.8.1] and the union bound, simultaneously for all
$I\subset[N],|I|=cN$, with probability at least $1-\delta$, it holds for some
absolute constant $c_{4}>0$ that

$\displaystyle\sum\limits_{i\in[N]\setminus I}\|Y_{i}\|^{2}$
$\displaystyle\leqslant
N\operatorname{Tr}(\Sigma)+c_{4}\|\Sigma\|\left(\sqrt{N(\log(1/\delta)+cN\log(e/c))}+\log(1/\delta)+cN\log(e/c)\right)$
$\displaystyle\leqslant 10N\operatorname{Tr}(\Sigma).$

The last inequality holds provided that $c$ is small enough and $c_{1}$ is
large enough. We choose $I$ to be the set corresponding to the set of
contaminated points. For this set $I$, on the same event, we have

$\frac{1}{N}\sum\limits_{i\in[N]\setminus
I}\|X_{i}\|^{2}=\frac{1}{N}\sum\limits_{i\in[N]\setminus
I}\|Y_{i}\|^{2}\leqslant 10\operatorname{Tr}(\Sigma).$

The claim follows by the union bound. ∎

#### Estimating $\alpha(H)$.

We conclude by the analysis of a real-valued parameter $\alpha=\alpha(H)$,
defined in (7). In what follows, $H$ is a known positive semi-definite matrix.
When allowing slightly sub-optimal dependence on $\varepsilon$, we can use the
analysis of the trimmed mean estimator in $\mathbb{R}$ (see [LM21, Theorem
1]). Unfortunately, the analysis becomes more complicated when the linear
dependence on the contamination level is of interest. Recall that we are
interested in finding $\alpha=\alpha(H)$ such that, with probability at least
$1-\delta$,

$|\alpha-\operatorname{Tr}(\Sigma H)|\leqslant c\operatorname{Tr}(\Sigma
H)\left(\sqrt{\frac{\mathbf{r}(\Sigma)+\log(1/\delta)}{N}}+\varepsilon\right).$

We present an estimator that achieves this error rate in almost any
interesting regime. More precisely, we will either make an additional
assumption that $\delta\geqslant\exp(-\sqrt{\mathbf{r}(\Sigma)})$, or that
$\log d\leqslant\mathbf{r}(\Sigma)$. In what follows, $e_{1},\ldots,e_{d}$
denotes the standard basis in $\mathbb{R}^{d}$.

######  Proposition 2.

There are absolute constants $c,c_{1},c_{2}>0$ such that the following holds.
Assume that $X$ is a Gaussian zero mean vector in $\mathbb{R}^{d}$ with
covariance $\Sigma$. Let $X_{1},\ldots,X_{N}$ be an $\varepsilon$-contaminated
set of independent copies of $X$. Fix $\delta\in(0,1)$. Assume that
$\varepsilon\leqslant c$ and $N\geqslant c_{1}\log(1/\delta)$. Then, with
probability at least $1-\delta$, it holds

$\displaystyle\left|(\Phi^{-1}(3/4))^{-2}\sum\limits_{i=1}^{d}\operatorname{Med}\left(\langle
e_{i},H^{1/2}X_{1}\rangle^{2},\ldots,\langle
e_{i},H^{1/2}X_{N}\rangle^{2}\right)-\operatorname{Tr}(\Sigma H)\right|$
$\displaystyle\qquad\qquad\leqslant c_{1}\operatorname{Tr}(\Sigma
H)\left(\frac{\log(1/\delta)}{\sqrt{N}}+\varepsilon\right).$ (18)

If additionally, $N\geqslant c_{2}(\log d+\log(1/\delta))$, then on the same
event, it holds

$\displaystyle\left|(\Phi^{-1}(3/4))^{-2}\sum\limits_{i=1}^{d}\operatorname{Med}\left(\langle
e_{i},H^{1/2}X_{1}\rangle^{2},\ldots,\langle
e_{i},H^{1/2}X_{N}\rangle^{2}\right)-\operatorname{Tr}(\Sigma H)\right|$
$\displaystyle\qquad\qquad\leqslant c_{1}\operatorname{Tr}(\Sigma
H)\left(\sqrt{\frac{\log d+\log(1/\delta)}{N}}+\varepsilon\right).$

###### Proof.

Let $Y_{1},\ldots,Y_{N}$ denote the uncontaminated sample, and let $Y$ be a
zero mean Gaussian in $\mathbb{R}^{d}$ with covariance $\Sigma$. Since
$\operatorname{Tr}(\Sigma H)=\operatorname{Tr}(H^{1/2}\Sigma H^{1/2})$, by
triangle inequality and the arguments of the proof of Theorem 1, we have

$\displaystyle\left|\sum\limits_{i=1}^{d}\operatorname{Med}\left(\langle
e_{i},H^{1/2}X_{1}\rangle^{2},\ldots,\langle
e_{i},H^{1/2}X_{N}\rangle^{2}\right)-(\Phi^{-1}(3/4))^{2}\operatorname{Tr}(\Sigma
H)\right|$
$\displaystyle\leqslant\sum\limits_{i=1}^{d}\left|\operatorname{Med}\left(\langle
e_{i},H^{1/2}X_{1}\rangle^{2},\ldots,\langle
e_{i},H^{1/2}X_{N}\rangle^{2}\right)-(\Phi^{-1}(3/4))^{2}\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\right|$
$\displaystyle\leqslant\sum\limits_{i=1}^{d}\left|\operatorname{Quant}_{1/2+\varepsilon}\left(\langle
e_{i},H^{1/2}Y_{1}\rangle^{2},\ldots,\langle
e_{i},H^{1/2}Y_{N}\rangle^{2}\right)-(\Phi^{-1}(3/4))^{2}\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\right|$
$\displaystyle\qquad+\sum\limits_{i=1}^{d}\left|\operatorname{Quant}_{1/2-\varepsilon}\left(\langle
e_{i},H^{1/2}Y_{1}\rangle^{2},\ldots,\langle
e_{i},H^{1/2}Y_{N}\rangle^{2}\right)-(\Phi^{-1}(3/4))^{2}\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\right|.$

We only consider the first sum, as the second sum is analyzed similarly.
Observe that by the spherical symmetry the random variable $\langle
e_{i},H^{1/2}Y\rangle^{2}/\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}$ is distributed
according to the $\chi^{2}_{1}$ distribution. Denote

$S_{N,i}=\\{\langle
e_{i},H^{1/2}Y_{1}\rangle^{2}/\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2},\ldots,\langle
e_{i},H^{1/2}Y_{N}\rangle^{2}/\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\\}.$

Using the notation from the previous display, triangle inequality and the fact
that $F^{-1}_{\chi^{2}_{1}}(1/2)=(\Phi^{-1}(3/4))^{2}$, we arrive at

$\displaystyle\left|\operatorname{Quant}_{1/2+\varepsilon}\left(\langle
e_{i},H^{1/2}Y_{1}\rangle^{2},\ldots,\langle
e_{i},H^{1/2}Y_{N}\rangle^{2}\right)-(\Phi^{-1}(3/4))^{2}\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\right|$
$\displaystyle\leqslant\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\left|\operatorname{Quant}_{1/2+\varepsilon}\left(S_{N,i}\right)-{\mathbf{E}}\operatorname{Quant}_{1/2+\varepsilon}\left(S_{N,i}\right)\right|$
$\displaystyle\qquad+\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\left|{\mathbf{E}}\operatorname{Quant}_{1/2+\varepsilon}\left(S_{N,i}\right)-F^{-1}_{\chi^{2}_{1}}(1/2+\varepsilon)\right|$
$\displaystyle\qquad+\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\left|F^{-1}_{\chi^{2}_{1}}(1/2+\varepsilon)-F^{-1}_{\chi^{2}_{1}}(1/2)\right|$
$\displaystyle=(I)_{i}+(II)_{i}+(III)_{i}.$

By Lemma 2 we have for some $c_{1}>0$,

$\left\|(I)_{i}+(II)_{i}\right\|_{\psi_{1}}\leqslant\frac{c_{1}\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}}{\sqrt{N}},\quad\textrm{and
therefore,}\quad\left\|\sum\limits_{i=1}^{d}\bigl{(}(I)_{i}+(II)_{i}\bigr{)}\right\|_{\psi_{1}}\leqslant\frac{c_{1}\operatorname{Tr}(\Sigma
H)}{\sqrt{N}},$

where the last expression follows from the triangle inequality. Using the
exact form of the inverse cumulative distribution function of the
$\chi^{2}_{1}$ distribution and the same technique used to bound the
difference of quantiles of half-normal distribution, one can verify that for
any $\varepsilon\leqslant 1/4$ we have
$\left|F^{-1}_{\chi^{2}_{1}}(1/2\pm\varepsilon)-F^{-1}_{\chi^{2}_{1}}(1/2)\right|\leqslant
c_{2}\varepsilon$, where $c_{2}>0$ is an absolute constant. This readily
yields

$(III)_{i}=\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\left|F^{-1}_{\chi^{2}_{1}}(1/2+\varepsilon)-F^{-1}_{\chi^{2}_{1}}(1/2)\right|\leqslant
c_{2}\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\varepsilon.$

Therefore, for $\varepsilon\leqslant 1/4$, we have
$\sum\nolimits_{i=1}^{d}(III)_{i}\leqslant
c_{2}\varepsilon\operatorname{Tr}(\Sigma H)$. Combing the above computations
and using the tail bound of [Ver16, Proposition 2.7.1], we prove the
inequality (18).

To prove the second part of the bound we propose a slightly different analysis
for the term $(I)_{i}$. Denote

$S^{\prime}_{N,i}=\\{|\langle
e_{i},H^{1/2}Y_{1}\rangle|/\|\Sigma^{1/2}H^{1/2}e_{i}\|,\ldots,|\langle
e_{i},H^{1/2}Y_{N}\rangle|/\|\Sigma^{1/2}H^{1/2}e_{i}\|\\},$

and observe that $S^{\prime}_{N,i}$ consists of independent half-normal random
variables. We have

$\displaystyle(I)_{i}$
$\displaystyle=\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\left|\operatorname{Quant}_{1/2+\varepsilon}\left(S^{\prime}_{N,i}\right)-\sqrt{{\mathbf{E}}\operatorname{Quant}_{1/2+\varepsilon}\left(S_{N,i}\right)}\right|$
$\displaystyle\qquad\times\left|\operatorname{Quant}_{1/2+\varepsilon}\left(S^{\prime}_{N,i}\right)+\sqrt{{\mathbf{E}}\operatorname{Quant}_{1/2+\varepsilon}\left(S_{N,i}\right)}\right|.$

We first bound the second multiplier of the expression from the last display.
By Lemma 3, with probability at least $1-\delta$, we have

$\displaystyle\operatorname{Quant}_{1/2+\varepsilon}(S^{\prime}_{N,i})+\sqrt{{\mathbf{E}}\operatorname{Quant}_{1/2+\varepsilon}\left(S_{N,i}\right)}$
$\displaystyle\leqslant
2\Phi_{\operatorname{H}}^{-1}(1/2+\varepsilon)+c_{3}\sqrt{\frac{\log(1/\delta)}{N}}.$
(19)

Now observe that the last expression from the previous line is bounded by some
absolute constant given that $\varepsilon\in[0,1/4]$ and $N\geqslant
c_{4}\log(1/\delta)$. Using Lemma 3 once again together with union bound, we
bound the term $(I)_{i}$, with probability at least $1-\delta$, as follows

$\displaystyle(I)_{i}\leqslant
c_{4}\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\sqrt{\frac{\log(1/\delta)}{N}}.$ (20)

By the union bound, for all $i\in[d]$ we have, with probability at least
$1-\delta$,

$\displaystyle(I)_{i}\leqslant
c_{3}\|\Sigma^{1/2}H^{1/2}e_{i}\|^{2}\sqrt{\frac{\log d+\log(1/\delta)}{N}},$
(21)

whenever $N\geqslant c_{2}(\log d+\log(1/\delta))$. Taking the sum over all
$i\in[d]$ concludes the proof. ∎

##  7 Concluding remarks

Several natural questions follow. The first is on the existence of
computationally efficient estimators achieving our bounds. It is known that
getting a polynomial time algorithm with a linear dependence on $\varepsilon$
in the strong contamination model matching the bound of Theorem 1 is a
challenging problem, even when the covariance matrix is identity. Covariance
estimation is an even harder problem from the computational perspective. To
the best of our knowledge, it is unknown if there is a polynomial time
algorithm achieving the statistical performance of Theorem 2 even with the
much weaker $\sqrt{\varepsilon}$-dependence on the contamination level.

One simpler question is if our bounds can be generalized beyond the Gaussian
case. The answer is yes, and we opted for explicit Gaussian computations only
to make our proofs more reader-friendly. In particular, the proof of Theorem 1
only uses the following properties of the distribution:

  1. 1.

The distribution of $X-\mu$ is symmetric around the origin.

  2. 2.

The distribution is spherically symmetric. That is, for any $v\in S^{d-1}$,
the distribution of $\langle X-\mu,v\rangle/\sqrt{v^{\mathsf{T}}\Sigma v}$
does not depend on $v$. Denote the density function of this distribution by
$f$.

  3. 3.

The inverse of the cumulative distribution function $F$ that corresponds to
the density $f$ satisfies $|F^{-1}(1/2\pm\varepsilon)-F^{-1}(1/2)|\leqslant
c\varepsilon$ for some $c>0$ and small enough $\varepsilon$.

  4. 4.

The density function $f$ is separated from zero by some absolute constant for
all $x\in[F^{-1}(1/2-\varepsilon),F^{-1}(1/2+\varepsilon)]$.

  5. 5.

The distribution corresponding to the density function $f$ is sub-Gaussian.
That is, for $Y$ distributed according to this distribution we have
$\|Y\|_{\psi_{2}}\leqslant c$ for some $c>0$.

Following the lines of our proof almost verbatim, one can analyze these more
general distributions. It will be interesting to understand if the sub-
Gaussian tails assumption (Property 5) can be avoided. In fact, assuming
Properties $1$-$4$, and additionally that $\Sigma=I_{d}$, combining our
techniques and the analysis in [DK19, Proposition 1.3], one can build an
estimator $\widehat{\mu}$ satisfying, with probability at least $1-\delta$,

$\|\widehat{\mu}-\mu\|\leqslant
c\left(\sqrt{\frac{d+\log(1/\delta)}{N}}+\varepsilon\right),$

whenever $N\geqslant c_{1}(d+\log(1/\delta))$. Here $c,c_{1}>0$ are some
absolute constants. A similar bound without the sub-Gaussian assumption is
also given by Chen, Gao, and Ren [CGR18, Section 4]. In our case, the sub-
Gaussian assumption (Property 5) is needed to control the moment generating
function when applying Lemma 4, while the proof in [DK19, Proposition 1.3] is
based on the union bound over the $\varepsilon$-net for which we do not need
sub-Gaussian tails in the large deviation regime.

Finally, some of the parameters of Theorem 2 are rather hard to estimate
without making additional assumptions on the sample size, and confidence
level. One can possibly adapt other approaches, such as, for example,
Lepskii’s method [Lep90]. This could provide an alternative way of tuning
these parameters.

#### Acknowledgments.

The authors would like to thank Ankit Pensia for a discussion on differences
between contamination models, Vladimir Ulyanov for a discussion on asymptotic
laws for sample quantiles, and Arnak Dalalyan for many insightful discussions
and useful comments. The work of AM was supported by the grant Investissements
d’Avenir (ANR-11-IDEX-0003/Labex Ecodec/ANR-11-LABX-0047) and by the FAST
Advance grant.

## References

  * [ABM19] Jason Altschuler, Victor-Emmanuel Brunel, and Alan Malek.  Best arm identification for contaminated bandits.  Journal of Machine Learning Research, 20(91):1–39, 2019.
  * [AC11] Jean-Yves Audibert and Olivier Catoni.  Robust linear least squares regression.  The Annals of Statistics, 39(5):2766–2794, 2011.
  * [AZ22] Pedro Abdalla and Nikita Zhivotovskiy.  Covariance estimation: Optimal dimension-free guarantees for adversarial corruption and heavy tails.  arXiv preprint arXiv:2205.08494, 2022.
  * [Bah66] Raj Bahadur.  A note on quantiles in large samples.  The Annals of Mathematical Statistics, 37(3):577–580, 1966.
  * [BL19] Sergey Bobkov and Michel Ledoux.  One-dimensional empirical measures, order statistics, and Kantorovich transport distances, volume 1259 of Memoirs of the American Mathematical Society.  American Mathematical Society, 2019.
  * [BLM13] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  Concentration Inequalities: A Nonasymptotic Theory of Independence.  Oxford University Press, 2013.
  * [BMD22] Amir-Hossein Bateni, Arshak Minasyan, and Arnak S. Dalalyan.  Nearly minimax robust estimator of the mean vector by iterative spectral dimension reduction.  arXiv preprint arXiv:2204.02323, 2022.
  * [Bor75] Christer Borell.  The Brunn-Minkowski inequality in Gauss space.  Inventiones Mathematicae, 30(2):207–216, 1975.
  * [BT12] Stéphane Boucheron and Maud Thomas.  Concentration inequalities for order statistics.  Electronic Communications in Probability, 17:1–12, 2012.
  * [Bur97] Marat Burnashev.  Asymptotic expansions for median estimate of a parameter.  Theory of Probability & Its Applications, 41(4):632–645, 1997\.
  * [Cat16] Olivier Catoni.  PAC-Bayesian bounds for the Gram matrix and least squares regression with a random design.  arXiv preprint arXiv:1603.05229, 2016.
  * [CDG19] Yu Cheng, Ilias Diakonikolas, and Rong Ge.  High-dimensional robust mean estimation in nearly-linear time.  In Proceedings of the thirtieth annual ACM-SIAM symposium on discrete algorithms, pages 2755–2771. SIAM, 2019.
  * [CG17] Olivier Catoni and Ilaria Giulini.  Dimension-free PAC-Bayesian bounds for matrices, vectors, and linear least squares regression.  arXiv preprint arXiv:1712.02747, 2017.
  * [CGR18] Mengjie Chen, Chao Gao, and Zhao Ren.  Robust covariance and scatter matrix estimation under huber’s contamination model.  The Annals of Statistics, 46(5):1932–1960, 2018.
  * [DG92] David L Donoho and Miriam Gasko.  Breakdown properties of location estimates based on halfspace depth and projected outlyingness.  The Annals of Statistics, pages 1803–1827, 1992.
  * [DHF07] Laurens De Haan and Ana Ferreira.  Extreme Value Theory: An Introduction.  Springer Science & Business Media, 2007.
  * [DK19] Ilias Diakonikolas and Daniel M Kane.  Recent advances in algorithmic high-dimensional robust statistics.  arXiv preprint arXiv:1911.05911, 2019.
  * [DKK+17] Ilias Diakonikolas, Gautam Kamath, Daniel M. Kane, Jerry Li, Ankur Moitra, and Alistair Stewart.  Being robust (in high dimensions) can be practical.  In Proceedings of the 34th International Conference on Machine Learning, ICML 2017, volume 70, pages 999–1008, 2017.
  * [DKLP22] Ilias Diakonikolas, Daniel M Kane, Jasper CH Lee, and Ankit Pensia.  Outlier-robust sparse mean estimation for heavy-tailed distributions.  arXiv preprint arXiv:2211.16333, 2022.
  * [DKP20] Ilias Diakonikolas, Daniel M Kane, and Ankit Pensia.  Outlier robust mean estimation with subgaussian rates via stability.  Advances in Neural Information Processing Systems, 33:1830–1840, 2020.
  * [DKS17] Ilias Diakonikolas, Daniel M Kane, and Alistair Stewart.  Statistical query lower bounds for robust estimation of high-dimensional Gaussians and Gaussian mixtures.  In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 73–84, 2017.
  * [DL22] Jules Depersin and Guillaume Lecué.  On the robustness to adversarial corruption and to heavy-tailed data of the Stahel–Donoho median of means.  Information and Inference: A Journal of the IMA, 12 2022.
  * [DM22] Arnak S Dalalyan and Arshak Minasyan.  All-in-one robust estimator of the Gaussian mean.  The Annals of Statistics, 50(2):1193–1219, 2022.
  * [DN04] Herbert A David and Haikady N Nagaraja.  Order Statistics.  John Wiley & Sons, 2004.
  * [DV75] Monroe D. Donsker and Srinivasa Varadhan.  Asymptotic evaluation of certain markov process expectations for large time, I.  Communications on Pure and Applied Mathematics, 28(1):1–47, 1975\.
  * [Giu18] Ilaria Giulini.  Robust dimension-free Gram operator estimates.  Bernoulli, 24(4B):3864–3923, 2018.
  * [Han22] Qiyang Han.  Exact spectral norm error of sample covariance.  arXiv preprint arXiv:2207.13594, 2022.
  * [HL19] Samuel B Hopkins and Jerry Li.  How hard is robust mean estimation?  In Conference on Learning Theory, pages 1649–1682. PMLR, 2019.
  * [HRRS80] Frank R Hampel, Peter J Rousseeuw, Elvezio M Ronchetti, and Werner A Stahel.  Robust statistics: the approach based on influence functions.  Wiley, 1980.
  * [Hub64] Peter J Huber.  Robust estimation of a location parameter.  The Annals of Mathematical Statistics, 35(1):73–101, 1964.
  * [Hub81] Peter J Huber.  Robust statistics.  Wiley Series in Probability and Mathematical Statistics, 1981.
  * [Kie67] Jack Kiefer.  On Bahadur’s representation of sample quantiles.  The Annals of Mathematical Statistics, 38(5):1323–1342, 1967.
  * [KL17] Vladimir Koltchinskii and Karim Lounici.  Concentration inequalities and moment bounds for sample covariance operators.  Bernoulli, 23(1):110–133, 2017.
  * [KMR+19] Yuan Ke, Stanislav Minsker, Zhao Ren, Qiang Sun, and Wen-Xin Zhou.  User-friendly covariance estimation for heavy-tailed distributions.  Statistical Science, 34(3):454–471, 2019.
  * [Kol31] Andrey Kolmogorov.  La méthode de la mediane dans la théorie des erreurs.  Rec. Math. Moscou, 38(3-4):47–50, 1931.
  * [KZ20] Yegor Klochkov and Nikita Zhivotovskiy.  Uniform Hanson-Wright type concentration inequalities for unbounded entries via the entropy method.  Electronic Journal of Probability, 25, 2020.
  * [Lep90] Oleg V Lepskii.  A problem of adaptive estimation in gaussian white noise.  Teoriya Veroyatnostei i ee Primeneniya, 35(3):459–470, 1990.
  * [LM00] Beatrice Laurent and Pascal Massart.  Adaptive estimation of a quadratic functional by model selection.  Annals of Statistics, pages 1302–1338, 2000.
  * [LM19a] Gábor Lugosi and Shahar Mendelson.  Mean estimation and regression under heavy-tailed distributions: A survey.  Foundations of Computational Mathematics, 19(5):1145–1190, 2019\.
  * [LM19b] Gábor Lugosi and Shahar Mendelson.  Near-optimal mean estimators with respect to general norms.  Probability theory and related fields, 175(3):957–973, 2019.
  * [LM21] Gábor Lugosi and Shahar Mendelson.  Robust multivariate mean estimation: The optimality of trimmed mean.  The Annals of Statistics, 49(1):393–410, 2021.
  * [LRV16] Kevin A. Lai, Anup B. Rao, and Santosh S. Vempala.  Agnostic estimation of mean and covariance.  In IEEE 57th Annual Symposium on Foundations of Computer Science, FOCS 2016, pages 665–674, 2016.
  * [MW22] Stanislav Minsker and Lang Wang.  Robust estimation of covariance matrices: Adversarial contamination and beyond.  arXiv preprint arXiv:2203.02880, 2022.
  * [MZ20] Shahar Mendelson and Nikita Zhivotovskiy.  Robust covariance estimation under ${L}_{4}-{L}_{2}$ norm equivalence.  The Annals of Statistics, 48(3):1648–1664, 2020.
  * [Oli16] Roberto I Oliveira.  The lower tail of random quadratic forms with applications to ordinary least squares.  Probability Theory and Related Fields, 166(3-4):1175–1194, 2016\.
  * [OR22] Roberto I Oliveira and Zoraida F Rico.  Improved covariance estimation: Optimal robustness and sub-Gaussian guarantees under heavy tails.  arXiv preprint arXiv:2209.13485, 2022.
  * [RL05] Peter J Rousseeuw and Annick M Leroy.  Robust Regression and Outlier Detection, volume 589.  John Wiley & Sons, 2005.
  * [Sha03] Jun Shao.  Mathematical Statistics.  Springer Science & Business Media, 2003.
  * [Tal14] Michel Talagrand.  Upper and Lower Bounds for Stochastic Processes: Modern Methods and Classical Problems, volume 60.  Springer Science & Business Media, 2014.
  * [TIS76] Boris S Tsirelson, Ildar A Ibragimov, and Vladimir N Sudakov.  Norms of gaussian sample functions.  In Proceedings of the Third Japan—USSR Symposium on Probability Theory, pages 20–41. Springer, 1976.
  * [Ver16] Roman Vershynin.  High-Dimensional Probability: An Introduction with Applications, volume 47 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press, 2016.
  * [Xia19] Dong Xia.  Non-asymptotic bounds for percentiles of independent non-identical random variables.  Statistics & Probability Letters, 152:111–120, 2019.
  * [Zhi21] Nikita Zhivotovskiy.  Dimension-free bounds for sums of independent matrices and simple tensors via the variational principle.  arXiv preprint arXiv:2108.08198, 2021.