Efficient median of means estimator

Let $X$ be a random variable with mean $\mu$ and variance $\sigma^{2}$. A sub-Gaussian estimator of $\mu$ based on a sample $\mathcal{X}=\{X_{1},\ldots,X_{N}\}$ of i.i.d. copies of $X$ is a measurable function $\widetilde{\mu}:=\widetilde{\mu}(\mathcal{X};t)$ such that $\mathbb{P}{\left(\left|\widetilde{\mu}-\mu\right|\geq C\sigma\sqrt{\frac{t}{N}}\right)}\leq ce^{-t}$ for a absolute constants $c,C>0$ and all $t\in\left[1,t_{\mathrm{max}}(N)\right]$. It is known (for instance, see the work by Catoni, (2012)) that $C\geq\sqrt{2}$. A natural question, posed previously by Devroye et al., (2016), is whether sub-Gaussian estimators with $C=\sqrt{2}+o(1)$, where $o(1)$ is a function that goes to $0$ as $N$ (and possibly $t$) tend to infinity, exist.

Several authors showed that such estimators can indeed be constructed under various additional assumptions. In one of the earliest works on the topic, Catoni, (2012) presented the first known example of sharp sub-Gaussian estimators for distributions with finite fourth moment and a known upper bound on the kurtosis, as well as for distributions with known variance. Construction by Devroye et al., (2016) similarly required the fourth moment to be finite. One of the strongest results is the one by Lee and Valiant, (2020): their estimator attains required guarantees uniformly over the class of distributions with finite variance, assuming just the finite second moment, albeit with $C=\sqrt{2}$ only in the limit as $t\to\infty$. Minsker, (2023) proposed a permutation-invariant version of the well known median of means (MOM) estimator (Nemirovski and Yudin,, 1983; Jerrum et al.,, 1986; Alon et al.,, 1996) and proved that it achieves desired guarantees for the class of distributions with more than $\frac{3+\sqrt{5}}{2}$ finite moments and “sufficiently regular” probability density functions.

The main goal of this essay is to present a modification of the “permutation-invariant” MOM estimator that attains sub-Gaussian guarantees with asymptotically optimal constants for distributions possessing $2+\varepsilon$ moments for some $\varepsilon>0$. This result could yield improvements for a variety of robust algorithms (e.g., see the survey by Lugosi and Mendelson, (2019)) that rely on the classical MOM estimator serves as a subroutine.

Let us recall the definition of the classical median of means estimator. Given an i.i.d. sample $\mathcal{X}=\{X_{1},\ldots,X_{N}\}$ from distribution $P$ with mean $\mu$ and variance $\sigma^{2}$, let $G_{1}\cup\ldots\cup G_{k}\subseteq[N]$ be a collection of $k$ disjoint subsets of cardinality $\lfloor N/k\rfloor$ each, $\bar{X}_{j}:=\frac{1}{|G_{j}|}\sum_{i\in G_{j}}X_{i}$ and $\widehat{\mu}_{\mathrm{MOM}}=\mbox{med}\left(\bar{X}_{1},\ldots,\bar{X}_{k}\right),$ where $\mbox{med}\left(\cdot\right)$ stands for the “median.” It is known that $\widehat{\mu}_{\mathrm{MOM}}$ satisfies the inequality $\mathbb{P}{\left(\left|\widehat{\mu}_{\mathrm{MOM}}-\mu\right|\geq C\sigma\sqrt{\frac{t}{N}}\right)}\leq 2e^{-t}$ with $C=\sqrt{\pi}+o(1)$, where $o(1)$ goes to $0$ as $k,N/k\to\infty$. Minsker, (2023) proved that allowing the overlapping subsets of data improves the constant: given $J\subseteq[N]$ of cardinality $|J|=\lfloor N/k\rfloor$, set $\bar{X}_{J}:=\frac{1}{|J|}\sum_{j\in J}X_{j}$ and define $\widehat{\mu}_{U}=\mbox{med}\left(\bar{X}_{J},\ |J|=\lfloor N/k\rfloor\right)$, where $\left\{\bar{X}_{J},\ |J|=\lfloor N/k\rfloor\right\}$ denotes the set of sample averages computed over all possible subsets of $[N]$ of cardinality $\lfloor N/k\rfloor$. Then $\widehat{\mu}_{U}$ attains sub-Gaussian deviations with $C=\sqrt{2}+o(1)$ under the assumptions described in section 1. Essentially, $\widehat{\mu}_{U}$ is a function of the order statistics which are complete and sufficient for the family of all distributions with finite variance.

