Optimal Sub-Gaussian Mean Estimation in $\mathbb{R}$

Given a set of i.i.d. samples from a real-valued distribution, the goal is to return, with extremely high probability, an accurate estimate of the distribution’s mean. Specifically, given a sample set $X$ of size $n$ consisting of independent draws from a real-value distribution $D$, an $(\epsilon,\delta)$-estimator of the mean is a function $\hat{\mu}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ such that, except with failure probability $\leq\delta$, the estimate $\hat{\mu}(X)$ is within $\epsilon$ of the true mean $\mu(D)$. Namely,

The goal is to find the optimal tradeoff between the sample size $n$, and the error parameters $\epsilon$ and $\delta$, for the distribution $D$. Fixing any two of the three parameters and minimizing the third yields essentially equivalent reformulations of the problem: we can fix $\epsilon,\delta$ and minimize the sample complexity $n$; we can fix $\delta,n$ and minimize error $\epsilon$; or we can fix $\epsilon,n$ and minimize the failure probability $\delta$ (maximizing the robustness $1-\delta$).

Perhaps the most standard and well-behaved setting for mean estimation is when the distribution $D$ is a Gaussian. The sample mean (the empirical mean) is a provably optimal estimator in our sense when $D$ is Gaussian: for any $\epsilon,\delta>0$, the sample mean $\mu(X)$ is an $(\epsilon,\delta)$-estimator when given a sample set of size $n=(2+o(1))\frac{\sigma^{2}(D)\cdot\log\frac{1}{\delta}}{\epsilon^{2}}$ (all logarithms will be base $e$); and there is no $(\epsilon,\delta)$- estimator for Gaussians if the constant 2 in the previous expression for the sample size is changed to any smaller number.

The main result of this paper is an estimator that performs as well, on any distribution with finite variance, as the sample mean does on a Gaussian, without knowledge of the distribution or its variance:

We have parameterized the above theorem in terms of fixing the sample size $n$ and the robustness parameter $\delta$ and asking for the minimum error $\epsilon$; however, because of the simple functional form of the bounds of Theorem 1, we can equivalently rephrase it as saying that, for any $\epsilon,\delta$, (a reparameterized) Algorithm 1 is an $(\epsilon,\delta)$ estimator using $(2+o(1))\frac{\sigma(D)^{2}}{\epsilon^{2}}\log\frac{1}{\delta}$ samples; or for any $n,\epsilon$, Algorithm 1 gives an estimate that is $\delta=\exp(-\frac{n\epsilon^{2}}{(2+o(1))\cdot\sigma^{2}(D)})$-robust. For each of these formulations, the performance is optimal up to the $1+o(1)$ factor, as evidenced by the well-known Gaussian case, as explained above.

We make the following observations regarding the main (minimal) assumption in the theorem, namely the finiteness of the variance of the unknown distribution. First, imposing further assumptions about the finiteness of higher moments will not yield any improvements to the result, since matching lower bounds are provided by Gaussians, for which all moments are finite. Second, as shown by Devroye et al. [7], relaxing the finite variance assumption by only assuming, say, the finiteness of the $(1+\beta)^{\text{th}}$ moment for some $\beta<1$ will yield strictly worse sample complexity. In particular, the sample complexity will have an $\epsilon$-dependence that is $\omega(1/\epsilon^{2})$. Thus, our result shows that mean estimation can be performed at a sub-Gaussian rate, with the optimal multiplicative constant of $2$ in the sample complexity, if and only if the variance of the underlying distribution is finite.

We also contrast with previous works that attain optimal sub-Gaussian convergence but make additional assumptions such as the finiteness of the kurtosis ($4^{\textrm{th}}$ moment) [4, 7]. The gap in assumptions between those works and this work is not only theoretical, but also of practical consequence: power law distributions are known to be good models of certain real-world phenomena, and for exponents in the range $(3,5]$, the variance exists, but not the kurtosis.

We briefly describe the main features of our estimator, as a setting for what follows, and to distinguish it from prior work. At the highest level: in order to return a $\delta$-robust estimate of the mean, our estimator “throws out the $\frac{1}{3}\log\frac{1}{\delta}$ most extreme points in the samples”, and returns the mean of what remains. More specifically, outliers are thrown out in a weighted manner, where we throw out a fraction of each data point, with the fraction proportional to the square of its distance from a median-of-means initial guess for the mean, where the fraction is capped at 1, and the proportionality constant is chosen so that the total weight thrown out equals exactly $\frac{1}{3}\log\frac{1}{\delta}$. See Estimator 1 for full details, but we stress here that the estimator is simple to implement—it may be computed in linear time—and therefore applicable in practice.

