Sub-Gaussian estimators of the mean of a random matrix with heavy-tailed entries

Introduction

Because the only assumption on $X$ is the existence of a second moment, it is natural to call such an estimator “robust” For the classical treatment of robust estimators based on the notion of a breakdown point, we refer the reader to .: it admits strong deviation bounds even for the heavy-tailed distributions that can be used to model outliers in the data. Ideas behind these results have also been extended to empirical risk minimization methods which cover a wide range of statistical applications. Let us emphasize that the aforementioned estimators do not require any assumptions on the “shape” of the distribution, such as unimodality or elliptical symmetry.

Generalizations of univariate results to the case of random vectors and random matrices are not straightforward since element-wise deviation inequalities do not always translate into desired bounds. In some cases, element-wise bounds yield inequalities for the “wrong” norm: for example, estimating each entry of the covariance matrix results in a deviation inequality for the Frobenius norm, while we are frequently interested in the bounds for the operator norm that can be much smaller. An approach which often yields “dimension - free” bounds was proposed in and (using generalizations of the median in higher dimensions); however, to the best of our knowledge, results of these papers are still not sufficient to obtain deviation guarantees in the operator norm that we are mainly interested in. Under more restrictive assumptions on the sequence of random matrices $Y_{1},\ldots,Y_{n}$ (such as $\|Y_{j}\|\leq M$ almost surely for some fixed $M>0$, $j=1,\ldots,n$, where $\|\cdot\|$ stands for the operator norm), behavior of the sample mean $\bar{Y}=\frac{1}{n}\sum_{j=1}^{n}Y_{j}$ has been analyzed with the help of matrix concentration inequalities .

A closely related covariance matrix estimation problem has been extensively studied in the past decades. A comprehensive review is beyond the scope of this introduction, so we will just mention few classical results and more recent work related to the current line of research. Statistical properties of the sample covariance matrix for Gaussian and sub-Gaussian observations have been investigated in detail, see and references therein; under weaker moment assumptions, sample covariance estimator has been studied in . Some popular robust estimators of scatter are discussed in , including the Minimum Covariance Determinant (MCD) estimator and the Minimum Volume Ellipsoid estimator (MVE). However, rigorous results for these estimators are available only for elliptically symmetric distributions; see for results on MCD and for results on MVE. Popular Maronna’s and Tyler’s M-estimators of scatter also admit theoretical guarantees for the family of elliptically symmetric distributions, but we are unaware of any results extending beyond this case.

Finally, let us mention that the problem of robust matrix recovery (that is discussed as an example below) has also received attention recently: for instance, the work investigates robust matrix completion under the “low rank + sparse” model. In , authors study low-rank matrix recovery under the assumption that the additive noise has only $(2+\varepsilon)$ moments, and obtain strong results via truncation argument. We propose a different approach based on general techniques developed in this paper and achieve similar results for the matrix completion problem while requiring only the finite variance of the noise.

2 Organization of the paper

Section 2 contains definitions, notation and background material. Our main results are introduced in section 3. After presenting core results, we discuss applications to covariance estimation and low-rank matrix completion in section 4, and illustrate the role of various quantities involved in the general bounds through these examples. Sections 5 and 6 discuss adaptation to unknown parameters that appear in our construction, and contain longer proofs.

Appendix contains proofs of several technical lemmas and results that were omitted in the main text.

Preliminaries

In this section, we introduce main notation and recall several useful facts from linear algebra, matrix analysis and probability theory that we rely on in the subsequent exposition.

Given two self-adjoint matrices $A$ and $B$, we will write $A\succeq B\ (\text{or }A\succ B)$ iff $A-B$ is nonnegative (or positive) definite.

2 Tools from linear algebra

In this section, we collect several facts from linear algebra, matrix analysis and probability theory that are frequently used in our arguments.

Additionally, we will often use the following facts:

Matrix logarithm is operator monotone: if $A\succ 0,\ B\succ 0$ and $A\succeq B$, then $\log(A)\succeq\log(B)$.

Given a fixed self-adjoint matrix $H$, the function

is concave on the cone of positive definite matrices.

See and . Let us mention that Lieb’s theorem is one of the key tools for proving matrix concentration inequalities, and its power in this context was first demonstrated by J. Tropp . ∎

This is a consequence of Peierls inequality, see Theorem 2.9 in and the comments following it. ∎

Since $\mathcal{H}(A)^{2}=\begin{pmatrix}AA^{\ast}&0\\ 0&A^{\ast}A\end{pmatrix},$ it is easy to see that $\|\mathcal{H}(A)\|=\|A\|$. Another tool useful in dealing with rectangular matrices is the following lemma:

Main results

See remark 1 below for examples of such functions. Given $\theta>0$, let $\hat{\mu}_{\theta}$ be such that

(clearly, $\hat{\mu}_{\theta}$ always exists due to monotonicity). Set $\eta=v\sqrt{\frac{2t}{n(1-2t/n)}}$ and $\theta_{\ast}=\sqrt{\frac{2t}{n(v^{2}+\eta^{2})}}$. Assuming that $n>2t$, it is shown in that $|\hat{\mu}_{\theta_{\ast}}-\mu|\leq\eta$ with probability $\geq 1-2e^{-t}$.

where $\theta>0$ is an appropriate constant. It follows from Fact 2.6 that $\widehat{T}^{\ast}_{\theta}$ exists, moreover, it is unique if $\psi(x)$ is strictly increasing. It is also not hard to see that (3.3) is equivalent to

Indeed, if $F_{\psi}(S):=\mbox{tr\,}\sum_{j=1}^{n}\Psi\left(\theta(Y_{j}-S)\right)$, then (3.4) simply states that the gradient of $F_{\psi}$ evaluated at $\widehat{T}^{\ast}_{\theta}$ is equal to zero; see Lemma A.1 in the appendix for more details.

To understand the properties of the estimator defined via (3.3) and (3.4), we will first consider another estimator $\widehat{T}^{(0)}_{\theta}$ that shares many important properties with $\widehat{T}^{\ast}_{\theta}$ but is easier to analyze.

The “preliminary estimator” $\widehat{T}^{(0)}_{\theta}$ is constructed as follows: given $\theta>0$ and a function $\psi$ satisfying (3.1), set $X_{j}:=\psi\left(\theta Y_{j}\right),\ j=1,\ldots,n$ and

Assume that $n$ is large enough and $\theta$ is chosen properly. Then the estimator $\widehat{T}_{\theta}^{\ast}$ defined via (3.4) satisfies the inequality

Most of our results do not depend on the concrete choice of the function $\psi$. One possibility is

Since the latter function is bounded, it can provide additional advantages (such as robustness) in applications. However, note that $\psi_{2}(x)$ does not satisfy (3.1); instead, it satisfies a slightly weaker inequality

hence all subsequent results hold for $\psi_{2}$ as well, albeit with slightly worse constant factors. We also note that both $\psi_{1}$ and $\psi_{2}$ are operator Lipschitz functions; see Lemma A.3 for details.

In this section, we will establish deviation inequalities for the estimator $\widehat{T}_{\theta}^{(0)}=\frac{1}{n\theta}\sum_{j=1}^{n}\psi(\theta Y_{j})$. The lemma below is the cornerstone of our results. As before, given $\theta>0$, let $X_{j}=\psi(\theta Y_{j})$.

It remains to note that by Fact 2.1 and the inequality $\log(1+x)\leq x$ (that holds $\forall\ x>-1$), for all $j=1,\ldots,n$

To establish the second inequality of the lemma, we use the relation $-X_{j}=-\psi(\theta Y_{j})\preceq\log\left(I-\theta Y_{j}+\frac{\theta^{2}}{2}Y_{j}^{2}\right)$ (which follows from (3.1) and Fact 2.1) together with the Fact 2.2 to deduce that

and apply inequality (3.8) to the sequence $-Y_{1},\ldots,-Y_{n}$ with

We are ready to state and prove the main result of this section.

In particular, setting $\theta=\frac{t\sqrt{n}}{\sigma_{n}^{2}}$, we get the “sub-Gaussian” tail bound $2d\exp\left(-\frac{t^{2}}{2\sigma_{n}^{2}/n}\right)$, for a given $t>0$. Alternatively, setting $\theta=\frac{\sqrt{n}}{\sigma_{n}^{2}}$ (independent of $t$), we obtain sub-exponential concentration with tail $2d\exp\left(-\frac{2t-1}{2\sigma_{n}^{2}/n}\right)$ for all $t>1/2$.

where $\widehat{T}_{\theta}^{(0)}$ was defined in (3.5).

As before, set $X_{j}:=\psi\left(\theta Y_{j}\right),\ j=1,\ldots,n$. Then

where we used the second inequality of Lemma 3.1 instead. The result follows by taking $s:=t\sqrt{n}$ since for a self-adjoint matrix $A$, $\|A\|=\max\left(\right.\lambda_{\mbox{\footnotesize{max}\,}}(A),$ $-\lambda_{\mbox{\footnotesize{min}\,}}(A)\left.\right)$. ∎

Sub-Gaussian guarantees provided by Theorem 3.1 hold for a given confidence parameter $t>0$ that has to be fixed a priori: in particular, the optimal value of $\theta$ depends it. However, as it was noted in , this is sufficient to construct (via Lepski’s method ) estimators that admit sub-Gaussian tails uniformly over $t$ in a certain range.

2 Bounds depending on the effective dimension

The bound obtained in Theorem 3.1 explicitly depends on the dimension $d$ of random matrices. Example is subsection 3.2.1 below shows that the dimensional factor in the right-hand side of the inequality is unavoidable in general. However, it is possible to prove a similar inequality which only includes the “effective dimension” defined as

As before, we can set $\theta=\frac{t\sqrt{n}}{\sigma_{n}^{2}}$ to get

For the values of $t\geq\sqrt{\sigma_{n}^{2}/n}$ (when the bound becomes useful), it further simplifies to

For the “sub-exponential regime” with $\theta=\frac{\sqrt{n}}{\sigma_{n}^{2}}$, we get that for all $t\geq\frac{1}{2}\vee\sigma_{n}^{2}/n$ simultaneously,

The argument is similar in spirit to the proof of Theorem 3.1. Details are included in appendix C. ∎

with $f(d)\leq Cd$ for some absolute constant $C$. Since $\left\|\sum_{j=1}^{d}\gamma_{j}e_{j}e_{j}^{T}\right\|=\max\left(|\gamma_{1}|,\ldots,|\gamma_{d}|\right)$, it follows from Lemma A.5 that

for any $0<\tau<1/2$ and some constant $c(\tau)>0$. This shows that the dimensional factor $f(d)$ can not grow slower than $d^{1/2-\tau}$ for any $\tau>0$.

3 Bounds for arbitrary rectangular matrices

and the first inequality follows. To obtain the second inequality, it is enough to use Theorem 3.2 instead of Theorem 3.1 and note that

4 Bounds under weaker moment assumptions

The argument repeats the steps of Lemma 3.1 and Theorem 3.1, the only difference being that application of Fact 2.4 is replaced by Lemma A.2. ∎

Note that for $\alpha=2$, we recover (3.11).

Before we proceed with discussion or further improvements and adaptation issues, let us demonstrate applications of developed techniques to popular problems in statistics and highlight the advantages over existing results.

Examples

We present two examples which highlight the potential improvements obtained via our technique in popular scenarios: estimation of the covariance matrix in Frobenius and operator norms, and low-rank matrix completion problem.

Note that for any matrix $X=\lambda UU^{T}$ of rank $1$ (where $\|U\|_{2}=1$),

Of course, the initial assumption that $\mu$ is known is often unrealistic, hence we modify the estimator as follows. Given $\theta>0$, set

Before presenting the proof, let us make several additional remarks.

It is not hard to show that (see Corollary A.1) that

Construction of $\widehat{\Sigma}_{2n}(\theta)$ essentially halves the effective sample size. While the loss of a constant factor can be deemed insignificant in non-asymptotic theoretical bounds, it is undesirable in applications. A more natural version of the estimator based on a sample of size $2n$ is the U-statistic

Another possibility to avoid “halving” the sample size is to center the data using a robust estimator of location, such as the spatial median or the median-of-means estimator . Analysis of the estimators of these types is not covered in the present paper, and requires a slightly different set of technical tools to deal with dependent summands; see for results in this direction.

2 Estimation of the covariance matrix in Frobenius norm

Next, we present an estimator which achieves strong deviation guarantees in the Frobenius norm. Estimation of the covariance matrix with respect to this norm has been previously investigated in the literature, for instance, see , and references therein; Frobenius norm is a natural choice when one wants to understand the effect of the rank of an unknown covariance matrix on the estimation error . Let $\hat{S}_{2n}$ be the sample covariance estimator based on $Z_{1},\ldots,Z_{2n}$:

The following “soft thresholding” estimator has been studied in ; here, $\tau>0$ is a fixed threshold parameter:

We propose to replace the sample covariance $\hat{S}_{2n}$ by $\widehat{\Sigma}_{2n}$, and consider

It is not hard to see (e.g., see the proof of Theorem 1 in ) that $\widehat{\Sigma}_{2n}^{\tau}$ can be written explicitly as

where $\lambda_{j}(\widehat{\Sigma}_{2n})$ and $v_{j}(\widehat{\Sigma}_{2n})$ are the eigenvalues and corresponding eigenvectors of $\widehat{\Sigma}_{2n}$. The following result holds:

Result stated above mimics the (almost) optimal rates obtained in (in the situation when no data is missing) under significantly weaker assumptions on the underlying distribution.

The proof is based on the following lemma:

Inequality (4.3) holds on the event $\mathcal{E}=\left\{\tau\geq 2\left\|\widehat{\Sigma}_{2n}-\Sigma\right\|\right\}$.

To verify this statement, it is enough to repeat the steps of the proof of Theorem 1 in , replacing each occurrence of the sample covariance $\hat{S}_{2n}$ by its robust counterpart $\widehat{\Sigma}^{\tau}_{2n}$. Result then follows from corollary 4.1 that $\Pr(\mathcal{E})\geq 1-e^{-t}$ whenever $\tau\geq 4\hat{\sigma}\sqrt{\frac{t+\log(2d)}{2n}}$. ∎

3 Matrix completion

To incorporate the structural (low-rank) assumption on $A_{0}$, the following estimator has been considered in the literature: let $\tau>0$, and define

However, strong theoretical guarantees for this estimator exist only when the “noise term” $\xi_{j}$ is either bounded with probability 1, or has sub-exponential tails. We propose to replace $\widehat{A}_{s}$ with a robust estimator

The reasoning behind this choice of $\theta$ is explained below. Consider

Assume that $\xi_{j}$ is independent of $X_{j},\ j=1,\ldots,n$, and that $\mbox{Var}(\xi)<\infty$. For any

By the definition of $\widehat{R}^{\tau}$, we see that

If we replace $\frac{1}{d_{1}d_{2}}\widehat{R}$ by $\frac{1}{d_{1}d_{2}}\widehat{A}_{s}=\frac{1}{n}\sum_{j=1}^{n}Y_{j}\mathcal{H}(X_{j})$, the result follows from Theorem 1 in immediately. To obtain the current statement, it is enough to repeat the argument of Theorem 1 in , replacing each occurrence of the matrix $\frac{1}{d_{1}d_{2}}\widehat{A}_{s}$ by $\frac{1}{d_{1}d_{2}}\widehat{R}$. ∎

To complete the proof, we will estimate each side of the inequality of Lemma 4.2. First, it is obvious from the definition of the Frobenius norm that

It remains to estimate the probability of the event $\mathcal{E}=\left\{\tau\geq 2\|M\|\right\}$. Let

Assume that $\xi_{j}$ is independent of $X_{j}$, $j=1,\ldots,n$. Then

with probability $\geq 1-e^{-t}$. Final result now follows from the combination of this inequality with (4.4), (4.3) and Lemma 4.2.

Optimal choice of θ𝜃\theta and adaptation to the unknown second moment

Parameters $\sigma_{\mbox{\footnotesize{min}\,}}$ and $\sigma_{\mbox{\footnotesize{max}\,}}$ are “crude” preliminary bounds that can differ from $\sigma_{n}/\sqrt{n}$ by several orders of magnitude. Let $\sigma_{j}=\sigma_{\mbox{\footnotesize{min}\,}}2^{j}$ and

be a set of cardinality $|\mathcal{J}|\leq 1+\log_{2}(\sigma_{\mbox{\footnotesize{max}\,}}/\sigma_{\mbox{\footnotesize{min}\,}})$, and for each $j\in\mathcal{J}$ set $\theta_{j}=\theta(j,t)=\sqrt{\frac{2t}{n}}\frac{1}{\sigma_{j}}$. Define

where $\psi(\cdot)$ satisfies (3.1). Finally, set

Next result shows that adaptation is possible at the cost of an additional multiplicative constant factor $6$ in the deviation bound.

The following inequality holds for any $t>0$:

Let $\bar{j}=\min\left\{j\in\mathcal{J}:\ \sigma_{j}\geq\frac{\sigma_{n}}{\sqrt{n}}\right\}$ (hence $\sigma_{\bar{j}}\leq 2\frac{\sigma_{n}}{\sqrt{n}}$). First, we will show that $j_{\ast}\leq\bar{j}$ with high probability. Indeed,

where we used Theorem 3.1 to bound each of the probabilities in the sum. The display above implies that the event

of probability $\geq 1-2d\log_{2}\left(\frac{2\sigma_{\mbox{\footnotesize{max}\,}}}{\sigma_{\mbox{\footnotesize{min}\,}}}\right)e^{-t}$ is contained in $\mathcal{E}=\left\{j_{\ast}\leq\bar{j}\right\}$. Hence, on $\mathcal{B}$ we have that

Let $\sigma_{\mbox{\footnotesize{min}\,}},\sigma_{0,\mbox{\footnotesize{min}\,}}$ and $\sigma_{\mbox{\footnotesize{max}\,}},\sigma_{0,\mbox{\footnotesize{max}\,}}$ be known constants such that

We will first discuss the simplest (but not the most efficient) approach based on splitting the sample $Y_{1},\ldots,Y_{n}$ into two disjoint subsets $G_{1}$ and $G_{2}$ of cardinality $\geq\lfloor n/2\rfloor$ each, and performing one step of the steepest descent. The main advantage of this approach is the fact that it requires very mild assumptions. The idea is to apply Lepski’s method (as discussed in section 5) twice: on the first step, we obtain an estimator $\hat{T}_{0}$ based on subsample $G_{1}$, and on the second step we apply Lepski’s method again to the subsample $\left\{Y_{j}-\hat{T}_{0}:\ 1\leq j\leq n,\ Y_{j}\in G_{2}\right\}$.

Here is the more detailed description: set $\sigma_{j}=2^{j}\sigma_{\mbox{\footnotesize{min}\,}}$,

and $\sigma_{0,j}=2^{j}\sigma_{0,\mbox{\footnotesize{min}\,}}$,

and let $\hat{T}_{0}$ be the “Lepski-type” adaptive estimator based on the subsample $G_{1}$ defined as

$\theta_{j}=\sqrt{\frac{2t}{n/2}}\frac{1}{\sigma_{j}}$, $\psi(\cdot)$ satisfies (3.1) and

$\hat{T}_{1}$ is then defined as follows:

The main feature of this result is the variance term $\sigma_{0}+12\sigma\sqrt{\frac{t}{n}}$ that can be much smaller compared to $\sigma$ as long as $t\ll n$.

We will next show how to design an estimator with deviations controlled by “correct” variance term without sample splitting (however, subject to the condition that the sample size is sufficiently large). In what follows, we will make an additional assumption about the function $\psi$:

For example, we may take $\psi=\psi_{1}$ or $\psi=\psi_{2}$ (see Lemma A.3 for details). As before, let $t>0$ be fixed, set $\sigma_{0,j}=2^{j}\sigma_{0,\mbox{\footnotesize{min}\,}}$,

For all $j\in\mathcal{J}$, define $\delta^{(0)}_{j}=\sigma_{\mbox{\footnotesize{max}\,}}\sqrt{\frac{2t}{n}}$ and

for $k\geq 1$. Next, for each $j\in\mathcal{J}$, we define

for $k\geq 1$. Finally, we apply Lepski’s method to the collection of estimators $\left\{T^{(k)}_{n,j}:j\in\mathcal{J}\right\}$. To this end, define $\hat{T}_{k}:=T^{(k)}_{n,j_{k}^{\ast}}$, where

Note that the estimator $\hat{T}_{k}$ is completely data-dependent. We are ready to state the main result of this section:

where $K>0$ is an absolute constant, and assume that $\tau\leq 1/6$. Moreover, assume that

with probability $\geq 1-8d\left(1+2\log_{2}\left(\frac{12\sigma_{\mbox{\footnotesize{max}\,}}}{5\sigma_{0,\mbox{\footnotesize{min}\,}}}\right)\right)\log_{2}\left(\frac{2\sigma_{0,\mbox{\footnotesize{max}\,}}}{\sigma_{0,\mbox{\footnotesize{min}\,}}}\right)e^{-t}$.

The next corollary easily follows from the preceding result. Let $\mathcal{A}$ be the event of probability

defined in Theorem 6.2. Since by the properties of the steepest descent scheme $T_{n,j}^{(k)}$ converges to the solution (denoted $\widehat{T}_{\theta_{j}}^{\ast}$) of the problem (3.3), we can easily deduce the following inequality.

Let $\left\{\widehat{T}_{\theta_{j}}^{\ast}\right\}_{j\in\mathcal{J}}$ satisfy the equations

One can further apply Lepski’s method (see section 5) to the collection $\left\{\widehat{T}_{\theta_{j}}^{\ast}\right\}_{j\in\mathcal{J}}$ to obtain a completely data-dependent estimator $\widehat{T}^{\ast}$ that satisfies

with high probability (in particular, on event $\mathcal{A}$).

Numerical simulation results

The goal of numerical experiment was to evaluate the quality of estimation of the covariance matrix $\Sigma$ as well as its first eigenvector $e_{1}$ corresponding to $\lambda_{1}=10$. We tested two scenarios with sample sizes equal $n$ to $100$ and $1000$. In both cases, we generated $Z_{1},\ldots,Z_{n}$, i.i.d. copies of $Z$, and centered the data via the spatial (or geometric) median defined as

We compared two estimators, $\widehat{S}_{n}$ and $\widehat{\Sigma}_{n}$ constructed as follows: set $Z_{j}^{0}:=Z_{j}-\widehat{M}_{n}$ for brevity, and

which is the analogue of sample covariance with “robust centering”.

Next, $\widehat{\Sigma}_{n}$ was constructed using a version of Lepski’s method described in section 5. We provide details for completeness: set

and let $\psi(\cdot)$ be the function defined in (3.6). Let $t=\log 10$, and for $j\in\mathcal{J}$, set $\theta_{j}=\sqrt{\frac{2t}{n}}\frac{1}{1.3^{j}}$ and $\hat{\Sigma}_{n,j}=\frac{1}{n\theta_{j}}\sum_{i=1}^{n}\psi\left(\theta_{j}Z_{i}^{0}{Z_{i}^{0}}^{T}\right).$ Finally, define

(note that we modified some constants compared to the “theoretical” version), and finally set $\widehat{\Sigma}_{n}:=\hat{\Sigma}_{n,j_{\ast}}$.

Quality of covariance estimation was evaluated via comparing $\frac{\|\widehat{S}_{n}-\Sigma\|}{\|\Sigma\|}$ with $\frac{\|\widehat{\Sigma}_{n}-\Sigma\|}{\|\Sigma\|}$ over 500 runs of simulations. We also compared errors of estimation of projectors onto the first principal component,

where $u_{1}(\cdot)$ denotes the eigenvector corresponding to the largest eigenvalue of a matrix. Histograms illustrating performance of both estimators are presented in figures 1a and 1b (for the sample size $n=100$), and in figures 2a and 2b (for the sample size $n=1000$). It is clear from the graphs that in all scenarios, $\widehat{\Sigma}_{n}$ performs significantly better than $\widehat{S}_{n}$.

Acknowledgements

I want to thank L. Goldstein, A. Juditsky, A. Nemirovski, as well as the anonymous Referees and the Associate Editor for their insightful suggestions that helped to improve the quality of presentation.

References

Appendix A Supplementary technical results

hence the claim holds for monomials. By linearity, it also holds for arbitrary polynomials. It remains to extend the claim to arbitrary continuously differentiable function via a standard approximation argument (for instance, see [4, chapter 5, section 3]).

Let $1<\alpha\leq 2$ and $c_{\alpha}=\frac{\alpha-1}{\alpha}\vee\sqrt{\frac{2-\alpha}{\alpha}}$. Then $1+y+c_{\alpha}|y|^{\alpha}>0$ and

To check the first claim, it is enough to note that $f(y)=1+y+c_{\alpha}|y|^{\alpha}$ is convex and its minimum is attained for $y_{m}=-\left(\frac{1}{\alpha c_{\alpha}}\right)^{1/(\alpha-1)}$. It is easy to check that $f(y_{m})=1-y_{m}+\frac{y_{m}}{\alpha}$, which implies that $f(y_{m})>0\iff c_{\alpha}>\frac{\alpha-1}{\alpha^{2}}$ which always holds since $c_{\alpha}\geq\frac{\alpha-1}{\alpha}$ and $\alpha>1$.

Choosing $p:=\frac{\alpha}{2(\alpha-1)}$, $q:=\frac{\alpha}{2-\alpha}$, we get $y^{2}\leq\frac{2(\alpha-1)}{\alpha}y^{\alpha}+\frac{2-\alpha}{\alpha}y^{2\alpha}$ which is further bounded above by $2c_{\alpha}y^{\alpha}+c_{\alpha}^{2}y^{2\alpha}$ for $c_{\alpha}=\frac{\alpha-1}{\alpha}\vee\sqrt{\frac{2-\alpha}{\alpha}}$. ∎

Functions $\psi_{1}(x)$ and $\psi_{2}(x)$ defined in Remark 1 are operator Lipschitz, with Lipschitz constants independent of the dimension.

Lipshitz property of $\psi_{1}(x)$ follows from Theorem 1.6.1 in . Result for $\psi_{2}(x)$ follows from Theorem 1.1.1 in the same paper. ∎

For a self-adjoint matrices $R,Q$, $\|R\|\geq\|Q\|$ iff $\|R^{2}\|\geq\|Q^{2}\|$. Clearly,

It implies that $\left\|\begin{pmatrix}S&A\\ A^{\ast}&T\end{pmatrix}^{2}\right\|\geq\left\|S^{2}+AA^{\ast}\right\|\geq\|AA^{\ast}\|$ and $\left\|\begin{pmatrix}S&A\\ A^{\ast}&T\end{pmatrix}^{2}\right\|\geq\left\|T^{2}+A^{\ast}A\right\|\geq\|A^{\ast}A\|$. Since $\begin{pmatrix}0&A\\ A^{\ast}&0\end{pmatrix}^{2}=\begin{pmatrix}AA^{\ast}&0\\ 0&A^{\ast}A\end{pmatrix}$, we obtain

The following lemma is a generalization of Chebyshev’s association inequality.

where the last inequality follows from Lemma A.4 inequality by setting $f\left(V_{1},\ldots,V_{d}\right):=V_{1}^{2}$ and $g\left(V_{1},\ldots,V_{d}\right):=\|V\|_{2}^{2}$. ∎

$\Pr\left(\max\left(|\gamma_{1}|,\ldots,|\gamma_{n}|\right)\geq\left(\frac{1}{2}-\tau\right)\log n\right)\geq c(\tau)>0$ for every $0<\tau<1/2$.

Appendix B Tools from probability theory and linear algebra

We recall several useful results that we will need in the proofs below.

Then $\left\|\sum_{j=1}^{n}Z_{j}\right\|\leq t$ with probability $\geq 1-2d\exp\left(-\frac{t^{2}}{8\sum_{j=1}^{n}M_{j}^{2}}\right)$.

We conclude this section by recalling the notion of Talagrand’s generic chaining complexity (see ) and several related results. Given a metric space $(T,\rho)$, let $\left\{\Delta_{n}\right\}$ be a nested sequence of partitions of $T$ such that ${\rm card}\,\Delta_{0}=1$ and ${\rm card}\,\Delta_{n}\leq 2^{2^{n}}$. For $s\in T$, let $\Delta_{n}(s)$ be the unique subset of $\Delta_{n}$ containing $s$. The generic chaining complexity $\gamma_{2}(T,\rho)$ is defined as

See Theorem 3.2 in for a more general statement. ∎

Appendix C Proof of Theorem 3.2

Define $\phi(x)=\max(e^{x}-1,0)$ and $X_{j}=\psi(\theta Y_{j})$. Proceeding as in the proof of Theorem 3.1, we get that for $s\geq 0$,

Here we have used the fact that $A\preceq B$ implies $SAS^{\ast}\preceq SBS^{\ast}$ for $S=S^{\ast}:=\sqrt{B_{n}^{2}}$, and the equality $\frac{e^{x}-1}{x}=\sum_{j=1}^{\infty}\frac{x^{j-1}}{j!}$. We have shown that

where we used an elementary inequality $\frac{e^{\theta s}}{e^{\theta s}-1}\leq 1+\frac{1}{\theta s}$ on the last step.

Combining the same steps with Fact 2.4 and the equality $-\lambda_{\mbox{\footnotesize{min}\,}}(A)=\lambda_{\mbox{\footnotesize{max}\,}}(-A)$, we get

Finally, replace $s$ by $t\sqrt{n}$ to get the bound in the required form.

Appendix D Proof of Theorem 6.1

Let $\mathcal{E}_{1}$ be the event defined by

In particular, on event $\mathcal{E}_{1}$,

Here, we use the definition of $\widetilde{\Pr}(\cdot)$ on the first step and (D.3) on the second step. The last inequality follows from independence of $G_{2}$ from $\hat{T}_{0}$ (under $\widetilde{\Pr}(\cdot)$) and Theorem 5.1 applied conditionally on $\hat{T}_{0}$: indeed, this can be done since (D.3) holds on $\mathcal{E}_{1}$. It remains to combine the last bound with (D.4) and (D.1). ∎

Appendix E Proof of Theorem 6.2

We will first state several technical results that are required in the proof. Let $j\in\mathcal{J}$ be such that $\sigma_{0,j}=2^{j}\sigma_{0,\mbox{\footnotesize{min}\,}}\geq\sigma_{0}$, and define

which is a consequence of scalar inequality $\log(1+x)\leq x,\ x>-1$ and fact M.2.1, hence we can deduce from fact M.2.2 that

by the definition of $\psi(\cdot)$, we conclude that

Result follows from Theorem M.3.1 and the inequality

which is a consequence of lemma E.2. Indeed,

We are ready to proceed with the proof of the theorem. Let

where $T_{n}^{(0)}$ was defined in M.6.2 as $T_{n}^{(0)}=\frac{1}{n\theta}\sum_{i=1}^{n}\psi\left(\theta Y_{i}\right)$, and note that $\Pr(\mathcal{E}_{0})\geq 1-2de^{-t}$ by Theorem M.3.1. Let

and note that $k_{\max}\leq 1+\left\lfloor\frac{\log_{2}\left(\frac{12\sigma_{\mbox{\footnotesize{max}\,}}}{5\sigma_{0,\mbox{\footnotesize{min}\,}}}\right)}{\log_{2}1.1}\right\rfloor\leq 1+8\log_{2}\left(\frac{12\sigma_{\mbox{\footnotesize{max}\,}}}{5\sigma_{0,\mbox{\footnotesize{min}\,}}}\right)$. Define

Expression under the supremum in (E.2) can be decomposed as follows:

We will treat 3 terms separately: first, it follows from Lemma E.2 that on $\Omega_{j}$

Putting the bounds (E.3),(E.4),(E.5) together, we can estimate the supremum in (E.2) as

Note that we have used bounds $\theta_{j}\sigma_{0}^{2}\leq\sigma_{0}\sqrt{\frac{2t}{n}}$ and $\theta_{j}\left(\delta^{(k-1)}_{j}\right)^{2}\leq\theta_{j}\delta^{(k-1)}_{j}$ (indeed, inequality (6.3) implies that $\delta_{j}^{(m)}\leq 1$ for all $j$ and $m$) to get the second inequality above. Since $j$ was chosen such that $\sigma_{0,j}\geq\sigma_{0}$ and $\tau=1.1K\sqrt{\frac{d^{2}+Lt}{n}}+\sqrt{\frac{2t}{n}}\frac{1}{2\sigma_{0}}\leq\frac{1}{6}$ by assumption, we have shown that

where the last equality follows from the fact that the sequence $\delta_{j}^{(k)}$ defined in (6.1) satisfies the recursive relation

To complete the proof, it is enough to follow the steps of the proof Theorem 5.1 applied to the collection of estimators $\left\{T_{n,j}^{(k)}:\ j\in\mathcal{J}\right\}$: first, let $\bar{j}=\min\left\{j\in\mathcal{J}:\sigma_{0,j}\geq\sigma_{0}\right\}$, and note that the event

has probability $\geq 1-8d\left(1+2\log_{2}\left(\frac{12\sigma_{\mbox{\footnotesize{max}\,}}}{5\sigma_{0,\mbox{\footnotesize{min}\,}}}\right)\right)\log_{2}\left(\frac{2\sigma_{0,\mbox{\footnotesize{max}\,}}}{\sigma_{0,\mbox{\footnotesize{min}\,}}}\right)e^{-t}$. Moreover, on this event $j_{k}^{\ast}\leq\bar{j}$, hence

where we used the fact that $\sigma_{0,\bar{j}}\leq 2\sigma_{0}$ in the last inequality. ∎

To this end, we will use a chaining argument. Recall that the function $\psi(\cdot)$ is operator Lipschitz with Lipschitz constant $L$ by assumption. Recall that $X_{j,i}(S):=\psi\left(\theta_{j}(Y_{i}-S)\right),\ i=1,\ldots,n$. It follows from Assumption 2 (see also Lemma A.3) that for any Hermitian $S_{1},S_{2}$ and $1\leq i\leq n$,

Matrix Hoeffding’s inequality (Lemma B.1) applies with

and $M_{i}=\frac{2L}{n}\|S_{1}-S_{2}\|$, and yields that

where $|C|$ denotes the Lebesgue measure of a set $C$, and $A+C$ stands for the Minkowski sum of the sets $A$ and $C$. For $A=B(r)$, we get $N(B(r),\varepsilon)\leq\frac{\left|B(r+\varepsilon/2)\right|}{\left|B(\varepsilon/2)\right|}$. The volume of the unit ball is given by

with probability $\geq 1-2de^{-t}$. Recall the Dudley’s entropy integral bound (B.1):

Noting that $D(T(\delta_{k-1}),\rho_{d})=2L\,\delta_{k-1}$ and combining Dudley’s bound with the estimate of Lemma E.4, we get

where $C_{1}=\frac{2}{2\sqrt{2}-1}\int_{0}^{1}\log^{1/2}(1+4/\varepsilon)d\varepsilon$. Bound (E.7) implies that with probability $\geq 1-2de^{-t}$,