On the estimation of the mean of a random vector

We study the problem of estimating the mean of a multivariatedistribution based on independent samples. The main result is the proof of existence of an estimator with a non-asymptotic sub-Gaussian performance for all distributions satisfying some mild moment assumptions.


Introduction
Let X be a random vector taking values in R d . We assume throughout the paper that the mean vector µ = EX and covariance matrix Σ = (X − µ)(X − µ) T exist. Given n independent, identically distributed samples X 1 , . . . , X n drawn from the distribution of X, one wishes to estimate the mean vector.
A natural and popular choice is the sample mean (1/n) n i=1 X i that is known to have a near-optimal behavior whenever the distribution is sufficiently light tailed. However, whenever heavy tails are a concern, the sample mean is to be avoided as it may have a suboptimal performance. While the one-dimensional case (i.e., d = 1) is quite well understood (see [3], [5]), various aspects of the multidimensional problem are still to be revealed. This paper aims at contributing to the understanding of the multi-dimensional case.
Before stating the main results, we briefly survey properties of some mean estimators of real-valued random variables. Some of these techniques serve as basic building blocks for the estimators we propose for the vector-valued case.

Estimating the mean of a real-valued random variable
When d = 1, the simplest and most popular mean estimator is the sample mean µ n = (1/n) n i=1 X i . The sample mean is unbiased and the central limit theorem guarantees an asymptotically Gaussian distribution. However, unless the distribution of X has a light (e.g., sub-Gaussian) tail, there are no non-asymptotic sub-Gaussian performance guarantees for µ n . We refer the reader to Catoni [3] for details. However, perhaps surprisingly, there exist estimators of µ with much better concentration properties, see Catoni [3] and Devroye, Lerasle, Lugosi, and Oliveira [5].
A conceptually simple and quite powerful estimator is the so-called median-of-means estimator that has been proposed, in different forms, in various papers, see Nemirovsky and Yudin [14], Hsu [8], Jerrum, Valiant, and Vazirani [10], Alon, Matias, and Szegedy [1]. The median-of-means estimator is defined as follows. Given a positive integer b and x 1 , . . . , x b ∈ R, let q 1/2 denote the median of these numbers, that is, (If several i fit the above description, we take the smallest one.) For any fixed δ ∈ [e 1−n/2 , 1), first choose b = ⌈ln(1/δ)⌉ and note that b ≤ n/2 holds. Next, partition [n] = {1, . . . , n} into b blocks B 1 , . . . , B b , each of size |B i | ≥ ⌊n/b⌋ ≥ 2. Given X 1 , . . . , X n , we compute the sample mean in each block and define the median-of-means estimator by µ (δ) One can show (see, e.g., Hsu [8]) that for any n ≥ 4, where Var(X) denotes the variance of X.
Note that the median-of-means estimator µ (δ) n does not require any knowledge of the variance of X. However, it depends on the desired confidence level δ and the partition B 1 , . . . , B b . Any partition satisfying ∀i, |B i | ≥ ⌊n/b⌋ is valid in order to get (1). Hence, we do not keep the dependence on the partition B 1 , . . . , B b in the notation µ (δ) n . Devroye, Lerasle, Lugosi, and Oliveira [5] introduce estimators that work for a large range of confidence levels under some mild assumptions. Catoni [3] introduces estimators of quite different flavor and gets a non-asymptotic result of the same form as (1). Bubeck, Cesa-Bianchi and Lugosi [2] apply these estimators in the context of bandit problems.

Estimating the mean of random vectors
Consider now the multi-dimensional case when d > 1. The sample mean µ n = (1/n) n i=1 X i is still an obvious choice for estimating the mean vector µ.
If X has a multivariate normal distribution with mean vector µ and covariance matrix Σ, then µ n is also multivariate normal with mean µ and covariance matrix (1/n)Σ and therefore, for δ ∈ (0, 1), with probability at least 1 − δ, where Tr(Σ) and λ max denote the trace and largest eigenvalue of the covariance matrix, respectively (Hanson and Wright [7]). For non-Gaussian and possibly heavy-tailed distributions, one cannot expect such a sub-Gaussian behavior of the sample mean. The main goal of this paper is to investigate under what conditions it is possible to define mean estimators that reproduce a (non-asymptotic) sub-Gaussian performance similar to (2). Lerasle and Oliveira [11], Hsu and Sabato [9], and Minsker [13] extend the median-ofmeans estimator to more general spaces. In particular, Minsker's results imply that for each δ ∈ (0, 1) there exists a mean estimator µ (δ) n and a universal constant C such that, with probability at least 1 − δ, While this bound is quite remarkable-note that no assumption other than the existence of the covariance matrix is made-, it does not quite achieve a sub-Gaussian performance bound that resembles (2). An instructive example is when all eigenvalues are identical and equal to λ max . If the dimension d is large, (2) is of the order of (λ max /n)(d + log(δ −1 )) while (3) gives the order (λ max /n)(d log(δ −1 )). The main result of this paper is the construction of a mean estimator that, under some mild moment assumptions, achieves a sub-Gaussian performance bound in the sense of (2). More precisely, we prove the following.
Theorem 1 For all δ ∈ (0, 1) there exists a mean estimator µ (δ) n and a universal constant C such that if X 1 , . . . , X n are i.i.d. random vectors in R d with mean µ ∈ R d and covariance matrix Σ such that there exists a constant K > 0 such that, for all v ∈ R d with v = 1, then for all n ≥ CK log d (d + log(1/δ)), The theorem guarantees the existence of a mean estimator whose performance matches the sub-Gaussian bound (2), up to the additional term of the order of (1/n)λ max log log d for all distributions satisfying the fourth-moment assumption given above. The additional term is clearly of minor importance. (For example, it is dominated by the first term whenever Tr(Σ) > λ max log log d.) With the estimator we construct, this term is inevitable. On the other hand, the inequality of the theorem only holds for sample sizes that are at least a constant times d log d. This feature is not desirable for truly high-dimensional problems, especially taking into account that Minsker's bound is "dimension-free".
The fourth-moment assumption can be interpreted as a boundedness assumption of the kurtosis of (X − µ) T v. The same assumption has be used in Catoni [4] and Giulini [6] for the robust estimation of the Gram matrix. The fourth-moment assumption may be weakened to an analogous "(2 + ε)-th moment assumption" that we do not detail for the clarity of the exposition.
We prove the theorem by constructing an estimator in several steps. First we construct an estimator that performs well for "spherical" distributions (i.e., for distributions whose covariance matrix has a trace comparable to dλ max ). This estimator is described in Section 2. In the second step, we decompose the space-in a data-dependent way-into the orthogonal sum of O(log d) subspaces such that all but one subspaces are such that the projection of X to the subspace has a spherical distribution. The last subspace is such that the projection has a covariance matrix with a small trace. In each subspace we apply the first estimator and combine them to obtain the final estimator µ (δ) n . The proof below provides an explicit value of the constant C, though no attempt has been made to optimize its value.
The constructed estimator is computationally so demanding that even for moderate values of d it is hopeless to compute it in reasonable time. In this sense, Theorem 1 should be regarded as an existence result. It is an interesting an important challenge to construct estimators with similar statistical performance that can be computed in polynomial time (as a function of n and d). Note that the estimator of Minsker cited above may be computed by solving a convex optimization problem, making it computationally feasible, see also Hsu and Sabato [9] for further computational considerations.