The above description is rather different from the typical M-estimator/$\psi$-estimator approach of Catoni [4] and other works in this area. However, as we see in Section 3, our estimator can be reinterpreted as a 2-parameter $\psi$-estimator, and the proof of our main result will crucially rely on this reformulation.

Perhaps the most non-obvious part of our estimator is throwing out exactly $\frac{1}{3}\log\frac{1}{\delta}$ many samples. We motivate this quantity in this section, by considering the special case of estimating the mean of asymmetric—very biased—Bernoulli distributions, which is in some sense an extremal case for our setting.

Thus, the $\frac{1}{3}\log\frac{1}{\delta}$ term in our estimator arises essentially from a $3^{\textrm{rd}}$ order correction to the sample mean, at least in the special case of Bernoulli distributions. For additional intuition and motivation about the “$3^{\textrm{rd}}$ order correction” in our estimator, please refer to Appendix A.

In addition to settling the fundamental sample complexity question of mean estimation, we point out that the estimator construction and analysis may also be of independent interest. In particular, the analysis framework—as described below—is generalizable to other problem settings and estimator constructions.

Our overall analysis framework may be viewed as a Chernoff bound—showing exponentially small probability of estimation error via bounds on a moment generation function (expectation of an exponentiated real-valued random variable). However, since we seek to analyze our estimator to sub-constant accuracy, many standard approaches fail to yield the required resolution. We point out three crucial components of our approach.

First, our estimator (Estimator 1) is not a sum of independent terms, which is fundamental to standard Chernoff bound approaches, and thus we instead reformulate our estimator as a 2-parameter $\psi$-estimator (see Definition 2). This technique rewrites our estimate $\hat{\mu}$ as the first coordinate of the root $(\hat{\mu},\hat{\alpha})$ of a system of 2 equations $\psi_{\mu}(\hat{\mu},\hat{\alpha})=0$ and $\psi_{\alpha}(\hat{\mu},\hat{\alpha})=0$, where the functions $\psi_{\mu}(\hat{\mu},\hat{\alpha})=\sum_{i}\psi_{\mu}(x_{i},\hat{\mu},\hat{\alpha})$ and $\psi_{\alpha}(\hat{\mu},\hat{\alpha})=\sum_{i}\psi_{\alpha}(x_{i},\hat{\mu},\hat{\alpha})$ are explicitly sums of a corresponding function applied to each of the $n$ independent data points in the sample set. Thus we have bought independence at the price of making the estimator an implicit function, introducing two new variables. One-dimensional estimators of this form are standard: for example, Catoni’s [4] mean estimator in the case of known variance is a (1 parameter) $\psi$-estimator for which he proves finite sample concentration. However, adding another dimension—$\hat{\alpha}$, a new implicit variable whose value the estimator will ultimately discard—is less standard, without standard analysis techniques, yet significantly increases the expressive power of such estimators [24]. Our high-level approach is to find carefully chosen linear combinations of the functions $\psi_{\mu}$ and $\psi_{\alpha}$, each of which is now a sum of independent terms, and prove Chernoff bounds about these linear combinations.

Second, even after identifying these linear combinations of $\psi$ functions, the corresponding Chernoff bound analysis is difficult to directly tackle. The Chernoff bound analysis, as it turns out, is essentially equivalent to bounding a max-min optimization problem where the maximization is over the set of real-valued probability measures with mean 0 and variance 1. In other words, the max-min optimization problem can be interpreted as having uncountably infinitely many variables. In order to drastically simplify the problem and make it amenable to analysis, we use convex-concave programming and linear programming duality techniques to reduce the problem to a pure minimization problem with a small finite number of variables, which we can analyze tightly.

We believe that the above two ideas—1) reformulating an estimator as a multi-parameter $\psi$-estimator, so as to find a proxy of the estimator that is a sum of independent variables, and 2) viewing the corresponding Chernoff bound analysis as an optimization problems and applying relevant duality techniques—form a general analysis framework which expands the space of possible estimators that are amenable to tight analysis.

Consider the expression in Step 3 for the final returned value of the estimator, $\hat{\mu}=\kappa+\frac{1}{n}\sum_{i}(x_{i}-\kappa)(1-\min(\alpha(x_{i}-\kappa)^{2},1))$. Without the final $\min$ expression, the expression $\kappa+\frac{1}{n}\sum_{i}(x_{i}-\kappa)\cdot 1$ computes exactly the sample mean. The factor $(1-\min(\alpha(x_{i}-\kappa)^{2},1))$ may be thought of as a weight on the $i^{\textrm{th}}$ element, between 0 and 1, where a weight of 1 leaves that element as is, but a weight towards 0 essentially throws out part of the sample $x_{i}$ and instead defaults to the median-of-means estimate $\kappa$. Thus, rather than either keeping or discarding each entry, the weight $\min(\alpha(x_{i}-\kappa)^{2},1)$ specifies what fraction of the $i^{\textrm{th}}$ sample to discard.