Our construction, presented below, shows that it is not necessary to use all possible sample means, and that a much smaller collection of averages suffices: not only this makes computation easier, but the theoretical guarantees for the resulting estimator hold under weaker assumptions. The main idea is to split the data into subsets of size smaller than $\lfloor N/k\rfloor$, and construct all possible sample means using these subsets as “building blocks”. The size of the overlap is then naturally proportional to the size of the block. For instance, the estimator $\widehat{\mu}_{U}$ corresponds to the blocks of size $1$, resulting in the sample means over all possible subsets of a given size. Our results show that allowing the block size to be slowly growing with the the sample size could be beneficial. Formally, let $k,l$ be positive integers such that $\frac{\lfloor N/k\rfloor}{l}\in\mathbb{N}$. Assume that $G_{1}\cup\ldots\cup G_{lk}\subseteq[N]$ are disjoint subsets of cardinality $\lfloor\frac{N}{lk}\rfloor$ each, and $Z_{j}:=\bar{X}_{j}=\frac{1}{|G_{j}|}\sum_{i\in G_{j}}X_{i},\ j=1,\ldots,lk$. It will be convenient to set $n=lk$, $m=\lfloor\frac{N}{k}\rfloor$, and to view $Z_{1},\ldots,Z_{n}$ is a new i.i.d. sample; clearly, $Z_{1}$ has mean $\mu$ and variance $\frac{\sigma^{2}}{m/l}$. Given $J\subseteq[n]$ of cardinality $|J|=l$, set $\bar{Z}_{J}:=\frac{1}{l}\sum_{j\in J}Z_{j}$; note that $\bar{Z}_{J}$ is an average of $m$ observations from the original sample $\mathcal{X}_{N}$, same as in the definition of the standard MOM estimator. Define $\mathcal{A}_{n}^{(l)}=\left\{J\subset[n]:\ |J|=l\right\}$ and

where $\left\{\bar{X}_{J},\ J\in\mathcal{A}_{n}^{(l)}\right\}$ denotes the set of sample averages computed over all possible subsets of $[n]$ of cardinality $l$. In other words, $\widehat{\mu}_{N}$ is the median of means computed over overlapping subsets of data, where the size of the overlap is proportional to $\lfloor N/lk\rfloor$, the size of the block $G_{1}$. We remark here that all explicit, non-asymptotic deviations guarantees that are valid for the classical MOM estimator $\widehat{\mu}_{\mathrm{MOM}}$ automatically extend to $\widehat{\mu}_{N}$ in view of the so-called “Hoeffding representation” of U-statistics (Lee,, 2019) as the average of averages of independent random variables; pursuit of optimal constant however appears to require the bounds that include asymptotic terms. Everywhere below, it is assumed that $k,m,l$ and functions of the sample size $N$. We proceed with the statement of our main result. Denote

Feller, (1968) proved that $g(m)$ controls the rate of convergence in the central limit theorem, namely that $\sup_{t\in\mathbb{R}}\left|\Phi_{m}(t)-\Phi(t)\right|\leq g(m)$ where $\Phi_{m}$ and $\Phi$ are the distribution functions of $\frac{\sum_{j=1}^{m}X_{j}-\mu}{\sigma\sqrt{m}}$ and the standard normal law respectively. It is well known that $g(m)\to 0$ as $m\to\infty$ for distributions with finite variance. Moreover, $g(m)$ admits an upper bound of the form $g(m)\leq C\mathbb{E}\left|\frac{X_{1}-\mu}{\sigma}\right|^{2+\varepsilon}m^{-\varepsilon/2}$ whenever $\mathbb{E}|X_{1}-\mu|^{2+\varepsilon}<\infty$ for some $\varepsilon\in(0,1]$. In the context of the median of mean estimation, the role of $g(m)$ is to control the difference between the mean and the median corresponding to the distribution of $\frac{1}{m}\sum_{j=1}^{m}X_{j}$, which can be seen as the main contribution to the the bias of $\widehat{\mu}_{N}$.

Note that we can rewrite $h^{(1)}(Z_{1})$ as

Given an integer $r\geq 1$, let $\widetilde{\Phi}_{r}(t)$ be the cumulative distribution function of $\sum_{j=1}^{r}\tilde{X}_{j}$. Then

with $P_{Z_{1}\sqrt{m/l}}$ being the law of $Z_{1}\sqrt{m/l}$. Feller’s version of Berry-Esseen theorem implies that

where $\Phi$ is the distribution function of standard normal law. Therefore,

by assumption. At the same time,

where $q(x):=\left(\frac{m}{\sqrt{m-m/l}}\sqrt{\frac{t}{N}}-x\sqrt{\frac{m/l}{m-m/l}}\right)^{2}$ is such that $q(x)\to 0$ as $l\to\infty$ and $C(x,Z_{1})$ is a bounded function. Therefore, $h^{(1)}(Z_{1})-\sqrt{\frac{2}{\pi}}Z_{1}\sqrt{m/l}\to 0$ almost surely, assuming that $\sqrt{l}g(m)=o(1)$. Finally, note that

where the last inequality follows from the well known bound for the concentration function (Theorem 2.20 in the book by Petrov, (1995)); here, $C=C(P)>0$ is a constant that may depend on the distribution of $X_{1}$. We therefore conclude that the sequence $\left(h^{(1)}(Z_{1})-\sqrt{\frac{2}{\pi}}Z_{1}\sqrt{m/l}\right)^{2}$ is uniformly integrable (as $Z_{1}\sqrt{m/l}$ is), hence the claim follows.

The union bound together with Hoeffding’s decomposition entails that for any $t>0$ and $0<\varepsilon<1$ (to be chosen later as a decreasing sequence $\varepsilon(l)$),

where $h^{(j)},\ 2=1,\ldots,l$ are the degenerate kernels corresponding to the higher-order terms of Hoeffding’s decomposition (not to be confused with the derivatives!). Specifically,

where $\delta_{y}$ is the point measure concentrated at $y$; in particular, $\delta_{y}(h)=h(y)$. It is known that $h^{(j)}$ can be viewed geometrically as orthogonal projections of $h$ onto a particular subspace of $L_{2}(P_{Y}^{m})$. We refer the reader to the book by Lee, (2019) for futher details related to the background material on U-statistics and the Hoeffding’s decomposition. Bernstein’s inequality yields that

where $o(1)\to 0$ as $n/l\to\infty$ uniformly over $t$. It remains to control the expression involving higher order Hoeffding decomposition terms. To this end, we will show that it is bounded from above by $\operatorname{exp}\left(-\frac{t}{2\,\mbox{Var}\left(\sqrt{l}\,h^{(1)}(X_{1})\right)}\right)\cdot o(1)$ where $o(1)\to 0$ uniformly over the range of $t$. To this end, we will need concentration inequality for the U-statistics of growing order established in Minsker, (2023, Theorem 4.1). Set $t_{j,\varepsilon}=\left(\varepsilon\frac{\sqrt{t}}{j^{2}}\left(\frac{n}{l}\right)^{\frac{j-1}{2}}\right)^{2}$, and note that, in view of the union bound,

where the last inequality follows from the first bound of Theorem 4.1 in Minsker, (2023). Whenever $l\log^{2}(l)\ll k$ and $\varepsilon\gg\log^{-1/2}(l)$, the last expression is at most

In turn, it is bounded by $e^{-\frac{c_{2}t}{\varepsilon}}$ whenever $t<\frac{n}{l^{2}}\varepsilon^{4}$. Desired conclusion follows.

The author is grateful to the anonymous Referees for their insightful comments and constructive feedback that helped improve the quality of the presentation.