An estimator for spherical distributions
In this section we construct an estimator that works well whenever the distribution of X is sufficiently spherical in the sense that a positive fraction of the eigenvalues of the covariance matrix is of the same order as λ max . More precisely, for c ≥ 1, we call a distribution c-spherical if dλ max ≤ cTr(Σ).
For each δ ∈ (0, 1) and unit vector w n (w) as the median-of-means estimate (as defined in Section 1.1) of w T µ = Ew T X based on the i.i.d. sample w T X 1 , . . . , w T X n .
Let N 1/2 ⊂ S d−1 be a minimal 1/2-cover, that is, a set of smallest cardinality that has the property that for all u ∈ S d−1 there exists w ∈ N 1/2 with u − x ≤ 1/2. It is well known (see, e.g., [12,Lemma 13 Noting that Var(w T X) ≤ λ max , by (1) and the union bound, we have that, with probability at least 1 − δ, In other words, if, for λ > 0, we define the empirical polytope then with probability at least 1 − δ, µ ∈ P δ,λmax . In particular, on this event, P δ,λmax is nonempty. Suppose that an upper bound of the largest eigenvalue of the covariance matrix λ ≥ λ max is available. Then we may define the mean estimator Now suppose that µ ∈ P δ,λ and let y ∈ P δ,λ be arbitrary. Define u = (y −µ)/ y −µ ∈ S d−1 , and let w ∈ N 1/2 be such that w − u ≤ 1/2. (Such a w exists by definition of N 1/2 .) Then where we used Cauchy-Schwarz and the fact that y, µ ∈ P δ,λ . Rearranging, we obtain that, on the event that µ ∈ P δ,λ , provided that λ ≥ λ max . Summarizing, we have proved the following.
The bound we obtained has the same sub-Gaussian form as (2), up to a multiplicative constant, whenever the distribution is c-spherical. To make the estimator fully datadependent, we need to find an estimate λ that falls in the interval [λ max , 2λ max ], with high probability. This may be achieved by splitting the sample in two parts of equal size (assuming n is even), estimating λ max using samples from one part and computing the mean estimate defined above using the other part. In the next section we describe such a method as a part of a more general procedure.

Empirical eigendecomposition
In the previous section we presented a mean estimate that works well for "spherical" distributions. We will use this estimator as a building block in the construction of an estimator that has the desirable performance guarantee for distributions with any covariance matrix. In addition to finite covariances, we assume that there exists a constant K > 0 such that, for all v ∈ R d with v = 1, In this section we assume that n ≥ 2(400e) 2 K log 3/2 d d log 25 + log(2 log 3/2 d) + log(1/δ) . The basic idea is the following. We split the data into two equal halves. We use the first half in order to decompose the space into the sum of orthogonal subspaces such that the projection of X into each subspace is 4-spherical. Then we may estimate the projected means by the estimator of the previous section.
Next we describe how we obtain an orthogonal decomposition of the space based on n i.i.d. observations X 1 , . . . , X n . Let s = ⌈log 3/2 d 2 ⌉ and m = ⌊n/s⌋. Divide the sample into s blocks, each of size at least m. In what follows, we describe a way of sequentially decomposing R d into the orthogonal sum of s + 1 subspaces R d = V 1 ⊕ · · · ⊕ V s+1 . First we construct V 1 using the first block X 1 , . . . , X m of observations. Then we use the second block to build V 2 , and so on, for s blocks. The key properties we need are that (a) the random vector X, projected to any of these subspaces has a 4-spherical distribution; (b) the largest eigenvalue of the covariance matrix of X, projected on V i is at most λ max (2/3) i−1 .
To this end, just like in the previous section, let N γ ⊂ S d−1 be a minimal γ-cover of the unit sphere S d−1 for a sufficiently small constant γ ∈ (0, 1). The value γ = 1/100 is sufficient for our purposes and in the sequel we assume this value. Note that |N γ | ≤ (4/γ) d (see [12,Lemma 13.1.1] for a proof of this fact).
Initially, we use the first block X 1 , . . . , X m . We may assume that m is even. Using these observations, for each u ∈ N γ , we compute an estimate V (δ) m (u) of u T Σu = E(u T (X −µ)) 2 = (1/2)E(u T (X − X ′ )) 2 , where X ′ is an i.i.d. copy of X. We may construct the estimate by forming m/2 i.i.d. random variables (1/2)(u T (X 1 − X m/2+1 )) 2 , . . . , (1/2)(u T (X m/2 − X m )) 2 and estimate their mean by the median-of-means estimate V (δ) m (u) with parameter δ/(s(4/γ) d ). Then (1), together with assumption (4) implies that, with probability at least 1 − δ/s, Our assumptions on the sample size guarantee that ε m < 1/100. The event that the inequality above holds is denoted by By the argument above, Σ ∈ M δ,m on the event E 1 . In particular, on E 1 , M δ,m in nonempty. Define the estimated covariance matrix m and Σ are in M δ,m , on this event, we have Now compute the spectral decomposition where λ 1 ≥ · · · ≥ λ d ≥ 0 are the eigenvalues and v 1 , . . . , v d the corresponding orthogonal eigenvectors. Let u ∈ S d−1 be arbitrary and let v be a point in N γ with smallest distance to u. Then (by Cauchy-Schwarz and using the fact that u − v ≤ γ) In particular, on E 1 we have λ 1 ≤ βλ max where β = 1+εm 1−εm /(1 − 3γ) < 1.1.
By a similar argument, we have that for any u ∈ S d−1 , if v is the point in N γ with smallest distance to u, then on E 1 , In particular, λ max ≤ β λ 1 ≤ (4/3) λ 1 . Similarly, Let d 1 be number of eigenvalues λ i that are at least λ 1 /2 and let V 1 be the subspace of R d spanned by v 1 , . . . , v d 1 . Denote by Π 1 (X) the orthogonal projection of the random variable X (independent of the X i used to build V 1 ) onto V 1 . Then for any u ∈ V 1 ∩ S d−1 , on the event E 1 , by (7), and therefore In particular, the ratio of the largest and smallest eigenvalues of the covariance matrix of Π 1 (X) is at most 4 and therefore the distribution of Π 1 (X) is 4-spherical.
On the other hand, on the event E 1 , for any unit vector u ∈ V ⊥ 1 ∩S d−1 in the orthogonal complement of V 1 , we have u T Σu ≤ 2λ max /3. To see this, note that u T Σ (δ) m u ≤ λ 1 /2 and therefore, denoting by v the point in N γ closest to u, (6), and a similar argument for the last term) In other words, the largest eigenvalue of the covariance matrix of Π ⊥ 1 (X) (the projection of X to the subspace V ⊥ 1 ) is at most (2/3)λ max .
In the next step we construct the subspace V 2 ⊂ V ⊥ 1 . To this end, we proceed exactly as in the first step but now we replace R d by V ⊥ 1 and the sample X 1 , . . . , X m on the first block by the variables Π ⊥ 1 (X m+1 ), . . . , . Just like in the first step, with probability at least 1 − δ/s we obtain a (possibly empty) subspace V 2 , orthogonal to V 1 such that Π 2 (X), the projection of X on V 2 , has a 4-spherical distribution and largest eigenvalue of the covariance matrix of Π ⊥ 2 (X) (the projection of X to the subspace ( We repeat the procedure s times and use a union bound the s events. We obtain, with probability at least 1 − δ, a sequence of subspaces V 1 , . . . , V s , with the following properties: (i) V 1 , . . . , V s are orthogonal subspaces.
(ii) For each i = 1, . . . , s, Π i (X), the projection of X on V i , has a 4-spherical distribution.
(iii) The largest eigenvalue of the covariance matrix of Π i (X) is at most λ (iv) The largest eigenvalue λ (i) 1 of the estimated covariance matrix of Π i (X) satisfies Note that it may happen for some T < s, we have R d = V 1 ⊕ · · · ⊕ V T . In that case we define V T +1 = · · · = V s = ∅.

Putting it all together
In this section we construct our final multivariate mean estimator and prove Theorem 1. To simplify notation, we assume that the sample size is 2n. This only effects the value of the universal constant C in the statement of the theorem.
The data is split into two equal halves (X 1 , . . . , X n ) and (X n+1 , . . . , X 2n ). The second half is used to construct the orthogonal spaces V 1 , . . . , V s as described in the previous section. Let d 1 , . . . , d s denote the dimension of these subspaces. Recall that, with probability at least 1 − δ, the construction is successful in the sense that the subspaces satisfy properties (i)-(iv) described at the end of the previous section. Denote this event by E. In the rest of the argument we condition on (X n+1 , . . . , X 2n ) and assume that E occurs. All probabilities below are conditional.
Our final estimator is µ (δ) n = s+1 i=1 µ s+1 . By the union bound, we have that, with probability at least 1 − δ, ln(e(2 log 3/2 d + 1)/δ) n s i=1 λ (i) 1 +C λ max log((2 log 3/2 d + 1)/δ) n First notice that, by properties (iii) and (iv) at the end of the previous section, On the other hand, since and for i ≤ s each Π i (X) has a 4-spherical distribution, we have that This concludes the proof of Theorem 1.