The condition in Step 2 of Estimator 1 picks $\alpha$ so that the total, weighted, number of discarded samples equals $\frac{1}{3}\log\frac{1}{\delta}$. The expression $\min(\alpha(x_{i}-\kappa)^{2},1)$ specifying what fraction of each $x_{i}$ to discard says, essentially, that this fraction should be proportional to the square of the deviation of $x_{i}$ from the mean estimate $\kappa$, capped at 1 so that we do not discard “more than 100% of” any sample $x_{i}$.

We point out three basic structural properties of Estimator 1 that both shed light on the estimator itself, and will be crucial to its analysis. We formally state and prove these properties in Appendix B.

First, the estimator is “affine invariant” in the sense that, if its input samples $\{x_{i}\}$ undergo an affine map $x\rightarrow ax+b$ then its output will be mapped correspondingly. Second, as is well known, the median-of-means estimate $\kappa$ of Step 1, while not as accurate as what we will eventually return, is robust in the sense that, with probability at least $1-\delta/2$, the median-of-means estimate has additive error from the true mean that is at most $O(\sigma\cdot\sqrt{\frac{\log\frac{1}{\delta}}{n}})$—proportional to the eventual guarantees of our estimator, but with somewhat worse proportionality constant. Third, if we temporarily ignore Step 1, treating $\kappa$ as a free parameter, we show that the final output of the algorithm, $\hat{\mu}$, varies very little with $\kappa$. Combined with the accuracy guarantees of the median-of-means estimate, the difference in the final estimate between using the median-of-means as $\kappa$ versus using the true mean as $\kappa$ is inconsequential (a $o(1)$ factor) compared to the total additive error we aim for. Therefore, for the purposes of analysis, it suffices to assume that $\kappa$ takes the value of the true mean (though an algorithm could not do this in practice, as the true mean is unknown).

These structural properties allow us to drastically simplify the analysis: the affine invariance means it is sufficient to show our estimator works for the special case of distributions with mean 0 and variance 1; the second and third properties mean that errors in $\kappa$ effectively do not matter, and, for distributions with mean 0, it is sufficient to omit Step 1 and instead just analyze the case where $\kappa=0$.

We point out that Estimator 1 when modified to set $\kappa=0$ (independently of the samples) is no longer affine invariant, nor is its reformulation as a $\psi$-estimator in Section 3.3. The structural properties in this section show that, instead of analyzing the actual estimator (Estimator 1 which is affine invariant), it suffices to analyze this artificially simplified, although no longer affine invariant, estimator which sets $\kappa=0$, on distributions with mean 0 and variance 1. Explicitly, in the rest of the paper we will show Proposition 3 (Section 4), which analyzes the mean-0 variance-1 case of the $\psi$-estimator defined below in Definition 2; the discussion of this section—made formal in Appendix B—shows that this proposition implies our main result, Theorem 1.

As discussed in Section 1.4, our estimator, even its simplified version with $\kappa=0$, is not a sum of independent terms, making it difficult to tightly bound its moment generating function, and hence also difficult to prove its concentration around the true mean using a Chernoff-style bound. Our solution is to reformulate Estimator 1, with the simplifying assumption that $\kappa=0$, as a 2-parameter $\psi$-estimator, as defined in Definition 2. This reformulation defines our estimate $\hat{\mu}$ implicitly in terms of two new functions $\psi_{\mu}$ and $\psi_{\alpha}$ that are indeed sums of $n$ independent terms, each term depending on a single $x_{i}$. We will use this representation crucially for the concentration analysis of the estimator.

Interestingly, in the case that none of the samples are “truncated”, (and $\kappa=0$), the overall output of the estimator becomes $\frac{1}{n}\sum_{i}x_{i}-\alpha x_{i}^{3}=\frac{1}{n}\sum_{i}x_{i}-\frac{\log(1/\delta)}{3n\hat{v}}x_{i}^{3}$, namely, $\hat{\mu}$ is “the empirical mean, corrected by subtracting $\frac{1}{3n}\log\frac{1}{\delta}$ times the ratio of the empirical $3^{\textrm{rd}}$ moment over the empirical $2^{\textrm{nd}}$ moment.”

