# Affine Invariant Covariance Estimation for Heavy-Tailed Distributions

Dmitrii M. Ostrovskii, Alessandro Rudi

► To cite this version:

Dmitrii M. Ostrovskii, Alessandro Rudi. Affine Invariant Covariance Estimation for Heavy-Tailed Distributions. COLT 2019 - 32nd Annual Conference on Learning Theory, Jun 2019, Phoenix, United States. hal-02011464v2

HAL Id: hal-02011464

https://hal.science/hal-02011464v2

Submitted on 23 Sep 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Affine Invariant Covariance Estimation
for Heavy-Tailed Distributions

Dmitrii M. Ostrovskii *†

dostrovs@usc.edu

Alessandro Rudi ‡

alessandro.rudi@inria.fr

September 23, 2019

## Abstract

In this work we provide an estimator for the covariance matrix of a heavy-tailed multivariate distribution. We prove that the proposed estimator $\hat{S}$ admits an affine-invariant bound of the form

$$ (1 - \varepsilon)\mathbf{S} \preceq \hat{\mathbf{S}} \preceq (1 + \varepsilon)\mathbf{S} $$

in high probability, where $\mathbf{S}$ is the unknown covariance matrix, and $\preceq$ is the positive semidefinite order on symmetric matrices. The result only requires the existence of fourth-order moments, and allows for $\varepsilon = O(\sqrt{\kappa^4 d \log(d/\delta)/n})$ where $\kappa^4$ is a measure of kurtosis of the distribution, $d$ is the dimensionality of the space, $n$ is the sample size, and $1 - \delta$ is the desired confidence level. More generally, we can allow for regularization with level $\lambda$, then $d$ gets replaced with the degrees of freedom number. Denoting $\text{cond}(\mathbf{S})$ the condition number of $\mathbf{S}$, the computational cost of the novel estimator is $O(d^2n + d^3 \log(\text{cond}(\mathbf{S})))$, which is comparable to the cost of the sample covariance estimator in the statistically interesting regime $n \ge d$. We consider applications of our estimator to eigenvalue estimation with relative error, and to ridge regression with heavy-tailed random design.

## 1 Introduction

We are interested in estimating the covariance matrix $\mathbf{S} = \mathbb{E}[X \otimes X]$ of a zero-mean random vector $X \in \mathbb{R}^d$ from $n$ independent and identically distributed (i.i.d.) copies $X_1, ..., X_n$ of $X$. This task is crucial – and often arises as a subroutine – in some widely used statistical procedures, such as linear regression, principal component analysis, factor analysis, generalized methods of moments, and mean-variance portfolio selection, to name a few [FHT01, Jol02, Han82, Mar52]. In some of them, the control of $\|\hat{\mathbf{S}} - \mathbf{S}\|_*$, where $\hat{\mathbf{S}}$ is a covariance estimator and $\| \cdot \|_*$ is the spectral, Frobenius or trace norm, does not result in sharp theoretical guarantees. Instead, it might be necessary to estimate the eigenvalues of $\mathbf{S}$ in relative scale, ensuring that

$$ |\lambda_j(\hat{\mathbf{S}}) - \lambda_j(\mathbf{S})| \le \varepsilon\lambda_j(\mathbf{S}), \quad j \in \{1, ..., d\}, $$

holds for $\varepsilon > 0$ (this task arises in the analysis of the subspace iteration method, see [HMT11] and Sec. 6). More generally, one may seek to provide affine-invariant bounds of the form

$$ (1 - \varepsilon)\mathbf{S} \preceq \hat{\mathbf{S}} \preceq (1 + \varepsilon)\mathbf{S}, \quad (1) $$

<footnote>*Viterbi School of Engineering, University of Southern California, Los Angeles, USA.</footnote>

<footnote>†This work has been done while the first author was at the SIERRA Project-Team, of Inria, Paris, France.</footnote>

<footnote>‡SIERRA Project-Team, Inria and École Normale Supérieure, PSL Research University, Paris, France.</footnote>
as in the analysis of linear regression with random design [see HKZ12, and Sec. 2.3 for more details], where ≲ is the positive semidefinite partial order for symmetric matrices. In fact, the basic and very natural *sample covariance estimator*

$$ \tilde{\mathbf{S}} = \frac{1}{n} \sum_{i=1}^{n} X_i \otimes X_i $$

can be shown to satisfy Eq. (1) with probability at least $1 - \delta$, $\delta \in (0, 1]$, and accuracy $\varepsilon$ scaling as $O((\mathrm{rank}(\mathbf{S}) \log(d/\delta)/n)^{1/2})$, provided that $X$ is subgaussian (see Sec. 2.3 for a detailed discussion). However, the assumption of sub-gaussianity might be too strong in the above applications. Going beyond it and similar assumptions is particularly important in mathematical finance, where it is widely accepted that the prices of assets might have heavy-tailed distributions [KJ14, BT03].

We propose a simple variation of the sample covariance estimator (see Algorithm 1) in the form

$$ \hat{\mathbf{S}} = \frac{1}{n} \sum_{i=1}^{n} \alpha_i X_i \otimes X_i, $$

where the coefficients $\alpha_1, \dots, \alpha_n > 0$ are chosen in a data-driven manner. Our main result, stated informally below, shows that the proposed estimator enjoys *high-probability* bounds analogous to those for the sample covariance estimator, under a weak moment assumption on the distribution. Namely, we assume that for some $\kappa \ge 1$ it holds

$$ \mathbb{E}^{1/4}[\langle X, u \rangle^4] \le \kappa \mathbb{E}^{1/2}[\langle X, u \rangle^2], \quad \forall u \in \mathbb{R}^d. \qquad (\mathrm{HT}) $$

In other words, the *kurtosis* of $X$ is bounded by $\kappa$ in all directions.<fnref n="1" /> Kurtosis is an affine-invariant and unitless quantity, and it is uniformly bounded from above by a constant for many common families of multivariate distributions: for example, $\kappa = \sqrt[4]{3}$ for any Gaussian distribution, and $\kappa \le \sqrt[4]{9}$ for the multivariate Student-t distribution with at least 5 degrees of freedom.

Now we are ready to informally state our main result.

**Theorem 1.1** (Simplified version of Thm. 4.1). Under (HT), there exists an estimator $\hat{\mathbf{S}}$ that has computational cost $O(d^2n + d^3)$, and with probability at least $1 - \delta$ satisfies (1) with accuracy

$$ \varepsilon \le 48\kappa^2 \sqrt{\frac{\mathrm{rank}(\mathbf{S}) \log(4d/\delta)}{n}}. \qquad (2) $$

This result shows that the proposed estimator is a valid alternative to the sample covariance estimator: it has comparable accuracy and the same computational complexity, while requiring only boundedness of the fourth moment of $X$ instead of sub-gaussianity. More generally, by allowing a *regularization level* $\lambda > 0$, i.e., using $\hat{\mathbf{S}}_\lambda := \hat{\mathbf{S}} + \lambda \mathbf{I}$ instead of $\hat{\mathbf{S}}$ to estimate $\mathbf{S}_\lambda := \mathbf{S} + \lambda \mathbf{I}$ instead of $\mathbf{S}$ (as required in ridge regression [HKZ12], we can replace $\mathrm{rank}(\mathbf{S})$ with the *degrees of freedom* number

$$ \mathrm{df}_{\lambda}(\mathbf{S}) := \mathrm{Tr}(\mathbf{S}\mathbf{S}_{\lambda}^{-1}). \qquad (3) $$

This leads to a better bound, since $\mathrm{df}_{\lambda}(\mathbf{S})$ is never larger than $\min\{\mathrm{rank}(\mathbf{S}), \mathrm{Tr}(\mathbf{S})/\lambda\}$, and can be way smaller depending on the eigenvalue decay of $\mathbf{S}$: for example, if $\lambda_j(\mathbf{S}) \le j^{-b}$ with $b \ge 1$, then $\mathrm{df}_{\lambda}(\mathbf{S}) \le \lambda^{-1/b}$.

<footnote n="1">¹Note that we use a slightly non-standard definition of kurtosis, extracting the corresponding roots from the moments.</footnote>
**Paper Organization.** In Sec. 2 we recall the known results for the sample covariance matrix estimator under light-tailed assumptions, together with some recent high-probability results for an alternative estimator applicable to heavy-tailed distributions. The novel estimator is presented in Sec. 3 and analyzed and discussed in detail in Sec. 4. In order to achieve the best statistical performance, it requires the knowledge of the distribution parameters $κ$ and $df_λ(S)$ in advance; in Sec. 5 we extend the algorithm, via a variant of Lepskii's method [Lep91], to be adaptive to these quantities. Applications to eigenvalue estimation and ridge regression are discussed in Sec. 6.

**Notation and Conventions.** For $X \in \mathbb{R}^d$, $X \otimes X$ denotes the outer product $XX^T$. W.l.o.g. we assume that $S = E[X \otimes X]$ is full-rank (otherwise we can work on its range). To reduce the clutter of parentheses, we convene that powers and multiplication have priority over the expectation, and we denote the $1/p$-th power of expectation, $p \ge 1$, with $E^{1/p}[\cdot]$. We use $\| \cdot \|$ for the spectral norm of a matrix (unless specified otherwise), as well as for the $\ell_2$-norm of a vector. We shortand $\min(a, b)$ to $a \wedge b$. We use the $O(\cdot)$ notation in a conventional way, and occasionally replace generic constants with $O(1)$. We use the notation $A_λ := A + λI$, where $A \in \mathbb{R}^{d \times d}$ and $I$ is the identity matrix.

# 2 Background and Related Work

In this section we recall some relevant previous work on covariance estimators for light-tailed and heavy-tailed distributions, and provide more intuition about affine-invariant error bounds. Moreover, we introduce basic concepts that will be used later on in the theoretical analysis.

## 2.1 Relative Error Bounds for the Sample Covariance Estimator

When the estimation error is measured by $\|\tilde{S} - S\|$, where $\| \cdot \|$ is the spectral norm, the problem can be reduced, via the Chernoff bounding technique, to the control of the matrix moment generating function, for which one can apply some deep operator-theoretic results such as the Goldon-Thompson inequality [AW02, Oli10] or Lieb's theorem [Tro12]. Alternatively, one may reduce the task to the control of an underlying empirical process, and exploit advanced tools from empirical process theory such as generic chaining [KL17]. Both families of approaches have been focused on the sample covariance estimator and its direct extensions [Ver12, Tro12, Tro15], requiring stronger assumptions on the distribution of $X$ than (HT). In particular, consider the subgaussian moment growth assumption

$$E^{1/p} [ |\langle X, u \rangle|^p ] \le \bar{\kappa} \sqrt{p} E^{1/p} [ |\langle X, u \rangle|^2 ], \quad \forall u \in \mathbb{R}^d \text{ and } p \ge 2, \qquad (SG)$$

which implies (HT) with $\kappa = 2\bar{\kappa}$, where $\kappa$ is defined in Eq. (HT). Define the effective rank of $S$ by

$$\mathbf{r}(S) := \mathrm{Tr}(S) / \|S\|. \qquad (4)$$

The following result is known.

**Theorem 2.1** (Simplified version of [Lou14, Prop. 3]). Under (SG), the sample covariance estimator $\tilde{S}$ with probability at least $1 - \delta$, $\delta \in (0, 1]$, satisfies

$$\|\tilde{S} - S\| \le O(1) \bar{\kappa}^2 \|S\| \sqrt{\frac{\mathbf{r}(S) \log(2d/\delta)}{n}}, \qquad (5)$$

provided that $n \ge \tilde{O}(1)\mathbf{r}(S)$, where $\tilde{O}(1)$ hides polynomial dependency on $\log(2d/\delta)$ and $\log(2n/\delta)$.
It can be shown [LT13, Prop. 6.10] that Eq. (5) nearly optimally depends on $\mathbf{r}(\mathbf{S})$, $\bar{\kappa}$, and $1/\delta$.<fnref n="2" /> Another remarkable property of this bound is that it is almost dimension-independent: up to a logarithmic factor, the complexity of estimating $\mathbf{S}$ is independent of the ambient dimension $d$. Instead, it is controlled by the distribution-dependent quantity $\mathbf{r}(\mathbf{S})$, which always satisfies $\mathbf{r}(\mathbf{S}) \le \text{rank}(\mathbf{S}) \le d$, and can be much smaller than $\text{rank}(\mathbf{S})$ when the distribution of $X$ lies close to a low-dimensional linear subspace, i.e., when $\mathbf{S}$ has only a few relatively large eigenvalues.

## 2.2 Relative Error Bounds for Heavy-Tailed Distributions

It is possible to obtain relative error bounds of the form $\|\hat{\mathbf{S}} - \mathbf{S}\|$, including the ones in high probability, under weak moment assumptions such as (HT), considering other estimators than the sample covariance matrix. In particular, [WM17] propose an estimator based on the idea of clipping observations with large norm. Formally, they define the truncation map $\psi_{\theta}: \mathbb{R} \to \mathbb{R}$,

$$ \psi_{\theta}(x) := (|x| \wedge \theta) \operatorname{sign}(x), \qquad (6) $$

given a certain threshold $\theta > 0$, and consider the estimator

$$ \hat{\mathbf{S}}^{\text{WM}} := \frac{1}{n} \sum_{i=1}^{n} \rho_{\theta}(\|X_i\|) X_i \otimes X_i, \quad \text{where } \rho_{\theta}(x) := \psi_{\theta}(x^2)/x^2. \qquad (7) $$

In other words, one simply truncates observations with squared norm larger than $\theta$ prior to averaging. This estimator is a key ingredient in our Algorithm 1, and we now summarize its statistical properties.

**Theorem 2.2** ([WM17, Lem. 2.1 and Lem. 5.7]). Define the matrix second moment statistic $W := \|\mathbb{E}[\|X\|^2 X \otimes X]\|$. Let $\bar{W} \ge W$, and $\delta \in (0, 1]$. Then estimator $\hat{\mathbf{S}}^{\text{WM}}$, cf. (7), with $\theta = \sqrt{n\bar{W}/\log(2d/\delta)}$ with probability at least $1-\delta$ satisfies $\|\hat{\mathbf{S}}^{\text{WM}} - \mathbf{S}\| \le 2\sqrt{\bar{W}\log(2d/\delta)/n}$.

In contrast with Thm. 2.1, Thm. 2.2 claims subgaussian concentration for the spectral-norm loss under the weak moment assumption (HT). Moreover, we arrive at the relative error bound akin to (5):

$$ \|\hat{\mathbf{S}}^{\text{WM}} - \mathbf{S}\| \le 2\kappa^2 \|\mathbf{S}\| \sqrt{\frac{\mathbf{r}(\mathbf{S}) \log(2d/\delta)}{n}}, \qquad (8) $$

if we bound the second moment statistic as

$$ W \le \|\mathbb{E}\|X\|^2 X \otimes X\| \le \kappa^4 \|\mathbf{S}\|^2 \mathbf{r}(\mathbf{S}), \qquad (9) $$

see [WM17, Lem. 2.3 and Cor. 5.1], and choose the appropriate truncation level

$$ \theta = \kappa^2 \|\mathbf{S}\| \sqrt{n\mathbf{r}(\mathbf{S})/\log(2d/\delta)} \ge \sqrt{nW/\log(2d/\delta)}. $$

Since this choice depends on the unknown $W$, one can use a larger value, which will result in the inflation of the right-hand side of Eq. (8). An alternative is to adapt to the unknown $W$ via Lepskii's method [Lep91] as described in [WM17, Thm 2.1]. To conclude, the estimator $\hat{\mathbf{S}}^{\text{WM}}$ enjoys subgaussian relative error bounds under the fourth moment assumption (HT), while having essentially the same computation cost as the sample covariance estimator.

<footnote n="2">²In fact, [KL17, Theorem 9] replaces $\mathbf{r}(\mathbf{S}) \log(2d/\delta)$ by $\mathbf{r}(\mathbf{S}) + \log(1/\delta)$, making the bound dimension-independent, but does not specify the dependency on $\bar{\kappa}$. Similar results have been obtained under (HT) for robust estimators (see Sec. 2.2), which, however, are computationally intractable [MZ18].</footnote>
## 2.3 Affine-Invariant Bounds for the Sample Covariance Estimator

As we have seen previously, estimator $\hat{\mathbf{S}}^{\text{WM}}$ has favorable statistical properties compared to $\tilde{\mathbf{S}}$ when the goal is to estimate $\mathbf{S}$ in relative spectral-norm error as in Eq. (8). However, one can instead be interested in providing affine-invariant bounds in the form of Eq. (1). More generally, one may wish to estimate $\mathbf{S}$ only for the eigenvalues greater than some level $\lambda > 0$, that is, to guarantee that

$$ (1 - \varepsilon)\mathbf{S}_{\lambda} \preceq \tilde{\mathbf{S}}_{\lambda} \preceq (1 + \varepsilon)\mathbf{S}_{\lambda}. \qquad (10) $$

The need for such bounds arises, in particular, in random-design ridge regression, where the information about inferior eigenvalues is irrelevant, since it is anyway erased by regularization. Note that Eq. (10), for any $\lambda > 0$, can be reformulated in terms of the $\mathbf{S}_{\lambda}^{-1/2}$-transformed spectral norm:

$$ \|\mathbf{S}_{\lambda}^{-1/2}(\hat{\mathbf{S}} - \mathbf{S})\mathbf{S}_{\lambda}^{-1/2}\| \le \varepsilon. \qquad (11) $$

The task of obtaining such bounds, with arbitrary regularization level $0 \le \lambda \le \|\mathbf{S}\|$, will be referred to as calibrated covariance estimation. Generally, this task is harder than proving relative-error bounds in the spectral norm such as Eq. (8): the latter is equivalent, up to a constant factor loss of accuracy, to proving Eq. (11) with $\lambda = \|\mathbf{S}\|$. On the other hand, calibrated covariance estimation also subsumes Eq. (1) by taking $\lambda = O(\lambda_{\min}(\mathbf{S}))$, where $\lambda_{\min}(\mathbf{S})$ is the smallest eigenvalue of $\mathbf{S}$.

Now, one can make a simple observation that for the sample covariance estimator $\tilde{\mathbf{S}}$, calibrated bounds of the form (11) “automatically” follow from the dimension-free spectral-norm bounds akin to (5) or (8) due to its affine equivariance. Indeed, $\mathbf{J} := \mathbf{S}_{\lambda}^{-1/2}\mathbf{S}\mathbf{S}_{\lambda}^{-1/2}$ is precisely the covariance matrix of the “$\lambda$-decorrelated” observations $Z_i = \mathbf{S}_{\lambda}^{-1/2}X_i$, for which the sample covariance estimator is given by $\tilde{\mathbf{J}} := \frac{1}{n}\sum_{i=1}^n Z_i \otimes Z_i = \mathbf{S}_{\lambda}^{-1/2}\tilde{\mathbf{S}}\mathbf{S}_{\lambda}^{-1/2}$. Hence, we can apply the spectral-norm bound (5), replacing $\mathbf{S}$ and $\tilde{\mathbf{S}}$ with $\mathbf{J}$ and $\tilde{\mathbf{J}}$. Using the fact that $\|\mathbf{J}\| \le 1$ for any $\lambda \ge 0$, and that assumptions (HT), (SG) are themselves invariant under (non-singular) linear transforms, we obtain

$$ \mathbb{E}^{1/2} \left[ \| \mathbf{S}_{\lambda}^{-1/2} (\tilde{\mathbf{S}} - \mathbf{S}) \mathbf{S}_{\lambda}^{-1/2} \right]^2 \le O(1) \bar{\kappa}^2 \sqrt{\frac{\mathrm{df}_{\lambda}(\mathbf{S}) \log(2d)}{n}} \qquad (12) $$

once $n \ge \tilde{O}(1)\mathrm{df}_{\lambda}(\mathbf{S})$, where $\mathrm{df}_{\lambda}(\cdot)$ is defined in Eq. (3) and ranges from $O(\mathrm{tr}(\mathbf{S}))$ to $\mathrm{rank}(\mathbf{S}) \le d$ as $\lambda$ decreases from $\|\mathbf{S}\|$ to zero. In fact, when $\bar{\kappa}$ is a constant, this rate is known to be asymptotically minimax-optimal over certain natural classes of covariance matrices, e.g., Toeplitz matrices with spectra discretizing those of Hölder-smooth functions [BL08, CZZ10]. It is thus reasonable to ask whether one can extend Eq. (12) in the same manner as Eq. (8) extends Eq. (5). In other words, can one provide a high-probability guarantee for the calibrated error (cf. Eq. (11)) of the estimator $\hat{\mathbf{S}}^{\text{WM}}$ (cf. Eq. (7)) under assumption (HT)? The immediate difficulty is that $\hat{\mathbf{S}}^{\text{WM}}$ – in fact, the only estimator for which finite-sample high-probability guarantees under fourth moment assumptions are known to us – does not allow for the same reasoning as $\tilde{\mathbf{S}}$ because of the non-linearity introduced by the truncation map. On the other hand, the desired bounds are achieved by the “oracle” estimator that truncates the “$\lambda$-decorrelated” vectors $Z_i = \mathbf{S}_{\lambda}^{-1/2}X_i$ with accordingly adjusted $\theta$:

$$ \hat{\mathbf{S}}^{\circ} := \frac{1}{n} \sum_{i=1}^{n} \rho_{\theta}(\|\mathbf{S}_{\lambda}^{-1/2} X_i\|) X_i \otimes X_i = \mathbf{S}_{\lambda}^{1/2} \left[ \frac{1}{n} \sum_{i=1}^{n} \rho_{\theta}(\|Z_i\|) Z_i \otimes Z_i \right] \mathbf{S}_{\lambda}^{1/2}. \qquad (13) $$

cf. Eq. (7). Unfortunately, this estimator is unavailable since $Z_i$'s are not observable. In what follows, we present our main methodological contribution: an estimator that achieves the stated
goal, and moreover, has a similar complexity of computation and storage as the sample covariance matrix.

**Remark 2.1.** Some robust covariance estimators, such as MCD or MVE [Cam80, LR91, RD99], are affine equivariant, but to the best of our knowledge, the desired bounds are not known for them. On the other hand, [Oli16] shows that if one only seeks for the left-hand side bound in (10), the sample covariance estimator suffices under (HT).

## 3 Proposed Estimator

Our goal can be summarized as follows: given $\delta \in (0, 1]$, and $\lambda \ge 0$, provide an estimate $\hat{\mathbf{S}}$ satisfying

$$ \|\mathbf{S}_{\lambda}^{-1/2}(\hat{\mathbf{S}} - \mathbf{S})\mathbf{S}_{\lambda}^{-1/2}\| \le O(1)\kappa^2 \sqrt{\frac{\mathrm{df}_{\lambda}(\mathbf{S}) \log(2d/\delta)}{n}}, $$

where $\kappa$ is the kurtosis parameter of $X$ (cf. (HT)). Moreover, we can restrict ourselves to the case $\lambda \le \|\mathbf{S}\|$, since otherwise the task is resolved by the estimator $\hat{\mathbf{S}}^{\mathrm{WM}}$ as can be seen from Eq. (8).

As we have seen before, the oracle estimator $\hat{\mathbf{S}}^o$ introduced in the previous section (cf. Eq. (13)) achieves the stated goal, but is unavailable since it depends explicitly on $\mathbf{S}$. The key idea of our construction is to approximate $\hat{\mathbf{S}}^o$ in an iterative fashion – roughly, to start with $\hat{\mathbf{S}}^{(0)} = \hat{\mathbf{S}}^{\mathrm{WM}}$, which is already a good estimate for $\mathbf{S}_{\lambda}$ with the crudest regularization level $\lambda = \|\mathbf{S}\|$ due to Eq. (8), and then iteratively refine the estimate by computing

$$ \hat{\mathbf{S}}^{(t+1)} = \frac{1}{n} \sum_{i=1}^{n} \rho_{\theta}(\|[\hat{\mathbf{S}}_{\lambda}^{(t)}]^{-1/2} X_{i}\|) X_{i} \otimes X_{i}. \quad (14) $$

To make this simple idea work, we need to adjust it in two ways. Firstly, $\hat{\mathbf{S}}^{(t)}$ depends on the observations $X_i$, hence the random vectors $[\hat{\mathbf{S}}_{\lambda}^{(t)}]^{-1/2} X_i$ are not independent. To simplify the analysis, we split the sample into batches corresponding to different iterations, and at each iteration use observations of the new batch instead of $X_i$'s in Eq. (14). Since $\hat{\mathbf{S}}^{(t)}$ is independent from the new observations, we can apply Thm. 2.2 conditionally at each step.

Secondly, as discussed before, the estimator $\hat{\mathbf{S}}^{\mathrm{WM}}$ given by (7) already solves the problem for $\lambda = \|\mathbf{S}\|$. To achieve Eq. (11) for a given $\lambda$, the idea is to start with $\lambda_0 = \|\mathbf{S}\|$, and reduce $\lambda_t$ by a constant factor at each iteration, so that the error $\|\mathbf{S}_{\lambda}^{-1/2}(\hat{\mathbf{S}}^{(t)} - \mathbf{S})\mathbf{S}_{\lambda}^{-1/2}\|$ remains controlled for $\lambda = \lambda_t$ at each step. This way we also ensure that the total number of iterations is logarithmic in $\|\mathbf{S}\|/\lambda$, and, in particular, logarithmic in the condition number $\|\mathbf{S}\|/\lambda_{\min}(\mathbf{S})$ when $\lambda \ge \lambda_{\min}(\mathbf{S})$.

Algorithm 1 presented below implements these ideas. Note that the final batch of observations takes a half of the overall sample: this is needed to achieve the best possible accuracy (up to a constant factor) for the final regularization level, while at the previous levels it suffices to maintain the accuracy $\varepsilon = 1/2$, and one can use smaller batches taking up a half of the sample in total, see Lem. 4.1 in Sec. 4 for details. Once the coefficients $\alpha_i^{(t)}$ at the given step have been computed, the new estimate $\hat{\mathbf{S}}^{(t+1)}$ reduces to the sample covariance matrix of the weighted observations, which can be computed in time $O(d^2m)$ where $m$ is the size of the batch. The total cost of these computations in the course of the algorithm is thus $O(d^2n)$. As for $\alpha_i^{(t)}$, they are obtained by first performing the Cholesky decomposition [GVL12] of $\hat{\mathbf{S}}_{\lambda_t}^{(t)}$, i.e., finding the unique lower-triangular matrix $\mathbf{R}_t$ such that $\mathbf{R}_t\mathbf{R}_t^{\top} = \hat{\mathbf{S}}_{\lambda_t}^{(t)}$ which requires $O(d^3)$ in time and $O(d^2)$ in
space,<fnref n="3" /> and then computing each product $\mathbf{R}_t^{-1}X_i^{(t+1)}$ by solving the corresponding linear system in $O(d^2)$. The total complexity of Algorithm 1 is thus

$O(d^2n + d^3 \log(L/\lambda))$ in time, $O(d^2)$ in space.

Moreover, the time complexity becomes $O(d^2n)$ when $n \gg d \log(L/\lambda)$; as we show next, this is anyway required to obtain a statistical performance guarantee. Note moreover that it is possible to obtain the non-regularized version of Eq. (1) by choosing $\lambda = O(\lambda_{\min}(\mathbf{S}))$, in time $O(d^2n + d^3 \log(\text{cond}(\mathbf{S})))$, where $\text{cond}(\mathbf{S}) = \|\mathbf{S}\|/\lambda_{\min}$ is the condition number of $\mathbf{S}$.

**Remark 3.1.** In practice, sample splitting in Algorithm 1 could be avoided, and iterations could be performed on the same sample. We expect our statistical guarantees in Sec. 4 to extend to this setup.

Next we present a statistical guarantee for Algorithm 1, and suggest a way to select the parameters.

<caption>**Algorithm 1 Robust Calibrated Covariance Estimation**</caption>

**Input:** $X_1, \dots, X_n \in \mathbb{R}^d$, $\delta \in (0, 1]$, regularization level $\lambda \le \|\mathbf{S}\|$, $L \ge \|\mathbf{S}\|$, truncation level $\theta > 0$

1: $\lambda_0 = L$, $T = \lceil \log_2(L/\lambda) \rceil$, $m = \lfloor n/(2(T+1)) \rfloor$
2: $\alpha_i^{(0)} = L\rho_\theta(\|X_i\|/\sqrt{L})$, for $i \in \{1, \dots, m\}$, with $\rho_\theta(x) = \psi_\theta(x^2)/x^2$ and $\psi_\theta(\cdot)$ as in Eq. (6)
3: $\hat{\mathbf{S}}^{(0)} = \frac{1}{m} \sum_{i=1}^m \alpha_i^{(0)} X_i \otimes X_i$
4: **for** $t = 0$ **to** $T-1$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;$(X_1^{(t+1)}: X_m^{(t+1)}) = (X_{m(t+1)+1}: X_{m(t+1)+m})$ &nbsp;&nbsp;&nbsp;&nbsp;\# Obtain a new batch
6: &nbsp;&nbsp;&nbsp;&nbsp;**R**_t = Cholesky($\hat{\mathbf{S}}_{\lambda_t}^{(t)}$)
7: &nbsp;&nbsp;&nbsp;&nbsp;$$\alpha_i^{(t+1)} = \rho_\theta(\|\mathbf{R}_t^{-1}X_i^{(t+1)}\|), \text{ for } i \in \{1, \dots, m\}$$
8: &nbsp;&nbsp;&nbsp;&nbsp;$$\hat{\mathbf{S}}^{(t+1)} = \frac{1}{m} \sum_{i=1}^m \alpha_i^{(t+1)} X_i^{(t+1)} \otimes X_i^{(t+1)}$$
9: &nbsp;&nbsp;&nbsp;&nbsp;$$\lambda_{t+1} = \lambda_t/2$$
10: **end for**
11: $r = n - m(T+1)$ &nbsp;&nbsp;&nbsp;&nbsp;\# Size of the remaining sample, roughly $n/2$
12: $\theta_T = 2\theta(T+1)^{1/2} (1 + \log(T+1)/\log(4d/\delta))^{1/2}$ &nbsp;&nbsp;&nbsp;&nbsp;\# Final truncation level
13: $(X_1^*: X_r^*) = (X_{mT+1}: X_n)$ &nbsp;&nbsp;&nbsp;&nbsp;\# Remaining sample
14: $\mathbf{R}_T = \text{Cholesky}(\hat{\mathbf{S}}_{\lambda_T}^{(T)})$
15: $\alpha_i^* = \rho_{\theta_T}(\|\mathbf{R}_T^{-1}X_i^*\|), \text{ for } i \in \{1, \dots, r\}$
16: $\hat{\mathbf{S}}^* = \frac{1}{r} \sum_{i=1}^r \alpha_i^* X_i^* \otimes X_i^*$ &nbsp;&nbsp;&nbsp;&nbsp;\# Final estimate

**Output:** $\hat{\mathbf{S}}^*$

# 4 Statistical Guarantee

In Theorem 4.1 below, we show that the estimator produced by Algorithm 1 achieves a high-probability bound of the type (11) requiring only the existence of the fourth-order moments of $X$, and the correct choice of the truncation level $\theta$. We begin with the lemma that justifies the proposed update rule.

<footnote n="3">³Choletsky decomposition is known to work whenever the condition number (in our case $\|\mathbf{S}\|/\lambda$) is dominated by the inverse machine precision; when this condition does not hold, one could add some extra tricks such as pivoting, which still results in $O(d^3)$ complexity [GVL12].</footnote>
**Lemma 4.1.** Let $\hat{\mathbf{S}}$ be a symmetric estimate of $\mathbf{S}$ such that, for some $\lambda > 0$,

$$ \|\mathbf{S}_{\lambda}^{-1/2}(\hat{\mathbf{S}} - \mathbf{S})\mathbf{S}_{\lambda}^{-1/2}\| \le 1/2. \qquad (15) $$

Conditioned on $\hat{\mathbf{S}}$, let $X_1, \dots, X_m$ be i.i.d., have zero mean, covariance $\mathbf{S}$, and finite fourth-order moments. Let $\kappa$ be the associated (conditional) kurtosis as in Eq. (HT). Define $\mathbf{S}^{(+)}$ as

$$ \hat{\mathbf{S}}^{(+)} := \frac{1}{m} \sum_{j=1}^{m} \rho_{\theta}(\|\hat{\mathbf{S}}_{\lambda}^{-1/2} X_j\|) X_j \otimes X_j, \qquad (16) $$

with $\rho_{\theta}$ defined in Eqs. (6)-(7). Choose $\theta \ge 2\sqrt{2\kappa^2}\sqrt{m\mathrm{df}_{\lambda}(\mathbf{S})/\log(2d/\delta)}$, where $\mathrm{df}_{\lambda}(\mathbf{S})$ is defined by Eq. (3). Then with conditional probability at least $1 - \delta$ over $(X_1, \dots, X_m)$ it holds

$$ \| \mathbf{S}_{\lambda/2}^{-1/2} (\hat{\mathbf{S}}^{(+)} - \mathbf{S}) \mathbf{S}_{\lambda/2}^{-1/2} \| \le \frac{6\theta \log(2d/\delta)}{m}. \qquad (17) $$

Lemma 4.1 is proved in Appendix B. Its role is to guarantee the stability of the iterative process in Algorithm 1 when we pass to the next regularization level by $\lambda^{(t)} \leftarrow \lambda^{(t-1)/2}$. Indeed, if the size $m$ of the new batch is large enough, the right-hand side of Eq. (17) can be made smaller than $1/2$, which allows to apply Lemma 4.1 sequentially. We are now ready to present the guarantee for Algorithm 1.

**Theorem 4.1.** Let $X_1, \dots, X_n$ be i.i.d. zero-mean random vectors in $\mathbb{R}^d$ satisfying $\mathrm{E}[X_i \otimes X_i] = \mathbf{S}$ and the kurtosis assumption (HT). Let Algorithm 1 be run with $\delta \in (0, 1]$, $0 < \lambda \le \|\mathbf{S}\| \le L$, and

$$ \theta \ge \theta_* := 2\kappa^2 \sqrt{\frac{n\mathrm{df}_{\lambda}(\mathbf{S})}{\mathbf{q}\log(4\mathbf{q}d/\delta)}}, \quad \text{where } \mathbf{q} := \lceil \log_2(L/\lambda) \rceil + 1. \qquad (18) $$

Whenever the sample size satisfies

$$ n \ge 48\mathbf{q}\theta \log(4\mathbf{q}d/\delta), \qquad (19) $$

the resulting estimator $\hat{\mathbf{S}}$ with probability at least $1 - \delta$ satisfies

$$ \| \mathbf{S}_{\lambda}^{-1/2} (\hat{\mathbf{S}} - \mathbf{S}) \mathbf{S}_{\lambda}^{-1/2} \| \le \frac{24\theta \sqrt{\mathbf{q} \log(4\mathbf{q}d/\delta) \log(4d/\delta)}}{n}. \qquad (20) $$

The above theorem shows that when the conditions on $\theta$ and $n$ are met, the proposed estimator satisfies an affine-invariant error bound with accuracy of the same order as the one available for the sample covariance estimator (cf. Eq. (12)) under the more stringent sub-gaussian assumption. This is made explicit in the next corollary, where we simply put $\theta = \theta_*$ as suggested by Eq. (18), obtaining the bound (cf. Eq. (22)) that matches Eq. (12) up to a constant factor and the replacement of $\bar{\kappa}$ with $\kappa$.

**Corollary 4.1.** Under the premise of Theorem 4.1, assume that

$$ n \ge 96^2 \kappa^4 \mathrm{qdf}_{\lambda}(\mathbf{S}) \log(4\mathbf{q}d/\delta), \qquad (21) $$

and choose $\theta = \theta_*$, cf. Eq. (18). Then the estimator given by Algorithm 1 w.p. at least $1 - \delta$ satisfies

$$ \| \mathbf{S}_{\lambda}^{-1/2} (\hat{\mathbf{S}} - \mathbf{S}) \mathbf{S}_{\lambda}^{-1/2} \| \le 48\kappa^2 \sqrt{\frac{\mathrm{df}_{\lambda}(\mathbf{S}) \log(4d/\delta)}{n}}. \qquad (22) $$

We conclude with a remark on choosing $L$, while in Sec. 5 we will provide an adaptive version of the estimator based on a version of Lepskii's method [Lep91] in which $\theta$ is tuned automatically.

**Remark 4.1 (Choosing $L$).** One can simply put $L = 2\|\hat{\mathbf{S}}^{\mathrm{WM}}\|$, requiring that $\lambda \le \frac{2}{3}\|\hat{\mathbf{S}}^{\mathrm{WM}}\|$, with $\hat{\mathbf{S}}^{\mathrm{WM}}$ defined in Eq. (7) and using an independent subsample of size $n/(T+1)$. Indeed, under Eq. (19), the result of Theorem 2.2 ensures that $\frac{2}{3}\|\hat{\mathbf{S}}^{\mathrm{WM}}\| \le \|\mathbf{S}\| \le 2\|\hat{\mathbf{S}}^{\mathrm{WM}}\|$ with probability $\ge 1 - \delta$.
## 4.1 Proof of Theorem 4.1

Note that $q = T + 1$ is the number of batches processed by the end of the for-loop in Algorithm 1. Thus, using that $m = \lfloor n/(2q) \rfloor$ and $n \ge 4q$ (see Eq. (18)-(19) and use that $\kappa \ge 1$, $\text{df}_\lambda \ge 1$), we get

$$ n/(4q) \le m \le n/(2q). \qquad (23) $$

We will proceed by induction over the steps $0 \le t \le T$, showing that

$$ \| \mathbf{S}_{\lambda_t}^{-1/2} (\hat{\mathbf{S}}^{(t)} - \mathbf{S}) \mathbf{S}_{\lambda_t}^{-1/2} \| \le 1/2 \qquad (24) $$

holds for $t = 0, \dots, T$ with probability $\ge (1 - \delta/(2q))^{t+1}$. Then we will derive Eq. (20) as a corollary.

1.°. For the base, we can apply Thm. 2.2, exploiting that $\hat{\mathbf{S}}^{(0)} = L\hat{\mathbf{S}}^{\text{WM}}$ for the (renormalized) initial batch $\frac{1}{\sqrt{L}}(X_1, \dots, X_m)$. Thus, with probability at least $1 - \delta/(2q)$ over this batch, it holds

$$ \frac{1}{L} \| \hat{\mathbf{S}}^{(0)} - \mathbf{S} \| \le \frac{2\theta \log(4qd/\delta)}{m} \le \frac{8q\theta \log(4qd/\delta)}{n}, \qquad (25) $$

provided that (recall the condition $\theta \ge \sqrt{nW/\log(2d/\delta)}$ in Thm. 2.2 and combine it with Eq. (9)):

$$ \theta \ge \kappa^2 \|\mathbf{S}\|/L \cdot \sqrt{m\mathbf{r}(\mathbf{S})/\log(4qd/\delta)}. $$

But this follows from Eq. (18), since $\|\mathbf{S}\| \le L$, $m \le n/(2q)$, and $\mathbf{r}(\mathbf{S}) \le 2\text{df}_{\|\mathbf{S}\|}(\mathbf{S}) \le 2\text{df}_\lambda(\mathbf{S})$. Noting that $\|\mathbf{S}_L^{-1}\| \le 1/L$, from Eqs. (19) and (25) we get

$$ \| \mathbf{S}_L^{-1/2} (\hat{\mathbf{S}}^{(0)} - \mathbf{S}) \mathbf{S}_L^{-1/2} \| \le \frac{8q\theta \log(4qd/\delta)}{n} \le \frac{1}{6}. $$

Since $\lambda_0 = L$, the induction base is proved. Note that when $T = 0$, this already results in Eq. (24).

2.°. Let $T \ge 1$ and $0 \le t \le T - 1$. For the induction step, we apply Lemma 4.1 conditionally on the first $t$ iterations, with $\hat{\mathbf{S}}^{(t)}$ in the role of the current estimate, $\lambda_t$ in the role of the current regularization level, and $(X_1^{(t+1)}, \dots, X_m^{(t+1)})$ as the new batch (which is independent from $\hat{\mathbf{S}}^{(t)}$ by construction). By the induction hypothesis, we have Eq. (24) with conditional probability $\ge (1 - \delta/(2q))^{t+1}$. By Lemma 4.1, since $\lambda_t \ge \lambda$ (and thus $\text{df}_{\lambda_t}(\mathbf{S}) \le \text{df}_\lambda(\mathbf{S})$), Eq. (18), when combined with the upper bound in Eq. (23), guarantees that with conditional probability $\ge 1 - \delta/(2q)$ over the new batch,

$$ \| \mathbf{S}_{\lambda_{t+1}}^{-1/2} (\hat{\mathbf{S}}^{(t+1)} - \mathbf{S}) \mathbf{S}_{\lambda_{t+1}}^{-1/2} \| \le 6\theta \log(4qd/\delta)/m \le 24q\theta \log(4qd/\delta)/n \le 1/2. $$

Here in the second transition we used the lower bound of Eq. (23), and in the last transition we used Eq. (19). Thus, the induction claim is proved. In particular, we have obtained that the bound

$$ \| \mathbf{S}_{\lambda_T}^{-1/2} (\hat{\mathbf{S}}^{(T)} - \mathbf{S}) \mathbf{S}_{\lambda_T}^{-1/2} \| \le 1/2 $$

holds with probability at least $(1 - \delta/(2q))^q \ge 1 - \delta/2$ over the first $q = T + 1$ batches.

3.°. Finally, to obtain Eq. (20), we apply Lemma 4.1 once again, this time conditioning on $\hat{\mathbf{S}}^{(T)}$, and using the last batch $(X_{n-r+1}, \dots, X_n)$ with the final estimator $\hat{\mathbf{S}}$ in the role of $\hat{\mathbf{S}}^{(+)}$. Note that the first condition Eq. (15) of Lemma 4.1 follows from the just proved induction claim. On the other hand, by Eq. (18) the final truncation level $\theta_T$, cf. line 12 of Algorithm 1, satisfies

$$ \theta_T = 2\theta \sqrt{q \log(4qd/\delta)/\log(4d/\delta)} \ge 4\kappa^2 \sqrt{n \text{df}_\lambda(\mathbf{S})/\log(4d/\delta)}. $$
The number of degrees of freedom is a stable quantity: we can easily prove (see Lem. A.1 in Appendix) that $\mathrm{df}_{\lambda/2}(\mathbf{S}) \le 2\mathrm{df}_{\lambda}(\mathbf{S})$. On the other hand, $\mathrm{df}_{\lambda_T} \le \mathrm{df}_{\lambda/2}$ since $\lambda_T \ge \lambda/2$. Using that, we have

$$ \theta_T \ge 2\sqrt{2\kappa^2}\sqrt{n\mathrm{df}_{\lambda_T}(\mathbf{S})/\log(4d/\delta)} \ge 2\sqrt{2\kappa^2}\sqrt{r\mathrm{df}_{\lambda_T}(\mathbf{S})/\log(4d/\delta)}, $$

meeting the requirement on the truncation level imposed in Lemma 4.1. Applying the lemma, and using that $r \ge n/2$, we get that with conditional probability $\ge 1 - \delta/2$ over the last batch $(X_{n-r+1}, \dots, X_n)$,

$$ \|\mathbf{S}_{\lambda_T/2}^{-1/2}(\hat{\mathbf{S}} - \mathbf{S})\mathbf{S}_{\lambda_T/2}^{-1/2}\| \le \frac{6\theta_T \log(4d/\delta)}{r} \le \frac{24\theta\sqrt{\mathbf{q}\log(4\mathbf{q}d/\delta)\log(4d/\delta)}}{n}. $$

Since $\lambda_T \le \lambda$ implies $\|\mathbf{S}_{\lambda_T/2}\mathbf{S}_\lambda^{-1}\| \le 1$, by the union bound we arrive at Eq. (20). $\square$

# 5 Adaptive Estimator

One limitation of Algorithm 1 is that the truncation level $\theta$ has to be chosen in advance in order to obtain the optimal statistical performance (see Theorem 4.1), and the optimal choice $\theta_*$ (see Eq. (18)) depends on $\kappa$, $\mathrm{df}_\lambda(\mathbf{S})$ that are usually unknown. To address this, we propose an adaptive estimator (see Algorithm 2), in which Algorithm 1 is combined with a Lepskii-type procedure [Lep91], resulting in a near-optimal guarantee without the knowledge of $\theta_*$. Namely, let us be given a range $0 < \theta_{\min} \le \theta_{\max}$ known to contain $\theta_*$ but possibly very loose, and define the logarithmic grid

$$ \theta_j = 2^j \theta_{\min}, \quad \text{where } j \in \mathcal{J} := \{j \in \mathbb{Z} : \theta_{\min} \le \theta_j \le 2\theta_{\max}\}. \qquad (26) $$

Define also

$$ \varepsilon_j := \frac{24\theta_j \sqrt{\mathbf{q} \cdot \log(4\mathbf{q}d/\mathcal{J}/\delta) \cdot \log(4d/\mathcal{J}/\delta)}}{n}. \qquad (27) $$

Later on we will show (see Thm. 4.1) that $\varepsilon_j$ is the error bound, with probability at least $1 - \delta$, for the estimator produced by Algorithm 1 with truncation level $\theta = \theta_j$. In Algorithm 2, we first compute $\varepsilon_j$ and the basic estimators $\hat{\mathbf{S}}_j := \hat{\mathbf{S}}[\theta_j]$ for all truncation levels $\theta_j$, then select

$$ \hat{\mathcal{J}} = \min \left\{ j \in \mathcal{J} : \forall j' \ge j \text{ s.t. } j' \in \mathcal{J} \text{ it holds } \| \hat{\mathbf{S}}_{j',\lambda}^{-1/2} (\hat{\mathbf{S}}_{j'} - \hat{\mathbf{S}}_j) \hat{\mathbf{S}}_{j',\lambda}^{-1/2} \| \le 2(\varepsilon_{j'} + \varepsilon_j) \right\}, \quad (28) $$

and output $\hat{\mathcal{S}}_{\hat{\mathcal{J}}}$ as the final estimator. In Thm. 5.1 below, we show that this estimator admits essentially the same statistical guarantee (in the sense of Eq. (11)) as the “ideal” estimator which uses $\theta = \theta_*$. $\square$

**Algorithm 2 Robust Calibrated Covariance Estimation with Adaptive Truncation Level**

**Input:** $X_1, \dots, X_n \in \mathbb{R}^d$, $\delta \in (0, 1]$, regularization level $\lambda \le \|\mathbf{S}\|$, $L \ge \|\mathbf{S}\|$, range $[\theta_{\min}, \theta_{\max}]$

1: Form the grid $\mathcal{J} := \{j \in \mathbb{Z} : \theta_{\min} \le \theta_j \le 2\theta_{\max}\}$

2: **for** $j \in \mathcal{J}$ **do**

3: Compute the output $\hat{\mathbf{S}}_j$ of Algorithm 1 with truncation level $\theta_j = 2^j\theta_{\min}$; compute $\varepsilon_j$ by (27)

4: **end for**

**Output:** $\hat{\mathbf{S}}_{\hat{\mathcal{J}}}$ with $\hat{\mathcal{J}} \in \mathcal{J}$ selected according to (28)

Next we present a statistical performance guarantee for Algorithm 2. Its proof, given in Appendix C, hinges upon the observation that the matrix $\hat{\mathbf{S}}_{j',\lambda}$ in the error bound of Eq. (28) can essentially be replaced with its unobservable counterpart $\mathbf{S}_\lambda$; this makes the analyzed errors additive, so that the usual argument for Lepskii’s method could be applied.
**Theorem 5.1.** Assume (HT), and let Algorithm 2 be initialized with $\lambda \le \|\mathbf{S}\|$, $L \ge \|\mathbf{S}\|$, $\delta \in (0, 1]$, and a range $[\theta_{\min}, \theta_{\max}]$ containing the optimal truncation level $\theta_{*}$ given by Eq. (18). Moreover, let

$$n \ge 96q\theta_{\max} \log(4qd|\mathcal{J}|/\delta), \qquad (29)$$

where $q$ is defined in Eq. (18), and $|\mathcal{J}| \le 1 + \log_2(\theta_{\max}/\theta_{\min})$ is the cardinality of the grid defined in Eq. (26). Then the estimator $\hat{\mathbf{S}}_{\hat{\mathcal{J}}}$ produced by Algorithm 2 with probability at least $1 - \delta$ satisfies

$$\|\mathbf{S}_{\lambda}^{-1/2}(\hat{\mathbf{S}}_{\hat{\mathcal{J}}}-\mathbf{S})\mathbf{S}_{\lambda}^{-1/2}\| \le 720\kappa^2 \sqrt{\frac{(1+\rho)\mathrm{df}_{\lambda}(\mathbf{S})\log(4d|\mathcal{J}|/\delta)}{n}}, \quad \text{where } \rho := \frac{\log|\mathcal{J}|}{\log(4qd/\delta)}.$$

From the result of Thm. 5.1, we see that the adaptive estimator nearly attains the best possible stastistical guarantee that corresponds to the optimal value of the truncation level, up to the iterated logarithm of the ratio $\theta_{\max}/\theta_{\min}$. However, the premise of Thm. 5.1 requires $\theta_{\max}$ to be bounded both from above and below (cf. Eqs. (18) and (29)), and the two bounds are compatible only starting from a certain sample size. Next we state a corollary of Thm. 5.1 that explicitly specifies the required sample size (the requirement is similar to Eq. (21)), and provides a reasonables choice of $[\theta_{\min}, \theta_{\max}]$.

**Corollary 5.1.** Assume that we have

$$n \ge 192^2(1 + \rho_J)\kappa^4 \mathrm{qdf}_{\lambda}(\mathbf{S}) \log(4qdJ/\delta), \qquad (30)$$

where $J := 1 + \lceil \frac{1}{2} \log_2(n/96q) \rceil$ and $\rho_J := \log(J)/\log(4qd/\delta)$. Then the premise of Theorem 5.1 holds for the grid $\mathcal{J}$ with cardinality $J$ defined by $\theta_{\max} = n/96q \log(4qdJ/\delta)$, $\theta_{\min} = 2^{1-J}\theta_{\max}$.

# 6 Applications

## 6.1 Relative-Scale Bounds for Eigenvalues

Recall that the bounds obtained in Thms. 4.1 and 5.1 for the estimators $\hat{\mathbf{S}}$ given by Algorithms 1–2 read

$$(1 - \varepsilon)\mathbf{S}_{\lambda} \preceq \hat{\mathbf{S}}_{\lambda} \preceq (1 + \varepsilon)\mathbf{S}_{\lambda} \qquad (31)$$

for certain accuracy $\varepsilon < 1/2$ and regularization level $\lambda \ge 0$, provided that the sample size is large enough. Using that the positive-semidefinite order preserves the order of eigenvalues, we obtain the corollary of Thm. 4.1 for eigenvalue estimation (Thm. 5.1 has a similar corollary, which we omit).

**Corollary 6.1.** Assume that $n$ satisfies Eq. (21), and let $\hat{\mathbf{S}}$ be given by Algorithm 1 with the optimal choice of the truncation level $\theta = \theta_{*}$, cf. Eq. (18). Let also $\|\mathbf{S}\| = \lambda_1 \ge \dots \ge \lambda_d = \lambda_{\min}$ be the ordered eigenvalues of $\mathbf{S}$, and $\hat{\lambda}_1 \ge \dots \ge \hat{\lambda}_d$ those of $\hat{\mathbf{S}}$. Finally, assume that the regularization level in Algorithm 1 satisfies $\lambda \le \lambda_k$ for some $1 \le k \le d$. Then, with probability at least $1 - \delta$ it holds

$$(1 - 2\varepsilon)\lambda_i \le \hat{\lambda}_i \le (1 + 2\varepsilon)\lambda_i, \quad \text{for any } 1 \le i \le k, \qquad (32)$$

with $\varepsilon$ given by (22). As a consequence, we have $(1 - 2\varepsilon)^2 \cdot \lambda_i/\lambda_k \le \hat{\lambda}_i/\hat{\lambda}_k \le (1 - 2\varepsilon)^{-2} \cdot \lambda_i/\lambda_k$.

As a simple application of this result, consider the task of “noisy” principal component analysis (PCA), i.e., performing PCA for the unknown covariance matrix $\mathbf{S}$ from the observations $X_1, \dots, X_n$. A common way to approach it is by performing subspace iteration [HMT11, MCJ13, HP14, BDWY16] with the estimated covariance $\hat{\mathbf{S}}$: randomly choose $U^{(0)} \in \mathbb{R}^{d \times k}$, and then iteratively multiply $U^{(t)}$ by $\hat{\mathbf{S}}$ and orthonormalize the result until convergence. The iterate
converges to the projector on the subspace of the top k eigenvalues of $\hat{\mathbf{S}}$ (providing an estimate of the corresponding subspace for $\mathbf{S}$), and its rate of convergence is known to be controlled by the ratio $\hat{\lambda}_k/\hat{\lambda}_{k+1}$.<fnref n="4" /> Hence, if we use the estimate $\hat{\mathbf{S}}$ produced by Algorithm 1 or Algorithm 2, and if n is sufficient to guarantee that $\epsilon < 1/2$, the convergence rate to the top-k eigenspace of $\hat{\mathbf{S}}$ will essentially be the same as that for the exact method and the target subspace of $\mathbf{S}$.

## 6.2 Ridge Regression with Heavy-Tailed Observations

In random design linear regression [HKZ12], one wants to fit the linear model $Y = X^\top w$ from i.i.d. observations $(X_i, Y_i) \in \mathbb{R}^d \times \mathbb{R}$, $1 \le i \le n$. More previsely, the goal is to find a minimizer $w^* \in \mathbb{R}^d$ of the quadratic risk $L(w) := \mathbb{E}(Y - X^\top w)^2$, where the expectation is over the test pair $(X, Y)$ independent of the sample and coming from the same distribution. In ordinary ridge regression, one fixes the regularization level $\lambda \ge 0$, and estimates $w^*$ with the minimizer of the regularized empirical risk $\tilde{L}_\lambda(w) := \frac{1}{n} \sum_{i=1}^n (Y_i - X_i^\top w)^2 + \lambda \|w\|^2$. The case $\lambda = 0$ corresponds to the ordinary least-squares estimator, while $\lambda > 0$ allows for some bias.

Here we propose a couterpart of this estimator with a favorable statistical guarantee under fourth-moment assumptions on the design and response. Given the sample $X_1, \dots, X_{2n}$, we first compute the covariance estimator $\hat{\mathbf{S}}$ by feeding the hold-out sample $X_{n+1}, \dots, X_{2n}$ to Algorithm 1 with $\theta = \theta_*$ (one could also use Algorithm 2). Then, using the first half of the observations, we construct the “pseudo-decorrelated” observations $(\hat{Z}_1, \dots, \hat{Z}_n)$ with $\hat{Z}_i = \hat{\mathbf{S}}_\lambda^{-1/2} X_i Y_i$, and compute the estimator

$$ \bar{w}_{\lambda} = \hat{\mathbf{S}}_{\lambda}^{-1/2} \bar{Z}, \quad \text{where} \quad \bar{Z} = \frac{1}{n} \sum_{i=1}^{n} \rho_{\bar{\theta}}(\|\hat{Z}_i\|^{1/2}) \hat{Z}_i. \qquad (33) $$

Here, $\rho_{\bar{\theta}}(\cdot)$ is as in Eq. (6)-(7), and $\bar{\theta}$ is defined later. We prove the following result (see Appendix D).

**Theorem 6.1.** In the above setting, assume that $X$ satisfies $\mathbb{E}[X] = 0$, $\mathbb{E}[X \otimes X] = \mathbf{S}$, and assumption (HT), and that $Y$ has finite second and fourth moments: $\mathbb{E}[Y^2] \le v^2$, $\mathbb{E}[Y^4] \le \kappa^4 v^4$. Assume also that $n$ satisfies Eq. (21) from the premise of Thm. 4.1. Then, the estimator $\bar{w}_{\lambda}$ given by (33) with $\bar{\theta} = \sqrt{n\kappa^2\kappa^2v^2\mathrm{df}_{\lambda}(\mathbf{S})/\log(1/\delta)}$ with probability at least $1 - \delta$ satisfies

$$ L(\bar{w}_{\lambda}) - L(w^*) \le O(1) \left[ (\kappa^4 + \kappa^2\kappa^2) \frac{v^2 \mathrm{df}_{\lambda}(\mathbf{S}) \log(2d/\delta)}{n} + \lambda^2 \|\mathbf{S}_{\lambda}^{-1/2} w^*\|^2 \right]. \quad (34) $$

In the above result, the bias term is correct (leading to the minimax-optimal rates in the fixed design setting), and the stochastic term has the asymptotically optimal scaling $O(d_{\mathrm{eff}} \log(1/\delta)/n)$, see, e.g., [CDV07]. However, the obtained bound depends on the second moment $v^2$ of the response instead of its variance. We believe that this problem could be resolved, leading to the fully optimal result, by replacing the truncated estimator $\bar{Z}$ in Eq. (33) with median-of-means.

## 7 Conclusion

In this work, we have provided an estimator of the covariance matrix of a heavy-tailed multivariate distribution that admits a high-probability bound of the type

$$ (1 - \epsilon)\mathbf{S} \preceq \hat{\mathbf{S}} \preceq (1 + \epsilon)\mathbf{S}. $$

<footnote n="4">⁴One can use $U^{(t)} \in \mathbb{R}^{d \times r}$, $r \ge k$; the convergence rate is then controlled by the ratio of non-sequential eigenvalues.</footnote>
The novel estimator is computationally efficient, and has applications in principal component analysis, and in ridge regression with heavy-tailed random design. Let us now point out possible directions for future work.

First, one could investigate even weaker moment assumptions than (HT), for example, assuming the existence of the $2 + \varepsilon$ moment of the one-dimensional marginals of $X$. We envision no principal obstacles in extending our work in this direction.

Second, it would be interesting to reach full optimality in ridge regression with heavy-tailed design, replacing the second moment $v^2$ of the response with its variance $\sigma^2$ (see Eq. (34)). To the best of our understanding, a somewhat different expansion of the excess risk would be needed to achieve this; however, we believe that all tools needed to prove such a result are already in place, and it only remains to combine them in a right way.

## Acknowledgments

The first author has been supported by the ERCIM Alain Bensoussan Fellowship and the ERC grant SEQUOIA 724063. The second author acknowledges support from the ERC grant SEQUOIA 724063. We thank Zaid Harchaoui, Anatoli Juditsky, and Francis Bach for fruitful discussions. Finally, we thank Nikita Zhivotovskiy and anonymous reviewers for their insightful remarks, as well as for pointing out some relevant literature.
# A Degrees of Freedom Lemma

**Lemma A.1.** For any $\mathbf{S} \succcurlyeq 0$ and $\lambda \ge 0$, define $\mathrm{df}_{\lambda}(\mathbf{S}) := \mathrm{Tr}(\mathbf{S}\mathbf{S}_{\lambda}^{-1})$. Then, $\mathrm{df}_{\lambda/2}(\mathbf{S}) \le 2\mathrm{df}_{\lambda}(\mathbf{S})$.

*Proof.* We have

$$ 
\begin{aligned}
|\mathrm{df}_{\lambda/2}(\mathbf{S}) - \mathrm{df}_{\lambda}(\mathbf{S})| &= |\mathrm{Tr}[\mathbf{S}(\mathbf{S}_{\lambda}^{-1} - \mathbf{S}_{\lambda/2}^{-1})]| \\
&= |\mathrm{Tr}[\mathbf{S}_{\lambda}^{-1/2}\mathbf{S}\mathbf{S}_{\lambda}^{-1/2}(\mathbf{I} - \mathbf{S}_{\lambda}^{1/2}\mathbf{S}_{\lambda/2}^{-1}\mathbf{S}_{\lambda}^{1/2})]| \\
&\le \mathrm{df}_{\lambda}(\mathbf{S}) \|\mathbf{I} - \mathbf{S}_{\lambda}^{1/2}\mathbf{S}_{\lambda/2}^{-1}\mathbf{S}_{\lambda}^{1/2}\| \le \mathrm{df}_{\lambda}(\mathbf{S}),
\end{aligned}
$$

where we first used commutativity of the trace, and then the fact (following from the trace Hölder inequality) that $|\mathrm{Tr}(\mathbf{A}\mathbf{B})| \le \|\mathbf{B}\| \mathrm{Tr}(\mathbf{A})$ for $\mathbf{A} \succcurlyeq 0$ and $\mathbf{B}$ with compatible dimensions. The claim follows. □

# B Proof of Lemma 4.1

1°. We start by deriving the consequences of (15). First, note that

$$ \mathbf{S}_{\lambda}^{1/2} \hat{\mathbf{S}}_{\lambda}^{-1} \mathbf{S}_{\lambda}^{1/2} = \mathbf{S}_{\lambda}^{1/2} [\mathbf{S}_{\lambda} - (\mathbf{S} - \hat{\mathbf{S}})]^{-1} \mathbf{S}_{\lambda}^{1/2} = [\mathbf{I} - \mathbf{S}_{\lambda}^{-1/2}(\mathbf{S} - \hat{\mathbf{S}})\mathbf{S}_{\lambda}^{-1/2}]^{-1}. $$

Whence, using (15) and the similarity rules,

$$ \frac{2}{3} \mathbf{I} \preceq \mathbf{S}_{\lambda}^{1/2} \hat{\mathbf{S}}_{\lambda}^{-1} \mathbf{S}_{\lambda}^{1/2} \preceq 2\mathbf{I}. \qquad (35) $$

By the properties of the spectral norm, this implies

$$ 
\begin{aligned}
\| \mathbf{S}_{\lambda}^{1/2} \hat{\mathbf{S}}_{\lambda}^{-1/2} \|^2 &= \| \hat{\mathbf{S}}_{\lambda}^{-1/2} \mathbf{S}_{\lambda}^{1/2} \|^2 \le 2; \\
\| \mathbf{S}_{\lambda}^{-1/2} \hat{\mathbf{S}}_{\lambda}^{1/2} \|^2 &= \| \hat{\mathbf{S}}_{\lambda}^{1/2} \mathbf{S}_{\lambda}^{-1/2} \|^2 \le \frac{3}{2}.
\end{aligned}
\qquad (36) $$

Using that, and proceeding as in the proof of Lemma A.1, we can bound the degrees of freedom surrogate $\mathrm{Tr}(\mathbf{S}\hat{\mathbf{S}}_{\lambda}^{-1})$ in terms of the true quantity $\mathrm{df}_{\lambda}(\mathbf{S}) = \mathrm{Tr}(\mathbf{S}\mathbf{S}_{\lambda}^{-1})$:

$$ 
\begin{aligned}
\mathrm{Tr}(\mathbf{S}\hat{\mathbf{S}}_{\lambda}^{-1}) - \mathrm{df}_{\lambda}(\mathbf{S}) &= \mathrm{Tr}[\mathbf{S}(\mathbf{S}_{\lambda}^{-1} - \hat{\mathbf{S}}_{\lambda}^{-1})] \\
&= \mathrm{Tr}[\mathbf{S}\mathbf{S}_{\lambda}^{-1/2}(\mathbf{I} - \mathbf{S}_{\lambda}^{1/2}\hat{\mathbf{S}}_{\lambda}^{-1}\mathbf{S}_{\lambda}^{1/2})\mathbf{S}_{\lambda}^{-1/2}] \\
&= \mathrm{Tr}[\mathbf{S}_{\lambda}^{-1/2}\mathbf{S}\mathbf{S}_{\lambda}^{-1/2}(\mathbf{I} - \mathbf{S}_{\lambda}^{1/2}\hat{\mathbf{S}}_{\lambda}^{-1}\mathbf{S}_{\lambda}^{1/2})],
\end{aligned}
$$

where in the third line we used commutativity of the trace. Applying the trace Hölder inequality as in the proof of Lemma A.1, we obtain

$$ |\mathrm{Tr}(\mathbf{S}\hat{\mathbf{S}}_{\lambda}^{-1}) - \mathrm{df}_{\lambda}(\mathbf{S})| \le \|\mathbf{I} - \mathbf{S}_{\lambda}^{1/2}\hat{\mathbf{S}}_{\lambda}^{-1}\mathbf{S}_{\lambda}^{1/2}\| \mathrm{df}_{\lambda}(\mathbf{S}) \le 3\mathrm{df}_{\lambda}(\mathbf{S}), \qquad (37) $$

where we combined the triangle inequality with the right-hand side of (35).

2°. We now invoke the results of [WM17] (in what follows, the expectation is conditioned on $\hat{\mathbf{S}}_{\lambda}$). Note that conditionally on $\hat{\mathbf{S}}$, the random vectors $Z_j = \hat{\mathbf{S}}_{\lambda}^{-1/2}X_j$ are i.i.d. with mean zero and covariance $\hat{\mathbf{J}} = \hat{\mathbf{S}}_{\lambda}^{-1/2}\mathbf{S}\hat{\mathbf{S}}_{\lambda}^{-1/2}$. By (9), and using the linear invariance of (HT),

$$ \|\mathrm{E}[\|Z_j\|^2 Z_j \otimes Z_j]\| \le \kappa^4 \|\hat{\mathbf{J}}\|^2 \mathrm{Tr}(\hat{\mathbf{J}}) = \kappa^4 \|\hat{\mathbf{J}}\| \mathrm{Tr}(\hat{\mathbf{J}}). $$
Using (35) and (37) to bound $||\hat{\mathbf{J}}||$ and Tr($\hat{\mathbf{J}}$) correspondingly, this results in

$$ ||\mathbb{E}||Z_j||^2 Z_j \otimes Z_j|| \le 8\kappa^4 \mathrm{df}_\lambda(\mathbf{S}). \quad (38) $$

On the other hand, the estimator $\hat{\mathbf{S}}^{(+)}$ defined in (16) satisfies

$$ \hat{\mathbf{S}}_{\lambda}^{-1/2} \hat{\mathbf{S}}^{(+)} \hat{\mathbf{S}}_{\lambda}^{-1/2} = \frac{1}{m} \sum_{j=1}^{m} \rho_{\theta}(\|Z_j\|) Z_j \otimes Z_j, $$

that is, $\hat{\mathbf{S}}_{\lambda}^{-1/2} \hat{\mathbf{S}}^{(+)} \hat{\mathbf{S}}_{\lambda}^{-1/2}$ is precisely the Wei-Misnker estimator (cf. (7)) of $\hat{\mathbf{J}}$, computed from the sample $(Z_1, ..., Z_m)$. Hence, combining the result of Theorem 2.2 with (38), we see that whenever

$$ \theta \ge 2\sqrt{2\kappa^2} \sqrt{\frac{m \mathrm{df}_\lambda(\mathbf{S})}{\log(2d/\delta)}}, $$

with conditional probability at least $1 - \delta$ it holds

$$ \|\hat{\mathbf{S}}_{\lambda}^{-1/2}(\hat{\mathbf{S}}^{(+)} - \mathbf{S})\hat{\mathbf{S}}_{\lambda}^{-1/2}\| \le \frac{2\theta \log(2d/\delta)}{m}. $$

Finally, we arrive at (17) by writing

$$ \|\mathbf{S}_{\lambda/2}^{-1/2}(\hat{\mathbf{S}}^{(+)} - \mathbf{S})\mathbf{S}_{\lambda/2}^{-1/2}\| \le \|\mathbf{S}_{\lambda/2}^{-1/2}\mathbf{S}_{\lambda}^{1/2}\|^2 \|\mathbf{S}_{\lambda}^{-1/2}\hat{\mathbf{S}}_{\lambda}^{1/2}\|^2 \|\hat{\mathbf{S}}_{\lambda}^{-1/2}(\hat{\mathbf{S}}^{(+)} - \mathbf{S})\hat{\mathbf{S}}_{\lambda}^{-1/2}\|, $$

noting that $\|\mathbf{S}_{\lambda/2}^{-1/2}\mathbf{S}_{\lambda}^{1/2}\|^2 = \|\mathbf{S}_{\lambda/2}^{-1}\mathbf{S}_{\lambda}\| \le 2$, and bounding $\|\mathbf{S}_{\lambda}^{-1/2}\hat{\mathbf{S}}_{\lambda}^{1/2}\|^2 \le 3/2$ via (36). $\square$

## C Proof of Theorem 5.1 and Corollary 5.1

Let us call the truncation level $\theta_j = \theta_{\min}2^j$, with $j \in \mathcal{J}$, admissible if it satisfies the condition in (28), so that $\theta_j$ is the smallest such level. Let $j^* \in \mathcal{J}$ be the minimal $j \in \mathcal{J}$ such that $\theta_{j^*} \ge \theta^*$ for $\theta^*$ defined in Eq. (18); note that this is always possible by the definition (26), and we have

$$ \theta_{j^*} \le 2\theta^*. \quad (39) $$

1º. Let us prove that $\theta_{j^*}$ is admissible with probability at least $1 - \delta$. Indeed, due to (29), the premise (19) of Theorem 4.1 holds for any $\theta = \theta_j$ with $j \in \mathcal{J}$ (recall that $\theta_j \le 2\theta_{\max}$). On the other hand, the premise (18) of Theorem 4.1 holds whenever $\theta_j \ge \theta^*$. Hence the bound (20) of Theorem 4.1 holds for all $\theta = \theta_j$ with $j \ge j^*$, and by the union bound we get that with probability at least $1 - \delta$,

$$ \|\mathbf{S}_{\lambda}^{-1/2}(\hat{\mathbf{S}}_{j} - \mathbf{S})\mathbf{S}_{\lambda}^{-1/2}\| \le \varepsilon_j = \frac{24\theta_j \sqrt{\log(4qd|\mathcal{J}|/\delta) \log(4d|\mathcal{J}|/\delta)}}{n}, \quad \forall j \ge j^*, \quad j \in \mathcal{J}. \quad (40) $$

Moreover, from (29) we also obtain

$$ \varepsilon_j \le \frac{1}{2} \sqrt{\frac{\log(4d|\mathcal{J}|/\delta)}{\mathbf{q} \log(4qd|\mathcal{J}|/\delta)}} \le \frac{1}{2}, \quad j \in \mathcal{J}. $$

Whence for any $j \in \mathcal{J}$ such that $j \ge j^*$ we have, under the event (40), and denoting $\hat{\mathbf{S}}_{j,\lambda} = \hat{\mathbf{S}}_j + \lambda \mathbf{I}$,

$$ \begin{aligned} \|\hat{\mathbf{S}}_{j,\lambda}^{-1/2}(\hat{\mathbf{S}}_j - \hat{\mathbf{S}}_{j^*})\hat{\mathbf{S}}_{j,\lambda}^{-1/2}\| &\le \|\hat{\mathbf{S}}_{j,\lambda}^{-1/2}\mathbf{S}_{\lambda}^{1/2}\|^2 \cdot \|\mathbf{S}_{\lambda}^{-1/2}(\hat{\mathbf{S}}_j - \hat{\mathbf{S}}_{j^*})\mathbf{S}_{\lambda}^{-1/2}\| \\ &\le \|\hat{\mathbf{S}}_{j,\lambda}^{-1/2}\mathbf{S}_{\lambda}^{1/2}\|^2 \left( \|\mathbf{S}_{\lambda}^{-1/2}(\hat{\mathbf{S}}_j - \mathbf{S})\mathbf{S}_{\lambda}^{-1/2}\| + \|\mathbf{S}_{\lambda}^{-1/2}(\hat{\mathbf{S}}_{j^*} - \mathbf{S})\mathbf{S}_{\lambda}^{-1/2}\| \right) \\ &\le 2(\varepsilon_j + \varepsilon_{j^*}) \end{aligned} $$
where we used (40) to bound the terms in the parentheses, and also used (cf. (36) in Appendix B):

$$ \| \hat{\mathbf{S}}_{j,\lambda}^{-1/2} \mathbf{S}_{\lambda}^{1/2} \|^2 \le \frac{1}{1 - \varepsilon_j} \le 2, \quad \forall j \ge j_* $$

Thus, $j_*$ is indeed admissible with probability $\ge 1 - \delta$.

**2°.** Whenever $j_*$ is admissible, we have $\hat{j} \le j_*$, whence $\varepsilon_{\hat{j}} \le \varepsilon_{j_*}$ using that $\varepsilon_j$ increases in $j$. Thus, with probability at least $1 - \delta$ it holds

$$ \begin{aligned} & \| \mathbf{S}_{\lambda}^{-1/2} (\hat{\mathbf{S}}_{\hat{j}} - \mathbf{S}) \mathbf{S}_{\lambda}^{-1/2} \| \\ & \le \| \mathbf{S}_{\lambda}^{-1/2} \hat{\mathbf{S}}_{j_*, \lambda}^{1/2} \|^2 \cdot \| \hat{\mathbf{S}}_{j_*, \lambda}^{-1/2} (\hat{\mathbf{S}}_{\hat{j}} - \mathbf{S}) \hat{\mathbf{S}}_{j_*, \lambda}^{-1/2} \| \\ & \le \| \mathbf{S}_{\lambda}^{-1/2} \hat{\mathbf{S}}_{j_*, \lambda}^{1/2} \|^2 \cdot \left( \| \hat{\mathbf{S}}_{j_*, \lambda}^{-1/2} (\hat{\mathbf{S}}_{j_*} - \hat{\mathbf{S}}_{\hat{j}}) \hat{\mathbf{S}}_{j_*, \lambda}^{-1/2} \| + \| \hat{\mathbf{S}}_{j_*, \lambda}^{-1/2} (\hat{\mathbf{S}}_{j_*} - \mathbf{S}) \hat{\mathbf{S}}_{j_*, \lambda}^{-1/2} \| \right) \\ & \le \frac{3}{2} (2(\varepsilon_{\hat{j}} + \varepsilon_{j_*}) + \varepsilon_{j_*}) \le \frac{15}{2} \varepsilon_{j_*}, \end{aligned} \quad (41) $$

where in order to obtain the last line we used (cf. (36) in Appendix B) that

$$ \| \mathbf{S}_{\lambda}^{-1/2} \hat{\mathbf{S}}_{j_*, \lambda}^{1/2} \|^2 \le 1 + \varepsilon_{j_*} \le 3/2. $$

Finally, combining this with the expression for $\theta_*$ in (18), and using (39)-(41), we arrive at the claimed bound. □

**Proof of Corollary 5.1.** First, $\theta_{\max}$ defined in the premise satisfies Eq. (29) by construction; moreover, Eq. (29) is satisfied as an equality. On the other hand, by simple algebra Eq. (30) guarantees that $\theta_{\max} \ge \theta_*$ for $\theta_*$ defined in Eq. (18). Finally, verifying that $\theta_{\min} \le \theta_*$ is trivial using that $\kappa \ge 1$ and $\mathrm{df}_{\lambda}(\mathbf{S}) \ge 1$. □

# D Proof of Theorem 6.1

1°. First of all, note that $n$ satisfying (21) suffices to guarantee that

$$ \| \mathbf{S}_{\lambda}^{-1/2} (\hat{\mathbf{S}} - \mathbf{S}) \mathbf{S}_{\lambda}^{-1/2} \| \le 48\kappa^2 \sqrt{\frac{\mathrm{df}_{\lambda}(\mathbf{S}) \log(2d/\delta)}{n}} \le \frac{1}{2} \quad (42) $$

holds with probability at least $1 - \delta/2$, cf. (22). Note also that $\hat{\mathbf{S}}_{\lambda}$ is independent from $(X_1, \dots, X_n)$, hence the vectors $\hat{Z}_i = \hat{\mathbf{S}}_{\lambda}^{-1/2} X_i Y_i$, $1 \le i \le n$, are independent when conditioned on $(X_{n+1}, \dots, X_{2n})$. Finally, the conditional to the hold-out sample $X_{n+1}, \dots, X_{2n}$ expectation of $\hat{Z}_i$ is

$$ \hat{\mathbb{E}}[\hat{Z}_i] = \hat{\mathbf{S}}_{\lambda}^{-1/2} \mathbf{S} w_*, \quad (43) $$

where we used that the residual $\xi = Y - X^\top w_*$ satisfies $\mathbb{E}[\xi X] = 0$, which follows from the fact that $w_*$ minimizes $L(w)$.

**2°.** We now decompose the excess risk of $\bar{w}_{\lambda}$ as follows:

$$ \sqrt{L(\bar{w}_{\lambda}) - L(w^*)} \le \underbrace{\| \mathbf{S}^{1/2} (\bar{w}_{\lambda} - \hat{w}_{\lambda}) \|}_{E_1} + \underbrace{\| \mathbf{S}^{1/2} (\hat{w}_{\lambda} - w_{\lambda}) \|}_{E_2} + \underbrace{\| \mathbf{S}^{1/2} (w_{\lambda} - w^*) \|}_{E_3}, \quad (44) $$

where $w_{\lambda}$, given by

$$ w_{\lambda} = \mathbf{S}_{\lambda}^{-1} \mathbf{S} w^*, \quad (45) $$
is the minimizer of $L_{\lambda}(w) = L(w) + \lambda \|w\|^2$, and $\hat{w}_{\lambda} := \hat{\mathbb{E}}[\hat{\mathbf{S}}_{\lambda}^{-1/2} \hat{\mathbf{Z}}_1]$ can be calculated using (43):

$$ \hat{w}_{\lambda} = \hat{\mathbf{S}}_{\lambda}^{-1} \mathbf{S} w^{*}. \qquad (46) $$

The easiest to control in (44) is the term $E_3$ corresponding to the squared bias in the fixed-design setting:

$$ \|\mathbf{S}^{1/2}(w_{\lambda} - w^{*})\| \le \|\mathbf{S}_{\lambda}^{1/2}(w_{\lambda} - w^{*})\| = \lambda \|\mathbf{S}_{\lambda}^{-1/2} w^{*}\|, \qquad (47) $$

resulting in the second term in the brackets in (34).

3º. On the other hand, using (45)-(46) we have

$$ \begin{aligned} E_2 &\le \|\mathbf{S}_{\lambda}^{1/2}(\mathbf{S}_{\lambda}^{-1} - \hat{\mathbf{S}}_{\lambda}^{-1})\mathbf{S}w_{*}\| \\ &= \|\mathbf{S}_{\lambda}^{1/2}\hat{\mathbf{S}}_{\lambda}^{-1}(\mathbf{S} - \hat{\mathbf{S}})\mathbf{S}_{\lambda}^{-1}\mathbf{S}w_{*}\| \qquad (48) \\ &\le \|\mathbf{S}_{\lambda}^{1/2}\hat{\mathbf{S}}_{\lambda}^{-1/2}\|^2 \cdot \|\mathbf{S}_{\lambda}^{-1/2}(\mathbf{S} - \hat{\mathbf{S}})\mathbf{S}_{\lambda}^{-1/2}\| \cdot \|\mathbf{S}_{\lambda}^{-1/2}\mathbf{S}^{1/2}\| \cdot \|\mathbf{S}^{1/2}w^{*}\|, \end{aligned} $$

where the last inequality can be verified by removing the norms. Under the event (42), we can bound the first term by a constant (see the proof of Lemma B in Appendix B), and the second term by

$$ O(1)\kappa^2 \sqrt{\frac{\mathrm{df}_{\lambda}(\mathbf{S}) \log(2d/\delta)}{n}}, $$

cf. (42). The third term is at most one. Finally, we have $\|\mathbf{S}^{1/2}w^*\|^2 = \mathbb{E}[(X^\top w^*)^2] \le \mathbb{E}[Y^2] = v^2$. Collecting the above, under the event (42) we have

$$ E_2 \le O(1)\kappa^2 \sqrt{\frac{v^2 \mathrm{df}_{\lambda}(\mathbf{S}) \log(2d/\delta)}{n}}. $$

4º. Finally, let us estimate the term $E_1$ which corresponds to the additive noise, and delivers the first term in the brackets in (34). Note that we can bound

$$ \begin{aligned} E_1 &= \|\mathbf{S}^{1/2}(\bar{w}_{\lambda} - \hat{w}_{\lambda})\| \\ &\le \|\mathbf{S}_{\lambda}^{1/2}(\bar{w}_{\lambda} - \hat{w}_{\lambda})\| \\ &\le \|\mathbf{S}_{\lambda}^{1/2}\hat{\mathbf{S}}_{\lambda}^{-1/2}\| \cdot \|\hat{\mathbf{S}}_{\lambda}^{1/2}(\bar{w}_{\lambda} - \hat{w}_{\lambda})\| \\ &= \|\mathbf{S}_{\lambda}^{1/2}\hat{\mathbf{S}}_{\lambda}^{-1/2}\| \cdot \|\bar{Z} - \hat{\mathbb{E}}[\hat{Z}]\|, \quad \text{where } \hat{Z} = \hat{\mathbf{S}}_{\lambda}^{-1/2}XY, \end{aligned} $$

cf. (33) and (43). Recall that under the event (42), the first term in the product is bounded by a constant, and it remains to control the deviations of the estimator $\bar{Z}$ of $\hat{Z}$ from the (conditional) average $\hat{\mathbb{E}}[\hat{Z}]$. To this end, consider the following construction due to [Min18, Sec. 3.3]. For the general matrix $A \in \mathbb{R}^{d_1 \times d_2}$, define its Hermitian dilation

$$ \mathcal{H}(A) = \begin{pmatrix} 0_{d_1 \times d_1} & A \\ A^\top & 0_{d_2 \times d_2} \end{pmatrix}, \qquad (49) $$

and for any $\varphi : \mathbb{R} \to \mathbb{R}$, define the map on the space of $(d_1 + d_2) \times (d_1 + d_2)$ Hermitian matrices:

$$ \varphi(\mathbf{A}) := Q\varphi(\mathbf{\Lambda})Q^{\top} = Q \mathrm{diag}(\varphi(\lambda_1) \cdots \varphi(\lambda_d))Q^{\top}, $$

where $Q\mathbf{\Lambda}Q^{\top}$ is the eigendecomposition of $\mathbf{A}$. In this notation, consider the following estimator of $\mathbb{E}[A] \in \mathbb{R}^{d_1 \times d_2}$ from i.i.d. copies $A_1, \dots, A_n$ of $A$: compute the $(d_1 + d_2) \times (d_1 + d_2)$ Hermitian matrix

$$ \hat{\mathbf{T}} = \frac{1}{n} \sum_{i=1}^{n} \psi_{\bar{\theta}}(\mathcal{H}(A_i)), \qquad (50) $$
where $\psi_{\bar{\theta}}(\cdot)$ is the matrix map corresponding to (6) with truncation level $\bar{\theta}$, and then output the top right block of $\hat{T}$ (i.e., the one corresponding to $A$ in $\mathcal{H}(A)$) as the final estimate. As proved in [Min18, Cor. 3.1], the resulting estimate satisfies

$$ \|\bar{A} - \mathbb{E}[A]\| \le O(1) \frac{\bar{\theta} \log(1/\delta)}{n} $$

with probability at least $1 - \delta$, provided that $\delta \le 1/2$, and

$$ \bar{\theta} = \sqrt{\frac{n\bar{w}}{\log(1/\delta)}} \quad \text{for some } \bar{w} \ge \text{Tr} \mathbb{E}[A \otimes A]. $$

On the other hand, one can verify that this construction reduces to $\bar{Z}$ when estimating $\hat{\mathbb{E}}[\hat{Z}]$ from $\hat{Z}_1, \dots, \hat{Z}_n$ with the same $\bar{\theta}$. Thus, with (conditional) probability $\ge 1 - \delta/2$ it holds

$$ \|\bar{Z} - \hat{\mathbb{E}}[\hat{Z}]\| \le O(1) \frac{\bar{\theta} \log(2/\delta)}{n} $$

whenever $\bar{\theta}$ is taken to be

$$ \bar{\theta} = \sqrt{\frac{n\bar{w}}{\log(1/\delta)}} \quad \text{for some } \bar{w} \ge \text{Tr} \hat{\mathbb{E}}[\hat{Z} \otimes \hat{Z}], $$

which then results in the bound

$$ \|\bar{Z} - \hat{\mathbb{E}}[\hat{Z}]\| \le O(1) \sqrt{\frac{\bar{w} \log(2/\delta)}{n}}. $$

It remains to bound $\bar{w} = \text{Tr} \hat{\mathbb{E}}[\hat{Z} \otimes \hat{Z}]$. Using the trace Hölder inequality, we have

$$ \begin{aligned} \text{Tr} [\hat{\mathbb{E}}[\hat{Z} \otimes \hat{Z}]] &= \text{Tr} [\hat{\mathbf{S}}_{\lambda}^{-1} \mathbb{E}[Y^2 XX^{\top}]] \\ &\le \|\hat{\mathbf{S}}_{\lambda}^{1/2} \hat{\mathbf{S}}_{\lambda}^{-1/2}\| \cdot \text{Tr} [\hat{\mathbf{S}}_{\lambda}^{-1} \mathbb{E}[Y^2 XX^{\top}]], \\ &= \|\hat{\mathbf{S}}_{\lambda}^{1/2} \hat{\mathbf{S}}_{\lambda}^{-1/2}\| \cdot \mathbb{E} [\|Y \hat{\mathbf{S}}_{\lambda}^{-1/2} X\|^2], \end{aligned} $$

where the first term on the right is at most a constant under (42). Finally, under the fourth-moment assumptions in the premise of the theorem, we can bound the last term coordinatewise, using that each coordinate of $\hat{\mathbf{S}}_{\lambda}^{-1/2} X$ is simply the projection of $\hat{\mathbf{S}}_{\lambda}^{-1/2} X$ onto the corresponding coordinate vector, and proceeding via Cauchy-Schwarz:

$$ \mathbb{E} \left[ \| Y \hat{\mathbf{S}}_{\lambda}^{-1/2} X \|^2 \right] \le \varkappa^2 \kappa^2 v^2 \mathbb{E} \left[ \| \hat{\mathbf{S}}_{\lambda}^{-1/2} X \|^2 \right] = \varkappa^2 \kappa^2 v^2 \text{df}_{\lambda}(\mathbf{S}). $$

Combining the previous steps, we obtain the claimed result. $\square$

## References

[AW02] Rudolf Ahlswede and Andreas Winter. Strong converse for identification via quantum channels. *IEEE Transactions on Information Theory*, 48(3):569–579, 2002.

[BDWY16] Maria-Florina Balcan, Simon S. Du, Yining Wang, and Adams W. Yu. An improved gap-dependency analysis of the noisy power method. In *Conference on Learning Theory*, pages 284–309, 2016.
[BL08] Peter J. Bickel and Elizaveta Levina. Regularized estimation of large covariance matrices. *The Annals of Statistics*, 36(1):199-227, 2008.

[BT03] Brendan O. Bradley and Murad S. Taqqu. Financial risk and heavy tails. *Handbook of Heavy-Tailed Distributions in Finance, ST Rachev, ed. Elsevier, Amsterdam*, pages 35-103, 2003.

[Cam80] Norm A. Campbell. Robust procedures in multivariate analysis I: Robust covariance estimation. *Applied statistics*, pages 231-237, 1980.

[CDV07] Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares algorithm. *Foundations of Computational Mathematics*, 7(3):331-368, 2007.

[CZZ10] T. Tony Cai, Cun-Hui Zhang, and Harrison H. Zhou. Optimal rates of convergence for covariance matrix estimation. *Ann. Statist.*, 38(4):2118-2144, 2010.

[FHT01] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. *The elements of statistical learning*, volume 1. Springer series in statistics New York, NY, USA:, 2001.

[GVL12] Gene H. Golub and Charles F. Van Loan. *Matrix computations*, volume 3. JHU Press, 2012.

[Han82] Lars P. Hansen. Large sample properties of generalized method of moments estimators. *Econometrica: Journal of the Econometric Society*, pages 1029-1054, 1982.

[HKZ12] Daniel Hsu, Sham M. Kakade, and Tong Zhang. Random design analysis of ridge regression. *The Journal of Machine Learning Research*, 23(9):1-24, 2012.

[HMT11] Nathan Halko, Per-Gunnar Martinsson, and Joel A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. *SIAM review*, 53(2):217-288, 2011.

[HP14] Moritz Hardt and Eric Price. The noisy power method: A meta algorithm with applications. In *Advances in Neural Information Processing Systems*, pages 2861-2869, 2014.

[Jol02] Ian Jolliffe. Principal component analysis. In *International encyclopedia of statistical science*, pages 1094-1096. Springer, 2002.

[KJ14] Bryan Kelly and Hao Jiang. Tail risk and asset prices. *The Review of Financial Studies*, 27(10):2841-2871, 2014.

[KL17] Vladimir Koltchinskii and Karim Lounici. Concentration inequalities and moment bounds for sample covariance operators. *Bernoulli*, 23(1):110-133, 02 2017.

[Lep91] Oleg V. Lepskii. On a problem of adaptive estimation in Gaussian white noise. *Theory of Probability & Its Applications*, 35(3):454-466, 1991.

[Lou14] Karim Lounici. High-dimensional covariance matrix estimation with missing observations. *Bernoulli*, 20(3):1029-1058, 08 2014.

[LR91] Hendrik P. Lopuhaa and Peter J. Rousseeuw. Breakdown points of affine equivariant estimators of multivariate location and covariance matrices. *The Annals of Statistics*, 19(1):229-248, 1991.
[LT13] Michel Ledoux and Michel Talagrand. *Probability in Banach Spaces: isoperimetry and processes*. Springer Science & Business Media, 2013.

[Mar52] Harry Markowitz. Portfolio selection. *The Journal of Finance*, 7(1):77-91, 1952.

[MCJ13] Ioannis Mitliagkas, Constantine Caramanis, and Prateek Jain. Memory limited, streaming pca. In *Advances in Neural Information Processing Systems*, pages 2886-2894, 2013.

[Min18] Stanislav Minsker. Sub-gaussian estimators of the mean of a random matrix with heavy-tailed entries. *The Annals of Statistics*, 46(6A):2871-2903, 2018.

[MW17] Stanislav Minsker and Xiaohan Wei. Estimation of the covariance structure of heavy-tailed distributions. *arXiv:1708.00502*, 2017.

[MZ18] Shahar Mendelson and Nikita Zhivotovskiy. Robust covariance estimation under $L_4 - L_2$ norm equivalence. *arXiv:1809.10462*, 2018.

[Oli10] Roberto I. Oliveira. Sums of random Hermitian matrices and an inequality by Rudelson. *Electron. Commun. Probab.*, 15(26):203-212, 2010.

[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.

[RD99] Peter J. Rousseeuw and Katrien Van Driessen. A fast algorithm for the minimum covariance determinant estimator. *Technometrics*, 41(3):212-223, 1999.

[Tro12] Joel A. Tropp. User-friendly tail bounds for sums of random matrices. *Foundations of Computational Mathematics*, 12(4):389-434, 2012.

[Tro15] Joel A. Tropp. An introduction to matrix concentration inequalities. *Foundations and Trends in Machine Learning*, 8(1-2):1-230, 2015.

[Ver12] Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices. In *Compressed Sensing: Theory and Applications*, pages 210-268. Cambridge University Press, 2012.

[WM17] Xiaohan Wei and Stanislav Minsker. Estimation of the covariance structure of heavy-tailed distributions. In *Advances in Neural Information Processing Systems*, pages 2859-2868, 2017.