Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
80,595 characters · 14 sections · 52 citation commands
Tight Non-asymptotic Inference via Sub-Gaussian Intrinsic Moment Norm
Establishing rigorous error bounds and related finte--sample inferences for specific learning procedures is a popular topic in machine learning, statistics and econometrics communities, such as wainwright2019high,zhang2020concentration,yang2020non, romano2000finite,chassang2009non,arlot2010some,yang2020non,horowitz2020inference,armstrong2021finite,zheng2021finite,lucas2008rigorous,owhadi2013optimal,wang2020sub. In particular, concentration-based inferences have been developed when data is unbounded romano2000finite,auer2002finite,maurer2009empirical,hao2019bootstrapping,wang2021two,shiu2021efficient or Gaussian arlot2010some,duy2021exact,bettache2021fast,feng2021nonasymptotic. Notably, phan2021towards,austern2022efficient derived a result that is more refined than Hoeffding's inequality for bounded data, while bettache2021fast provided a concentration of the correlations for multivariate Gaussian data.
However, in reality, it is hard to learn the boundedness or distribution of data. If we mistakenly apply Hoeffding's inequality hoeffding1963probability to unbounded data, we will obtain very in-consistent confidence intervals (CIs); see Appendix SM1.1. Compromisely, it is a common practice to assume that data follows a sub-Gaussian distribution kahane1960. In this case, the Cramér–Chernoff method\footnote{For simplicity, we focus on zero-mean random variables (r.v.) throughout this paper when referencing all sub-Gaussian r.v.s.} boucheron2013concentration gives \[ \mathrm{P} (X \geq t) \leq \inf_{s > 0} [ \exp\{-st\} \mathrm{E} \{\exp (s X)\}] \] for all $t \ge 0$. Therefore, a tight control of the upper bound for the moment generating function (MGF) \(\mathrm{E} \{\exp (s X)\}\) is crucial in statistical inferences. This consideration leads to the concept of optimal variance proxy chow1966some,buldygin2000metric of sub-Gaussian distribution; see Appendix for more introductions of sub-Gaussianity.
The optimal proxy variance has a conservative lower bound: $\sigma^2_{opt}(X) \ge \operatorname{var}(X)$, as detailed in the Appendix SM1.2. When \( \sigma_{opt}^2(X) = \operatorname{var}(X) \), it denotes strict sub-Gaussianity arbel2020strict. Drawing from Theorems 1.5 in buldygin2000metric, one has $\mathrm{P} \left( X \ge t \right) \leq e^{{t^2}/\{{2 \sigma^2_{opt}(X)}\}}$ and
for independent sub-G r.v.s \( \{ X_{i}\}_{i = 1}^n \). By the Cramér–Chernoff method, ((ref)) offers the tightest upper bound in the form of \( \mathrm{P} (X > t) \leq e^{-C t^2}\) [or \( \mathrm{P} (|\sum_{i = 1}^n X_i | > t) \leq e^{-C t^2} \)], where $c$ a constant.
Given the unbounded data \( \{X_i\}_{i = 1}^n \stackrel{\rm{i.i.d.}} \sim \operatorname{subG}(\sigma^2_{opt}(X)) \), a direct application of ((ref)) yields a non-asymptotic \( 100(1-\alpha)\% \) CI:
A plug-in estimate\footnote{We note an inconsistent estimator, \(2 \inf_{t \in \mathbb{R}} \log [n^{-1}\sum_{i = 1}^n e^{t X_i}]/t^2\), from the statistical physics literature wang2020sub. It's important to stress that \(\widehat\sigma_{opt}^2(X)\) should take the supremum over \( {t \in \mathbb{R}} \), not the infimum.} for \( \sigma_{opt}^2(X) \) is given by $$\widehat\sigma_{opt}^2(X) := 2 \sup\limits_{t \in \mathbb{R}}{t^{-2}}\log M_n(t),$$ where $M_n(t):=n^{-1}\sum\nolimits_{i = 1}^n \exp(t X_i)$ is the empirical MGF. However, $\widehat\sigma_{opt}^2(X)$ is not practically useful for at least two reasons: (i) optimization outcome is unstable due to potential non-convexity of the objective function; (ii) variance of logarithm empirical MGF explodes as $t$ diverges, i.e., \[ \mathrm{var}\left\{\frac{\log (M_n(t))}{t^2}\right\}\approx \frac{\exp(2t\mu){\mathrm{var}\{X_1\}}}{nt^2M^2(t)}\to \infty\;\text{as}~t \to \infty,~\mu \ne 0; \] see Appendix SM1.3 for the proof. These concerns are empirically supported by the subsequence simulations in Sections (ref) & (ref), e.g., Figure (ref).
On the other hand, several other forms of variance-type parameters exist in the literature. For example, vaart1996weak introduced the Orlicz norm, defined as \(\|X\|_{w_2}:=\inf\{c>0: \mathrm{E} \{\exp (|X|^2/c^2) \}\leq 2\}\); vershynin2010introduction proposed a norm based on the scaling of moments, expressed as \(\|X\|_{\psi_2} := \max_{k \geq 2} k^{- 1/ 2} \mathrm{E} \{( |X|^k )^{1 / k}\}\). However, as illustrated in Table (ref) and further discussed in Appendix SM1.4, these two norms often fall short in providing {nice} probability bounds even for strict sub-Gaussian distributions, such as the standard Gaussian and symmetric Beta distributions.
Recently, horowitz2020inference derived confidence intervals for structural parameters in econometrics using a class of optimization problems based on Berry-Esseen (B-E) bounds. When compared to the traditional normal approximation underpinned by B-E bounds, our results offer notable advantages, especially for datasets with extremely small sample sizes. For instance, in certain experimental science measurements, the sample size ranges merely from \( n=4 \) to \( n=8 \) rousseeuw2002robust. We showcase this by examining Bernoulli observations, comparing two types of confidence intervals: one based on the B-E-corrected Central Limit Theorem (CLT) and the other on Hoeffding's inequality, as illustrated in Figure (ref).
In light of previous discussions, we advocate the use of the intrinsic moment norm (see Definition (ref)) in the construction of tight non-asymptotic CIs for two reasons: (i) it approximately recovers tight inequalities (ref); (ii) it can be easily estimated in a robust manner.
Note that \(\| X\|_G < \infty\) implies \(\sigma_{opt}(X) < \infty\) based on Theorem 2.6 in wainwright2019high. Thus, a random variable with a finite intrinsic moment norm is ineeded sub-Gaussian.
In this paper, our contributions are summarized as
This article is organized as follows: In Section 2, we present upper and lower bounds for the tail probabilities of sub-Gaussian variables and introduce the sub-Gaussian plot to determine whether data follows a sub-Gaussian distribution. In Section 3, we delve into the finite sample properties of the intrinsic moment norm, exploring its inherent characteristics and the concentration for summations' sharper bounds. In Section 4, we propose estimation techniques for this norm, including the plug-in estimator, the robust mean-of-median estimation, and methods for extremely small sample sizes, along with a discussion of their theoretical properties. In Section 5, we apply the insights derived from this norm to the non-asymptotic analysis of problems found in the multi-armed bandit. Detailed proofs and supplementary information can be found in the Appendix.
Before estimating $\|X\|_G$, the first step is to verify the data $X$ is indeed sub-G given its i.i.d. copies $\{X_i\}_{i=1}^n$. For bounded $X$, it is naturally sub-Gaussian by using Hoeffding's inequality. For unbounded $X$, Corollary 7.2 (b) in zhang2020concentration shows for r.v.s $X_i \sim \operatorname{subG}(\sigma_{opt}^2(X))$ (without independence assumption)
which implies $\max_{1 \leq i \leq j} X_i =O_{\mathrm{P}}(\sqrt{\log j})$. Moreover, we will show the rate $\sqrt{\log j}$ is indeed sharp for a class of unbound sub-G r.v.s characterized by the lower intrinsic moment norm below.
By the method in Theorem 1 of zhang2020non, we obtain the following tight rate result with a lower bound.
For instance, a mixture of two independent Gaussian variables satisfies the condition ${\left\| X \right\|_{\tilde G}}>0$ in Theorem (ref)(a). Refer to Example SM1.5 in Appendix for verification. For a Gaussian variable $X$, ${\left\| X \right\|_{G}}={\left\| X \right\|_{\tilde G}}$ and the constant $L(X)$ takes on the value of $1 / 2$ in the aforementioned lower bound. The upper bound follows from the proof of (ref) similarly. The proof of lower bound relies on the following sharp reverse Chernoff inequality, which is derived from Paley--Zygmund inequality (see paley1932note).
Theorem 1 of zhang2020non does not optimize the constant in Paley-Zygmund inequality. In contrast, our Lemma (ref) has an optimal constant; see Section (ref).
Motivated by Theorem (ref), we propose a novel sub-Gaussian plot to check whether i.i.d data $\{ X_i\}_{i = 1}^n$ follow a sub-G distribution.
Sub-Gaussian plot under unbounded assumption. Suppose that for each $j$, $\{ X_i^*\}_{i = 1}^j$ are independently sampled from $\mathbb{F}_n(x) = \frac{1}{n} \sum_{i = 1}^n \mathds{1}(X_i \leq x)$ as the empirical distribution of $\{ X_i\}_{i = 1}^n$. Specifically, we plot the order statistics $\{ \max_{1 \leq i \leq j} X_i^*\}_{j = 1}^n$ on the plane coordinate axis, where $x$ axis represents $\sqrt{\log j + 1}$ and $y$ axis the value of $ \max_{1 \leq i \leq j} X_i^*$. We check whether those points have a linear tendency at the boundary: the closer they are to the tendency of a straight line, the more we can trust the data are sub-Gaussian.
The Figure (ref) illustrates sub-Gaussian plot of $N(0, 1)$ and $\operatorname{Exp}(1)$. It can be seen that sub-Gaussian plot of $N(0, 1)$ shows a linear behavior near the boundary, while $\operatorname{Exp}(1)$ shows quadratic tendency at the boundary. For the quadratic tendency, we note that if $\{ X_i\}_{i = 1}^n$ have heavier tails such as sub-exponentiality, then $\max_{1 \leq i \leq j} X_i =O_{\mathrm{P}}( {\log j})$ instead of the order $O_{\mathrm{P}}(\sqrt{\log j})$; see Corollary 7.3 in zhang2020concentration.
Sub-Gaussian plots provide an insightful way to visualize sub-Gaussian distributions, serving a similar role to how Q-Q plots represent Gaussian distributions. While both types of plots exhibit linear trends, they achieve this in different ways: Q-Q plots rely on the actual distribution of the data, whereas sub-G plots make use of the MGF bound. To the best of the author's knowledge, this paper represents the first formal attempt to introduce a methodology for determining whether a dataset follows a sub-Gaussian distribution. It is worth emphasizing that the utility of sub-G plots is most pronounced for datasets with a large enough sample size. For small values of $n$, the data is approximately treated as bounded random variables, and there is no need to use a sub-G plot.
Another consideration is the data that exhibits heavy-tailed outliers, the next section will discuss employing mean-of-median estimation. This approach allows for the exclusion of outliers, ensuring that the inliers conform to a sub-Gaussian distribution. This will particularly suitable if the inliers are approximately bounded, a scenario often encountered in real-world datasets.
In this section, we characterize two important properties of the intrinsic moment norm that are used in constructing non-asymptotic confidence intervals.
Lemma (ref) below establishes that the intrinsic moment norm is estimable.
Lemma (ref) ensure that for any sub-G variable $X$, its intrinsic moment norm can be computed as \[\|X \|_G := \max_{m \in 2 \mathbb{N}} \left[ \frac{\mathrm{E} (X^{m})}{{(m - 1)!!}} \right]^{1 / m} = \max_{1 \leq k \leq k_{X}} \left[ \frac{\mathrm{E} (X^{2k})}{{(2k - 1)!!}} \right]^{1 / (2k)}\] with some finite $k_X < \infty$. This is an important property that other norms may not have. The $\|X\|_{\psi_2} := \max_{k \geq 2} k^{- 1/ 2} ( \mathrm{E} |X|^k )^{1 / k}$ for Gaussian $X$ achieves its optimal point at $k = \infty$; see Example (ref) in Appendix (ref).
As for $\sigma_{opt}^2(X):=2 \sup_{t \in \mathbb{R}} \frac{ \log \{\mathrm{E} [\exp(t X)]\}}{t^2}$, it is unclear that its value can be achieved at a finite $t$. Note that if $k_X = 1$, one has $\| X\|_G^2 =\operatorname{var}(X)$.
Next, we present an example in calculating the values of $k_X$. Denote $\left. \operatorname{Exp}(1) \right|_{[0, M]}$ as the truncated standard exponential r.v. on $[0, M]$ with the density as $f(x) = \frac{e^{-x}}{\int_0^M e^{-x} \, dx} 1_{\{x \in [ 0, M]\}}(x)$.
{In the following discussion, we demonstrate an additional property of $\|X \|_G$, illustrating that it nearly recovers the tight MGF bounds as defined in Definition (ref). More significantly, this property allows us to derive the sub-Gaussian version of Hoeffding's inequality, as shown in (ref).}
Theorem (ref)(a) can be found in Theorem 2.6 of wainwright2019high. As for Theorem (ref)(b), it's noteworthy that \((17/12)^{1/2}\) is approximately 1.19. In contrast, Lemma 1.5 in buldygin2000metric provides the bound $\mathrm{E} \{\exp\{tX_i\}\}\le \exp\big\{ {{\frac{t^2}{2}({3.1}^{1/4}{{\| X_i\|}_{G}})^2}}\big\}$ for \(t \in \mathbb{R}\), where \({3.1}^{1/4}\) is approximately 1.32. {Essentially, $(17/12)^{1/2}>1$ appears for asymmetric variables, since ${\| \cdot\|}_{G}$ is defined by comparing a Gaussian variable $G$ that is symmetric}. A technical improvement in our approach is due to the fact that \(\| \cdot\|_G\) bypasses the need for Stirling's approximation, thereby providing a more precise MGF bound when expanding the exponential function using Taylor's expansion. To further demonstrate the precision of Theorem (ref)(b), Figure (ref) offers a comparative analysis with various parameters such as \(\sigma_{o p t}(X)\), \((17/12)^{1/2}\| X \|_{G}\), \((10e)^{1/2}\|X\|_{\psi_2}\), \((5/2)^{1/2}\|X\|_{w_2}\), and \((\operatorname{var}X)^{1/2}\). The comparisons specifically explore the confidence lengths {under asymmetric Bernoulli and Beta distributions,} as detailed in Table (ref).
Ideally, suppose we have i.i.d. mean-zero sample generating from $X$. Then a straightforward approach to estimate $\| X\|_G$ is by the plug-in approach. Although $k_X$ is proven to be finite in Lemma (ref), its (possibly large) exact value is still unknown in practice. Hence, we suggest a feasible plug-in estimator with a non-decreasing index sequence $\{\kappa_n\}$ to replace $k_X$ {as follows:}
For example, $\kappa_n=\lceil \log n \rceil$ or a large constant.
Deriving the non-asymptotic property of the $\widehat{\|X\|}_G$ is not an easy task: the maximum point \[\hat k(\kappa_n):=\arg \max_{1 \leq k \leq \kappa_{n}} \left[ {\frac{{1}}{{(2k - 1)!!}} \frac{1}{n} \sum\limits_{i = 1}^n {X_i^{2k}}}\right]^{1/(2k)}\] will change with the sample size $n$ even $\kappa_n$ is fixed. To resolve this, we first examine the oracle estimator defined as
Here, based on Orlicz norm $\|Y\|_{\psi_{\theta}}:=\inf \{t>0: \mathrm{E} \exp (|Y|^{\theta} / t^{\theta}) \leq 2 \}$ of sub-Weibull r.v. $Y$ with $\theta>0$ hao2019bootstrapping, we present the non-asymptotic concentration of $\widetilde{\| X\|}_G$ {around it ture value $\|X \|_G$.}
In practice, datasets often contain a significant number of outliers; see hsu2016loss. The exponential-moment condition $\| X\|_{\psi_{1/k_X}}<\infty$ is too strong for the error bound of $ \widetilde{\|X \|}_G^{2 k_X} - \| X\|_G^{2 k_X} $ in Proposition (ref). {Additionally, the constants in (ref) may potentially be large, leading to overly loose CIs when the sample size $n$ is small. }
Apart from the direct plug-in estimator, we utilize the median-of-means (MOM) as described on Page 244 of nemirovskij1983problem as a robust alternative. Given a positive integer \( m \) and \( b \) such that \( n=mb \), let \( B_{1}, \ldots, B_{b} \) partition the set \( [n] \) into blocks of equal size \( m \). For every \( s \in[b] \), and define \( \mathbb{P}_m^{B_{s}} X \) as \( m^{-1} \sum_{i \in B_{s}} X_{i} \) for the independent dataset \( \{X_{i}\}_{i=1}^n \). We can then express the MOM-based intrinsic moment norm estimator as:
The naive plug-in estimator in Section (ref) can be just reviewed as a special case of MOM estimator with $b=1$. As stated in Proposition (ref), the naive plug-in estimator $\widehat{\|X \|}_G = \widehat{\| X \|}_{1,G}$ is not robust. MOM estimators (ref) with $b\gg 1$ have two merits: (a) it only needs finite moment conditions, but the exponential concentration bounds are still achieved; (b) it permits some outliers in the data. For non-asymptotic inference, we need to precisely bound \( \| X\|_{G} \) using a feasible estimator \( {\widehat{\| X \|}_{b,G}} \) while ensuring sharp constants. We proceed to establish a high-probability upper bound for the estimated norm under the following outlier assumptions labeled as \( O \cup I \):
$\bullet$ (M.1) Suppose that the dataset $\{X_{i}\}{i=1}^n$ consists of $n-n_o$ inliers $\{X_{i}\}{i \in I}$, each i.i.d. following a target sub-Gaussian distribution (for the $n_o$ outliers ${X{i}}_{i \in O}$, no distributional assumptions are made).
$\bullet$ (M.2) $b=b_O+b_S$, where $b_O$ is the number of blocks containing at least one outliers and $b_S$ is the number of sane blocks containing no outliers. Let $\varepsilon:=n_o/n$ be the fraction of the outliers and ${n_o}/{b}<{1}/{2}$. Assume here exists a fraction function $\eta(\varepsilon)$ for sane block such that $b_S \ge \eta(\varepsilon) b~\text{for a function}~ \eta(\varepsilon)\in (0,1].$
(M.1) is a reasonable assumption in our sub-Gaussian inference, since when the outliers are thrown away with a large proportion such that ${n_o}/{b}<{1}/{2}$ that are allowed in (M.2) the remaining inliers is approximately bounded data that are naturely sub-Gaussian. To serve for error bounds in the presence of outliers, (M.2) considers the specific fraction function of the polluted inputs; see laforgue2021generalization.
Define ${{\underline{g}_{k,m}}({\sigma _k})}$ and ${\bar g_{k,m}({\sigma _k})}$ as the sequences for any $m \in \mathbb{N}$ {and} $1 \le k \le \kappa_{n}$: $ {\overline g_{k,m}({\sigma _k})} : = 1 - \left[ \frac{{{\mathrm{E}}(X^{2k})}}{{(2k - 1)!!}} \right]^{-1/(2k)}\mathop {\max }\limits_{1 \le k \le \kappa_{n}} {\left[ {\frac{{{\mathrm{E}}(X^{2k})}}{{(2k - 1)!!}} -2\{\frac{m}{\eta(\varepsilon)}\}^{-1/2}\frac{\sigma_{k}^{k}}{{\mathrm{E}} (X^{2 k})} } \right]^{1/(2k)}}$ and ${ \underline{g}_{k,m}({\sigma _k})}:=[{2\{m/\eta(\varepsilon)\}^{-1/2}\sigma_{k}^{k}}/{{\mathrm{E}} (X^{2 k})}+1]^{1 /(2 k)}-1.$ A robust and non-asymptotic CI for $\left\| X \right\|_{G}$ is obtained as follow.
Theorem (ref) ensures the concentration of the estimator $\widehat{\| X\|}_{b,G}$ when $\kappa_{n} \ge k_X$ under enough sample. {If $\eta(\varepsilon)=1$ with $\varepsilon=0$, then the data are i.i.d., which have no outlier, and outlier assumptions in M.1-M.2 can be dropped in Theorem (ref)}. {When the data is i.i.d. Gaussianc variable, Proposition 4.1 in auer2002finite also gave a high-probability estimated upper bound for Gaussian standard deviations; our result is for the sub-G intrinsic moment norm and the Gaussian standard deviation is a special case of the sub-G norm.}
In practice, the block number $b$ can be taken by the adaptation method based on the Lepski method depersin2022robust. However, they only discuss Lepski's method in a theoretical context, and no simulation studies were conducted due to the method involving a large search radius with a constant of $808\times 1200$ (see Page 514). In practical terms, implementing Lepski's method can be computationally intensive, given the use of a dyadic grid for $b$. As an alternative, we propose the leave-one-out cross-validation (LOCV) method to determine $b$:
The high probability guarantee in Theorem (ref) requires that the index sequence $\kappa_{n}$ should not be very large for fixed $b$. The larger $\kappa_{n}$ needs larger $b$ in blocks $B_{1}, \ldots, B_{b}$. In the simulation, we will see that an increasing index sequence $\kappa_{n}$ with a slow rate will lead to a good performance.
We compare the estimated optimal variance proxy from lieber2022estimating with our estimated intrinsic moment norm. The latter is calculated using the direct method (equation (ref)), the MOM estimator (equation (ref)) with a fixed block number $b \equiv 5$, and the MOM estimator with the block number determined by equation (ref). We consider variables $X$ that follow standard Gaussian and Rademacher distributions. In both cases, we have $\| X \|_G^2 = \sigma_{opt}^2 (X) = \operatorname{var}(X) = 1$. Figure (ref) depicts the performance of the three estimators for sample sizes ranging from $n=10$ to $1000$, with $\kappa_n=\lceil \log n \rceil$. From Figure (ref), we observe that both variants of the MOM estimator (fixed block size and LOCV-determined block size) outperform others in terms of accuracy. In comparison, the plug-in estimator from lieber2022estimating shows a less consistent performance.
For high-quality datasets with a very limited sample size \( n \), specifically when \( n \) is extremely small, the leave-one-out Hodges-Lehmann method rousseeuw2002robust can be employed to effectively augment the sample size. In scenarios where the sample size is small (\( n \leq 20 \)), apart from the direct empirical moment method (DE), there are two alternative techniques to consider.
The first method, a variant of the Bootstrap approach, employs the ($n-1$)-out-of-$n$ Bootstrap technique. In this procedure, \( n - 1 \) Bootstrap estimators are generated from the data by iteratively resampling \( n - 1 \) times while leaving one observation out during each iteration. Once these estimators are constructed, the median of these estimators is taken as the final estimate. This method increases the robustness of the estimator, particularly in the context of small sample sizes. The second approach, known as the Leave-One-Out Hodges-Lehmann method (LOO-HL), has been introduced by rousseeuw2002robust. To elucidate, given a sample \( \{X_i \}_{i = 1}^n \), the LOO-HL empirical mean estimator can be explicitly defined as: \[ \widehat \mu_{\rm LOO} := \operatorname{med} \left\{ \frac{\operatorname{DE} (X_{(-i)}) + \operatorname{DE} (X_{(-j)})}{2} : 1 \leq i < j \leq n \right\}, \] Here, the function \( \operatorname{DE} \) computes the empirical mean of the dataset excluding a particular observation. Mathematically, \( \operatorname{DE} (X_{(-i)}) :=\sum_{k \ne i}X_k/(n-1) \) is the empirical mean estimator when data point \( i \) is excluded. This mean is calculated over the set $X_{(-i)} = \{ X_1, \ldots, X_{i - 1}, X_{i + 1}, \ldots, X_n\}$. Both of these methods aim to enhance the reliability and robustness of estimates derived from datasets with limited observations. By focusing on sub-samples or manipulating the data in a particular way, they mitigate some of the vulnerabilities associated with small sample statistics.
The nuances of the leave-one-out Hodges-Lehmann method differ significantly from the classic Hodges-Lehmann empirical mean estimator. The classic Hodges-Lehmann empirical mean estimator is represented as $\widehat \mu_{\rm HL} :=\operatorname{med}_{i \leq j} \frac{X_i+X_j}{2}$. In essence, this method computes the median of the means of all possible pairs of the sample data. Contrastingly, the leave-one-out version employs the leave-one-out mean, \( \operatorname{DE} (X_{(-i)}) \), as opposed to directly using individual data points \( X_i \). The departure from the classic Hodges-Lehmann method is driven by empirical findings, wherein the traditional method did not offer satisfactory performance under certain conditions.
To challenge these methods further, the underlying data gets perturbed with a \( 1\% \) independent Cauchy$(0, 5)$ noise, introducing outliers and asymmetry. The intention behind this perturbation is to assess the robustness of each method under non-ideal circumstances. The performance evaluation is in Figure (ref): while all three methods manage to yield reasonably good performance, the LOO-HL method seems to be in the lead. Its superiority becomes more pronounced for very small sample sizes, particularly when \( 2 \leq n \leq 4 \). This observation underscores the importance of selecting appropriate statistical techniques, especially the sample is extremely small.
In the multi-armed bandit problem (MAB), a player must choose between \( K \) slot machines, termed an \( K \)-armed bandit. Each machine provides a unique, unknown random reward with a random variable set \( \{Y_{k}\}_{k=1}^K \subseteq \mathbb{R} \). Every realization from a specific arm \( k \) is independent and identically distributed. Furthermore, we assume these rewards to be sub-Gaussian:
The primary objective is to identify the optimal arm with the highest expected reward, denoted as \( Y_{k^{\star}} \), through arm pulls. In every round \( t \in [T] \), the player selects an arm \( A_{t} \in [K] \). Given \( \{A_{t}=k\} \), the observed reward is defined as \( \{Y_{k,t}\}_{t \in [T]}\stackrel{\rm{i.i.d.}} \sim P_k \). The MAB exploration's aim is to minimize the cumulative regret after \( T \) rounds, represented as: \[ \operatorname{Reg}_{T}(Y,A) := \sum_{t=1}^{T}(\mu_{k^{\star}}-\mu_{A_{t}}). \] Without loss of generality, let \( {k^{\star}}=1 \). We seek to evaluate the expected bounds from the decomposition (see Lemma 4.5 in lattimore2020bandit),
{where \( \Delta_{k}:=\mu_1-\mu_k \) is the gap that represents the sub-optimality for arm \( k \) relative to the best arm. Here,\( \mathrm{E} \) is calculated based on the randomness in player's actions \( \{A_{t}\}_{t \in [T]} \).}
For each round \( t \), let \( T_{k}(t) \) denote the number of pulls for arm \( k \) up to time \( t \). The running average of rewards from arm \( k \) at time \( t \) is defined as: $\overline{Y}_{T_{k}(t)} := \frac{1}{T_{k}(t)} \sum_{\tau \leq t, A_{\tau}=k} Y_{k,\tau}.$ Given a tight concentration inequality, a \( 100(1-\delta)\% \) CI for \( \mu_k \) can be extracted as:
This confidence interval allows us to establish an upper bound for the reward of arm \( k \), \( \overline{Y}_{T_{k}(t)}+{c}_{k}(t) \), and subsequently play arm \( A_{t}=k \) to maximize the reward with a high probability for a given finite \( t \). This strategy is known as the Upper Confidence Bound (UCB) method; see auer2002finite.
Many works based on UCB methods appears; see lattimore2020bandit. However, many existent bandit algorithms contain unknown norms for the random rewards, they are actually hard to use. For example, hao2019bootstrapping use bootstrap method with the second order correction to give an algorithm with the explicit regret bounds for sub-Gaussian rewards. {However, the UCB algorithm in hao2019bootstrapping necessitates the use of the unknown Orlicz-norm (or its upper bound) of $\overline{Y}_k - \mu_k$. Consequently, employing this algorithm in practice proves to be infeasible, or it becomes risky to use an arbitrary value for the unknown norm (or its upper bound).}
Fortunately, our estimator can solve this problem. Suppose that $Y_{k}-\mu_{k}$ is symmetric around zero, by one-side version of Theorem (ref), the (ref) implies that for all $k$ and all $t$, \[ \mathrm{P} \left(\overline{Y}_{T_{k}(t)} > \mu_{k}+\|Y_{k}-\mu_{k}\|_{G}\sqrt{\frac{2}{T_{k}(t)}\log \frac{1}{\delta}} \right)\le \delta. \] Let sub-sample size $m_k$ and block size $b_k$ be positive integers such that $T_k(t) = m_k b_k$ for MOM estimators $\widehat{\|Y_{k}-\mu_{k}\|}_{b_k,G}$ in Section 3. Theorem (ref) {\rm{(a)}} guarantee that true norms can be replaced by MOM-estimated norms such that
as $\eta(\varepsilon)=1$ with $\varepsilon=0$.
In real-world scenarios, underlying data distribution details are mostly unknown. The multiplier bootstrap is a versatile method for uncertainty qualification, enabling the simulation of the non-asymptotic properties of the target statistic by reweighting its summands. The multiplier bootstrapped quantile for the i.i.d. observation \( \mathbf{Y}_n \) is: \[ q_{\alpha} (\mathbf{Y}_n - \overline{Y}_n, \mathbf{w}) = \inf \left\{ x \in \mathbb{R} \mid \mathrm{P}_{w} \left( n^{-1} \sum_{i = 1}^n w_i (Y_i - \overline{Y}_n) > x \right) \leq \alpha \right\}, \] where \( \mathbf{w} \) are bootstrap random weights independent of \( \mathbf{Y}_n \).
{Inspired by hao2019bootstrapping, we developed Algorithm (ref) that utilizes MOM estimators for the UCB. Our innovative regret bounds are articulated in Theorem (ref). These bounds assure comparatively minimal regret by incorporating a concentration-based second-order correction, $\widehat{\varphi}_G (\mathbf{Y}_{T_k(t)})$, as elucidated in Theorem (ref), into the bootstrapped threshold $q_{\alpha / 2} (\mathbf{Y}_{T_k(t)} - \overline{Y}_{T_k(t)}, \mathbf{w})$.}
For $\eta(\varepsilon)\in (0,1)$ where $\varepsilon\in(0,1)$, the regret bounds are comparable up to certain constants.
According to Theorem (ref), the regret of our approach attains the minimax rate of $\log T$ in a problem-dependent scenario, and $\sqrt{TK}$ for a problem-independent one. This is consistent up to a logarithmic factor, corroborated by tao2021optimal, demonstrating the {near-optimal property} of Algorithm (ref). In comparison to the traditional vanilla UCB, our technique refines the constants. To illustrate, when $Y_k \sim N(\mu_k, 1)$, the regret bound's constant factor in Theorem 4 of auer2002finite is $256$, which is significantly greater than the $16(2+\sqrt{2})^2$ presented in our theorem.
In scenarios where the UCB possesses unknown sub-G parameters, Theorem (ref) primarily investigates a practical UCB algorithm integrating sub-G parameter estimations. {In contrast, many earlier UCB algorithms operate under the assumption that the upper bound of the sub-G parameter is known beforehand.} One can refer to the algorithm detailed in hao2019bootstrapping for an example.
Next, we present a simulation of Algorithm (ref) for two sub-G cases to assess the performance of the estimated norms. Taking cues from hao2019bootstrapping,wang2020residual, we design the four UCB methods:
1. Use our method $\widehat{\varphi}_1(\mathbf{Y}_{T_k(t)})=\widehat{\varphi}_G (\mathbf{Y}_{T_k(t)})$ with Estimated Norm in Theorem (ref) with block number $b_k = 2 \lceil \log T_k(t) \rceil$;
2. Let
where $\widehat\sigma_{opt}^2(Y_k)$ is estimated using the formula in Remark 24 of lieber2022estimating.
3. Use Asymptotic Naive varphi $\widehat{\varphi}(\mathbf{Y}_{T_k(t)})$ satisfying \[\mathrm{P} \big(|\overline{Y}_{T_k(t)} - \mu_k| \leq \widehat{\varphi}_3(\mathbf{Y}_{T_k(t)}) \big) \rightarrow \alpha\]by CLT, i.e.
with $\widehat{\sigma}_k =\{\frac{1}{T_k(t)} \sum_{\tau \leq t, A_{\tau}=k} (Y_{k,\tau} - \overline{Y}_{T_k(t)})^2\}^{1/2}$ as the estimated standard deviation;
4. Regard all the unbounded rewards as bounded r.v. and use Hoeffding's inequality (wrongly use Hoeffding's inequality) to construct $\varphi$, i.e. \[ \widehat{\varphi}_4 (\mathbf{Y}_{T_k(t)}) = \big[\max\{ \mathbf{Y}_{T_k(t)} \} - \min \{\mathbf{Y}_{T_k(t)}\}\big] \left\{\frac{\log (2 / \alpha)}{2T_k(t)}\right\}^{1/2}. \] where $\max\{ \mathbf{Y}_{T_k(t)} \}:=\max\{Y_1, Y_2, \ldots, Y_{T_k(t)} \}$.
For our MAB simulation, we consider as follows, in each case, the number of arms is assigned as $K = 10$, and the mean reward from $k$-th arm $\mu_k$ is independently drawn from $\operatorname{Exp}(1)$. We consider two types of unbounded reward distributions: EG1. Gaussian $N(\mu_k, 1)$; EG2. Mix of two independent Gaussian
with $p_k \, \sim \, U(0, \frac{1}{2})$. Simultaneously, we account for potential data contamination due to Cauchy noise. For any given arm, the rewards generated in each round have a 15% chance of being affected by a $\operatorname{Cauchy}(0,0.01)$ (standard Cauchy distribution with a scale parameter of 1%). Additionally, we utilize Thompson Sampling agrawal2017near with Gaussian distributions for both rewards and priors. We tune the prior parameters on $\{ 2^{4- k} \}_{k = 1}^6$, using its optimal performance as a robust baseline for this simulation.
Both EG1 and EG2 represent sub-G reward structures. In the simulation, $\mu_k$ may not be bounded, complicating this problem. The simulation results are shown in Figure (ref). These results show that our method outperforms the other two methods under unbounded sub-Gaussian rewards and is even comparable to Thompson Sampling when sufficient correct prior knowledge is available. Furthermore, our method demonstrates the strong robustness of our estimated norm method.
Proof of Theorem (ref): Denote $Y_{k,i}$ represents the $i$-th value of reward of arm $k$ in its history, where $i \in [T_k(t)]$ for any round $t$. It is needed to bound $\mathrm{E} [T_k(t)]$ as one has \[\mathrm{E} (\operatorname{Reg}_T) \leq \sum_{k = 2}^K \Delta_k \mathrm{E} [T_k(t)],\] from (ref), so we can only focus on some fixed arm. Hence, we can just drop the subscript $k$ in $\{ Y_{k, i} \}_{i = 1}^{T_k(t)}$ as $\{ Y_{i}\}_{i = 1}^{T_k(t)}$.
We first give a lemma, which is crucial in the following proof.
Based on Lemma (ref), next we can prove Theorem (ref). We first state the assumptions for Theorem (ref) in detail.
Denote the population version of
with $\mathrm{E}(Y) = \mu$. Theorem (ref) gives
by the fact that $Y - \mu$ is symmetric around zero.
From Theorem (ref), we take $b \geq 8 \log (k_{Y} / \alpha)$ and define the event $\mathcal{E}_Y$ bellow associated with the MOM estimation of the intrinsic moment norm
with probability at least $1 - \alpha / 2$. Note that
where the last inequality is by taking $m$ big enough such that $\max_{1 \leq k \leq k_{Y} / 2} \overline{g}_{k, m}(\sigma_k) \leq 1 / \sqrt{n}$. Then, for $s \in \mathbb{N}_+$,
Now for any $k \in [K]$ and fixed $T_k(t) = s$, we know that
By the non-asymptotic second-order correction [see Theorem 2.2 in hao2019bootstrapping] and the assumption that $Y - \mu$ is symmetric around zero, one has \[ \mathrm{P} \left\{ \mu_k - \overline{Y}_s \geq q_{\alpha / 2} (\mathbf{Y}_s - \overline{Y}_s) + \left\{\frac{2 \log (4 / \alpha)}{s}\right\}^{1/2} \widehat{\varphi}_G(\mathbf{Y}_s) \right\} \leq 2 \alpha, \] where $q_{\alpha / 2} (\mathbf{Y}_{B_k} - \overline{Y}_{B_k}) : = q_{\alpha / 2}\left(\mathbf{Y}_s -\overline{Y}_s, \frac{1}{s} \mathbf{1}_s \right)$.
Denote the UCB index $\operatorname{UCB}_k (t) = \overline{Y}_{T_k(t)} + h_{\alpha} (\mathbf{Y}_{T_k(t)})$, and the good event for $k \in [K]$
where $B_k \in [T]$ is a constant to be chosen later. Following from the proof in (B.16)-(B.18) of hao2019bootstrapping, we can gives that $T_k(t) \leq B_k$ and
On the other hand, from Lemma (ref), then
Next, we need to the following assumptions for block $\{ b_k\}_{k \in [K]}$ in MOM estimator corresponding to Theorem (ref). Here we include the subscript $k$ to avoid confusion.
Under $\eta(\varepsilon)=1$ with $\varepsilon=0$, Theorem (ref) ensures for $b_k \geq 8 \log (k_{Y} / \alpha)$ and take $m$ large enough such that $\max_{1 \leq \kappa \leq k_{Y}} \underline{g}_{\kappa, m} (\sigma_\kappa) \leq n^{-1 / 2}$, then
Hence, we have
which implies with probability at least $1 - \alpha$,
where $B_k \geq 8 \log ( k_{Y_k}/ \alpha) \vee 2 \log (4 / \alpha)$ and $\max_{1 \leq \kappa \leq k_{Y_k} / 2} \overline{g}_{\kappa, m} (\sigma_{k, \kappa})\leq 1 / \sqrt{{B_k}}$ for each arm $k$.
Now, define the event
with $\Delta_k := \mu_1 - \mu_k$. Choose $B_k$ as
we have
Applying Theorem (ref) for the concentration of $\overline{Y}_{B_k} - \mu_k$ when $B_k$ is chosen as in (ref),
Taking account these results into (ref), we get that
by taking $\alpha = 4 / T^2$. Under the problem-dependent case, the regret is bounded by
To get the problem-independent bound, we let $\Delta > 0$ as an arbitrary threshold, then decompose $\operatorname{Reg}_T$, we get
by taking $\Delta = 8 (2 + \sqrt{2}) C \{(K \log T)/ T\}^{1/2}$. And finally, we take $C = \max_{k \in [K]}\|Y_k - \mu_k \|_G$.
{H. Zhang's research received support from the National Natural Science Foundation of China under Grant No. 12101630. G. Cheng's research was partly funded by the NSF ¨C SCALE MoDL, Grant No. 2134209. The authors extend their gratitude to Prof. Hongjun Li for initial discussions, and to Dr. Ning Zhang and Dr. Yanpeng Li for their invaluable insights and suggestions.}