In order to show Lemma 4, we aim to find bounds on the failure probability that are as strong as possible. Appealing to Lemma 6 that we have just proven, recall that, as in the standard Chernoff bound methodology, we are still free to choose the parameters $d_{\mu},d_{\alpha}$, which we do so as to minimize the resulting bound on the failure probability. Phrased abstractly, the goal is, for the $\hat{\mu},\hat{\alpha}$ of Lemma 4, to show that, for any distribution $D$ of mean 0 and variance 1, there is a choice $\boldsymbol{\mathbf{d}}=(d_{\mu},d_{\alpha})$ that makes Equation 4 sufficiently small. Phrased as an optimization problem, our goal is to evaluate (or tightly bound):

where $D$ ranges over distributions of mean 0 and variance 1.

We will use convex-concave programming and linear programming duality to significantly simplify the max-min program in Equation 5 before we dive into the part of analysis that is ad hoc for this problem. We wish to emphasize here again that the steps of 1) writing an estimator as a multi-parameter $\psi$-estimator and finding an analogous lemma to our Lemma 4, then 2) using mathematical programming duality to simplify the Chernoff bound analysis, are a framework generalizable for tightly analyzing other estimators.

For simplicity of exposition, assume that we restrict the support of $D$ to some sufficiently fine-grained finite set, meaning that the maximization in Equation 5 is now finite-dimensional, albeit an arbitrarily large finite number. For each support element $x$, let $D_{x}$ be a variable representing the probability of choosing $x$ under distribution $D$. The expectation component of Equation 5 may now be expressed as sum that is a linear function in the variables $D_{x}$:

Using the standard max-min inequality (a form of weak duality in optimization), we have that Equation 6 is upper bounded by swapping the maximization and minimization (Equation 7), meaning that the vector $\boldsymbol{\mathbf{d}}$ no longer depends on the distribution $D$.

Crucially, however, Equation 7 is not just an upper bound on Equation 6, but is in fact equal to it, due to Sion’s minimax theorem [23]. To apply Sion’s minimax theorem, it suffices to check that 1) both $\boldsymbol{\mathbf{d}}$ and $D$ are constrained to be in convex sets, at least one of which is compact, 2) the objective is convex in $\boldsymbol{\mathbf{d}}$ and 3) concave in the variables $D_{x}$. For the first condition, we note that the set of distributions on a finite domain is compact. The objective is convex in $\boldsymbol{\mathbf{d}}$ since the objective is the sum of exponentials that are each linear in $\boldsymbol{\mathbf{d}}$. And the objective is concave in $D_{x}$ because it is in fact a linear function of $D$.

The guarantee of Sion’s minimax theorem means that we may work with Equation 7 instead of Equation 6 without sacrificing tightness in our analysis. This justifies why we are free to choose $\boldsymbol{\mathbf{d}}=(d_{\mu},d_{\alpha})$ in Lemma 4 that does not depend on the distribution $D$.

To further simplify the problem in Equation 7, we note again that both the objective and the constraints on $D$ are linear in the variables $D_{x}$, meaning that the inner maximization is in fact a linear program. We can then apply linear programming (strong) duality to yield the following equivalent optimization (Equation 8). We note that, as above, for the purposes of upper bounding Equation 5, it suffices to only use weak duality. Strong duality however guarantees that this step does not introduce slack into the analysis.

The three variables $V,M,S$ in the inner minimization below are the dual variables corresponding to the three constraints on distribution $D$ originally: that $D$ has variance 1, mean 0, and total probability mass 1.

We have thus reduced the infinite-dimensional optimization problem of Equation 5 to the five-dimensional problem of Equation 8 (or six dimensions, if we include the universal quantification for $x\in\mathbb{R}$), a significant simplification. We bound Equation 8 by explicitly choosing values for $\boldsymbol{\mathbf{d}}=(d_{\mu},d_{\alpha}),V,M,S$ as functions of $\hat{\alpha},n,\log\frac{1}{\delta}$, and showing that they jointly satisfy the constraint of Equation 8, for all $x$. We factor out the terms in the exponential that do not depend on $x$; we make the variable substitutions $y\equiv\sqrt{\hat{\alpha}}x$ and $\hat{v}\equiv\frac{\log(1/\delta)}{3n\hat{\alpha}}$ to replace dependence on $\hat{\alpha},n,\log\frac{1}{\delta}$ with dependence on the single variable $\hat{v}$; taking the log of both sides (and swapping sides) yields an expression that is recognizable in the following lemma, where the multipliers of $1,y,y^{2}$ respectively on the right hand side are essentially our choices of $S,M,V$: