Subsampled Ensemble Can Improve Generalization Tail Exponentially

## Abstract

Abstract Ensemble learning is a popular technique to improve the accuracy of machine learning models. It hinges on the rationale that aggregating multiple weak models can lead to better models with lower variance and hence higher stability, especially for discontinuous base learners. In this paper, we provide a new perspective on ensembling. By selecting the best model trained on subsamples via majority voting, we can attain exponentially decaying tails for the excess risk, even if the base learner suffers from slow (i.e., polynomial) decay rates. This tail enhancement power of ensembling is agnostic to the underlying base learner and is stronger than variance reduction in the sense of exhibiting rate improvement. We demonstrate how our ensemble methods can substantially improve out-of-sample performances in a range of examples involving heavy-tailed data or intrinsically slow rates. Code for the proposed methods is available at https://github.com/mickeyhqian/VoteEnsemble .

## 1 Introduction

Ensemble learning (Dietterich, ; Zhou, ) is a class of methods to improve the accuracy of machine learning models. It comprises repeated training of models (the “base learners”), which are then aggregated through averaging or majority vote. In the literature, the main justification for ensemble methods, such as bootstrap aggregating (bagging) (Breiman, ) and boosting (Freund et al., ), pertains to bias/variance reduction or higher stability. This justification has been shown to be particularly relevant for certain U-statistics (Buja & Stuetzle, ) and models with hard-thresholding rules such as decision trees (Breiman, ; Drucker & Cortes, ).

Contrary to the established understanding, in this paper we present a new view of ensembling in offering an arguably stronger power than variance reduction: By suitably selecting the best base learners trained on random subsamples, ensembling leads to exponentially decaying excess risk tails. In particular, for general stochastic optimization problems that suffer from a slow, namely polynomial, decay in excess risk tails, ensembling can reduce these tails to an exponential decay rate. Thus, instead of the typical constant factor of improvement exhibited by variance reduction, our ensemble method offers a rate improvement, and moreover, the improvement is substantial.

In the following, we will first qualify our claims above by discussing how slow convergence can arise generically in machine learning and more general data-driven decision-making problems under heavy-tailed data. We then give intuition on our new ensembling perspective, proposed procedures, and the technicality involved in a full analysis.

#### Main results at a high level.

We begin by introducing a generic stochastic optimization problem

$$ $\min_{\theta\in\Theta}L(\theta):=\mathbb{E}\left[l(\theta,z)\right],$ (1) $$

where $\theta$ is the decision variable on space $\Theta$, $z\in\mathcal{Z}$ denotes the randomness governed by a probability distribution, and $l$ is the cost function. $n$ i.i.d. samples $\{z_{1},\ldots,z_{n}\}$ are available from the underlying distribution of $z$.
In machine learning, $\theta$ corresponds to model parameters, $\{z_{1},\ldots,z_{n}\}$ the training data, $l$ the loss function, and $L$ the population-level expected loss.
More generally, (1) encapsulates data-driven decision-making problems, namely the integration of data on $z$ into a downstream optimization task with overall cost function $l$ and prescriptive decision $\theta$.
These problems are increasingly prevalent in various industrial applications (Kamble et al., ; Bertsimas et al., ; Ghosal et al., ), such as in supply chain network design where $\theta$ may represent the decision to open processing facilities, $z$ the uncertain supply and demand, and $l$ the total cost of processing and transportation.

Given the data, we can train the model or decision with a learning algorithm that maps the data to an element in $\Theta$. This encompasses a wide range of methods, including machine learning training algorithms and data-driven approaches like sample average approximation (SAA) (Shapiro et al., ) and distributionally robust optimization (DRO) (Mohajerin Esfahani & Kuhn, ) in stochastic optimization. Our proposal and theory described below are agnostic to the choice of learning algorithm.

We characterize the generalization performance of a solution to (1), denoted by $\hat{\theta}$, via the tail probability bound on the excess risk or regret $L(\hat{\theta})-\min_{\theta\in\Theta}L(\theta)$, i.e., $\mathbb{P}(L(\hat{\theta})>\min_{\theta\in\Theta}L(\theta)+\delta)$ for some fixed $\delta>0$, where the probability is over both the data and training randomness. By a polynomially decaying generalization tail, we mean that

$$ $\mathbb{P}\Big{(}L(\hat{\theta})>\min_{\theta\in\Theta}L(\theta)+\delta\Big{)} \leq C_{1}n^{-\alpha}$ (2) $$

for some $\alpha>0$ and $C_{1}$ depends on $\delta$.
Such bounds are common under heavy-tailed data distributions (Kaňková & Houda, ; Jiang et al., ; Jiang & Li, ) due to slow concentration, which frequently arises in machine learning applications such as large language models (e.g, Jalalzai et al. (); Zhang et al. (); Cutkosky & Mehta () among others), finance (Mainik et al., ; Gilli & Këllezi, ) and physics (Fortin & Clusel, ; Michel & Chave, ), and are proved to be tight (Catoni, ) for empirical risk minimization (ERM) (Vapnik, ). As our key insight, our proposed ensembling methodology can improve the above to an exponential decay, i.e.,

$$ $\mathbb{P}\Big{(}L(\hat{\theta})>\min_{\theta\in\Theta}L(\theta)+\delta\Big{)} \leq C_{2}\gamma^{n/k},$ (3) $$

where $k$ is the subsampled data size and can be chosen at a slower rate in $n$, and $\gamma<1$ depends on $k,\delta$ such that $\gamma\to 0$ as $k\to\infty$.
Hence, when $k$ is properly chosen, the decay becomes exponential. This exponential acceleration is qualitatively different from the well-known variance reduction benefit of ensembling in several aspects.
First, variance reduction refers to the smaller variability in predictions from models trained on independent data sets, which has a more direct impact on the expected regret than the tail decay rate.
Second, the improvement by variance reduction is typical of a constant factor (e.g., Bühlmann & Yu () reported a reduction factor of $3$), thus affecting at best the constant $C_{1}$ in (2), whereas we obtain an order-of-magnitude improvement.

#### Main intuition.

To facilitate our explanation, let us first focus on discrete space $\Theta$. Our ensembling methodology uses a majority-vote mechanism at the model level: After repeatedly running the learning algorithm on subsamples to generate many models, we output the model that occurs most frequently. This implicitly solves a surrogate optimization problem over the same decision space $\Theta$ as (1) that maximizes the probability of being output by the learning algorithm. This conversion of the original general objective in (1) to a probability objective is the key: As an expectation of a random indicator function, the latter is uniformly bounded even if the original objective is heavy-tailed. Together with a bootstrap argument that establishes the closeness between subsample and full data, this in turn results in exponentially decaying tails for the regret.

For more general problems with continuous space, we replace the simple majority vote with a vote based on the likelihood of being $\epsilon$-optimal among all the generated models when evaluated on a random subsample. This avoids the degeneracy issue of using a simple majority vote for continuous problems while retaining similar (in fact, even stronger as we will see) guarantees. Regardless of discrete or continuous model space, our main insight on turning (2) into (3) applies. Moreover, in the discrete case, it turns out that not only the tail bound but also the average-case regret improves exponentially. This also explains why our improvement is particularly significant for discrete-decision problems in the experiments.

The rest of the paper is organized as follows. Section 2 presents our ensemble methods and their finite-sample bounds. Section 3 presents experimental results, and Section 4 discusses related work. Section 5 discusses limitations and concludes the paper.
A review of additional related work, technical proofs, and additional experimental results can be found in the appendix.

## 2 Methodology and theoretical guarantees

To solve (1) using data, we consider the generic learning algorithm in the form of a mapping

$$ $\mathcal{A}(z_{1},\ldots,z_{n};\omega):\mathcal{Z}^{n}\times\Omega\to\Theta$ $$

that takes in the training data $(z_{1},\ldots,z_{n})$ and outputs a model possibly under some algorithmic randomness $\omega$ that is independent of the data. Examples of $\omega$ include gradient sampling in stochastic first-order algorithms and feature/data subsampling in random forests. $\mathcal{A}(z_{1},\ldots,z_{n};\omega)$ serves as our base learner. For convenience, we omit $\omega$ to write $\mathcal{A}(z_{1},\ldots,z_{n})$ when no confusion arises.

### 2.1 A basic procedure

We first introduce a procedure called $\mathsf{MoVE}$ that applies to discrete solution or model space $\Theta$. $\mathsf{MoVE}$, which is formally described in Algorithm 1, repeatedly draws a total of $B$ subsamples from the data without replacement, learns a model via $\mathcal{A}$ on each subsample, and finally conducts a majority vote to output the most frequently subsampled model. Tie-breaking can be done arbitrarily.

Figure: Algorithm 1 Majority Vote Ensembling ($\mathsf{MoVE}$)

To understand $\mathsf{MoVE}$, we consider an optimization associated with the base learner $\mathcal{A}$

$$ $\max_{\theta\in\Theta}\ p_{k}(\theta):=\mathbb{P}\left(\theta=\mathcal{A}(z_{1 },\ldots,z_{k})\right),$ (4) $$

which maximizes the probability of a model being output by the base learner on $k$ i.i.d. observations. Here the probability $\mathbb{P}$ is with respect to both the training data and the algorithmic randomness. If $B=\infty$, $\mathsf{MoVE}$ essentially maximizes an empirical approximation of (4), i.e.

$$ $\max_{\theta\in\Theta}\ \mathbb{P}_{*}\left(\theta=\mathcal{A}(z^{*}_{1}, \ldots,z^{*}_{k})\right),$ (5) $$

where $\left(z^{*}_{1},\ldots,z^{*}_{k}\right)$ is a uniform random subsample from $\left(z_{1},\ldots,z_{n}\right)$, and $\mathbb{P}_{*}$ denotes the probability with respect to the algorithmic randomness and the subsampling randomness conditioned on $\left(z_{1},\ldots,z_{n}\right)$. With a finite $B<\infty$, extra Monte Carlo noises are introduced, leading to the following maximization problem

$$ $\max_{\theta\in\Theta}\ \frac{1}{B}\sum_{b=1}^{B}\mathbbm{1}(\theta=\mathcal{A }(z_{1}^{b},\ldots,z_{k}^{b})),$ (6) $$

which gives exactly the output of $\mathsf{MoVE}$. In other words, $\mathsf{MoVE}$ is a *bootstrap approximation* to the solution of (4). The following result materializes the intuition explained in the introduction on the conversion of the original potentially heavy-tailed problem (1) into a probability maximization (6) that possesses exponential bounds:

###### Theorem 1 (Finite-sample bound for Algorithm 1 )

Consider discrete decision space $\Theta$. Recall $p_{k}(\theta)$ defined in (4). Let $p_{k}^{\max}:=\max_{\theta\in\Theta}p_{k}(\theta)$, $\mathcal{E}_{k,\delta}:=\mathbb{P}(L(\mathcal{A}(z_{1},\ldots,z_{k}))>\min_{
\theta\in\Theta}L(\theta)+\delta)$ be the excess risk tail of $\mathcal{A}$, and

$$ $\eta_{k,\delta}:=p_{k}^{\max}-\mathcal{E}_{k,\delta}.$ (7) $$

For every $k\leq n$ and $\delta\geq 0$ such that $\eta_{k,\delta}>0$, the solution output by $\mathsf{MoVE}$ satisfies that

$$ $\displaystyle\ \mathbb{P}\left(L(\hat{\theta}_{n})>\min_{\theta\in\Theta}L( \theta)+\delta\right)$ (8) $\displaystyle\leq$ $\displaystyle\left|\Theta\right|\Bigg{[}\exp\left(-\dfrac{n}{2k}\cdot D_{ \operatorname{KL}}\left(p_{k}^{\max}-\dfrac{3\eta_{k,\delta}}{4}\Big{\|}p_{k}^ {\max}-\eta_{k,\delta}\right)\right)+2\exp\left(-\dfrac{n}{2k}\cdot D_{ \operatorname{KL}}\left(p_{k}^{\max}-\dfrac{\eta_{k,\delta}}{4}\Big{\|}p_{k}^{ \max}\right)\right)$ $\displaystyle+\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta_{k,\delta}^{2}}{\min \left(p_{k}^{\max},1-p_{k}^{\max}\right)+3\eta_{k,\delta}/4}\right)$ $\displaystyle+\mathbbm{1}\left(p_{k}^{\max}+\dfrac{\eta_{k,\delta}}{4}\leq 1 \right)\cdot\exp\Big{(}-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{ \max}+\dfrac{\eta_{k,\delta}}{4}\Big{\|}p_{k}^{\max}\right)-\dfrac{B}{24}\cdot \dfrac{\eta_{k,\delta}^{2}}{1-p_{k}^{\max}+\eta_{k,\delta}/4}\Big{)}\Bigg{]}.$ $$

In particular, if $\eta_{k,\delta}>4/5$, (8) is further bounded by

$$ $\left|\Theta\right|\left(3\min\left(e^{-2/5},C_{1}\max(1-p_{k}^{\max},\; \mathcal{E}_{k,\delta})\right)^{\frac{n}{C_{2}k}}+e^{-B/C_{3}}\right),$ (9) $$

where $C_{1},C_{2},C_{3}>0$ are universal constants, $\left|\Theta\right|$ denotes the cardinality of $\Theta$, and $D_{\operatorname{KL}}(p\|q):=p\ln\frac{p}{q}+(1-p)\ln\frac{1-p}{1-q}$ is the Kullback–Leibler divergence between two Bernoulli distributions with means $p$ and $q$.

Theorem 1 states that the excess risk tail of $\mathsf{MoVE}$ decays exponentially in the ratio $n/k$ and ensemble size $B$. The bound consists of three parts. The first part has two terms with the Kullback–Leibler (KL) divergences and arises from the bootstrap approximation of (4) with (5). The second part quantifies the Monte Carlo error in approximating (5) with a finite $B$. The third part comes from the interaction between the two sources of errors and is typically of higher order. The multiplier $\left|\Theta\right|$ in the bound is avoidable, e.g., via a space reduction as in our next algorithm.

The quantity $\eta_{k,\delta}$ plays two roles. First, it quantifies how suboptimality in the surrogate problem (4) propagates to the original problem (1) in that every $\eta_{k,\delta}$-optimal solution for (4) is $\delta$-optimal for (1). Second, $\eta_{k,\delta}$ is directly related to the excess risk tail $\mathcal{E}_{k,\delta}$ of the base learner, in addition to $p_{k}^{\max}$ that captures the concentration of the base learner on $\delta$-optimal solutions. Therefore, $\eta_{k,\delta}$ taking large values signals the situation where the base learner already generalizes well. In this case, (8) can be simplified to (9). The bound (9) suggests that our approach does not hurt the performance of an already high-performing base learner as its generalization power is inherited through the $\max(1-p_{k}^{\max},\;\mathcal{E}_{k,\delta})$ term in the bound. See Appendix B for a more detailed comparison.

The quantity $\eta_{k,\delta}$ also hints at how to choose the subsample size $k$. As long as $\eta_{k,\delta}$ is bounded away from 0 0, our bound decays exponentially fast. Therefore, $k$ can be chosen in such a way that the base learner outputs good models more often than bad ones in order for the exponential decay of our bound to take effect, but at the same time considerably smaller than $n$ to ensure the amount of acceleration. In the experiments, we choose $k=\max(10,n/200)$.

On the choice of $B$, note that the two KL divergences in the first part of the tail bound (8) are in general bounded below by $\mathcal{O}(\eta_{k,\delta}^{2})$ and so is the $\eta_{k,\delta}^{2}/(\min\left(p_{k}^{\max},1-p_{k}^{\max}\right)+3\eta_{k,
\delta}/4)$ in the second part as $\eta_{k,\delta}$ is no larger than $1$. Therefore using an ensemble size of $B=\mathcal{O}(n/k)$ is sufficient to control the Monte Carlo error to a similar magnitude as the data error.

###### Theorem 1 (Finite-sample bound for Algorithm 1 )

### 2.2 A more general procedure

We next present a more general procedure called $\mathsf{ROVE}$ that applies to continuous space where simple majority vote in Algorithm 1 can lead to degeneracy, i.e., all learned models appear exactly once in the pool. Moreover, this general procedure relaxes our dependence on $|\Theta|$ in the bound in Theorem 1.

Figure: Algorithm 2 Retrieval and $\epsilon$-Optimality Vote Ensembling ($\mathsf{ROVE}$ / $\mathsf{ROVEs}$)

$\mathsf{ROVE}$, displayed in Algorithm 2, proceeds initially the same as $\mathsf{MoVE}$ in repeatedly subsampling data and training the model using $\mathcal{A}$. However, in the aggregation step, instead of using a simple majority vote, $\mathsf{ROVE}$ outputs, among all the trained models, the one that has the highest likelihood of being $\epsilon$-optimal. This $\epsilon$-optimality avoids the degeneracy of the majority vote and, moreover, since we have restricted our output to the collection of trained models, the corresponding likelihood maximization is readily doable by simple enumeration. In addition, it helps reduce competition for votes among the best models as each subsample can now vote for multiple candidates, ensuring a high vote count for each of the top models even when there are many of them. This makes $\mathsf{ROVE}$ more effective than $\mathsf{MoVE}$ in the case of multiple (near) optima as our experiments will show. We have the following theoretical guarantees for Algorithm 2:

###### Theorem 2 (Finite-sample bound for Algorithm 2 )

Recall the tail $\mathcal{E}_{k,\delta}$ of the base excess risk from Theorem 1. Consider Algorithm 2 with data splitting, i.e., $\mathsf{ROVEs}$. Let $T_{k}(\cdot):=\mathbb{P}(\sup_{\theta\in\Theta}\lvert(1/k)\sum_{i=1}^{k}l(
\theta,z_{i})-L(\theta)\rvert>\cdot)$ be the tail function of the maximum deviation of the empirical objective estimate. For every $\delta>0$, if $\epsilon$ is chosen such that $\mathbb{P}\left(\epsilon\in[\underline{\epsilon},\overline{\epsilon}]\right)=1$ for some $0<\underline{\epsilon}\leq\overline{\epsilon}<\delta$ and $T_{k_{2}}\left((\delta-\overline{\epsilon})/2\right)+T_{k_{2}}\left(\underline
{\epsilon}/2\right)<1/5$, then

$$ $\displaystyle\mathbb{P}\Big{(}L(\hat{\theta}_{n})>\min_{\theta\in\Theta}L( \theta)+2\delta\Big{)}\leq$ $\displaystyle B_{1}\left[3\min\left(e^{-2/5},C_{1}T_{k_{2}}\left(\dfrac{\min( \underline{\epsilon},\delta-\overline{\epsilon})}{2}\right)\right)^{\frac{n}{2 C_{2}k_{2}}}+e^{-{B_{2}}/{C_{3}}}\right]$ (10) $\displaystyle+$ $\displaystyle\min\left(e^{-(1-\mathcal{E}_{k_{1},\delta})/C_{4}},C_{5}\mathcal {E}_{k_{1},\delta}\right)^{\frac{n}{2C_{6}k_{1}}}+e^{-B_{1}(1-\mathcal{E}_{k_{ 1},\delta})/C_{7}},$ $$

where $C_{1},C_{2},C_{3}$ are the same as those in Theorem 1, and $C_{4},C_{5},C_{6},C_{7}$ are universal constants.

Consider Algorithm 2 without data splitting, i.e., $\mathsf{ROVE}$, and discrete space $\Theta$. Assume $\lim_{k\to\infty}T_{k}(\delta)=0$ for all $\delta>0$. Then, for every fixed $\delta>0$, we have $\lim_{n\to\infty}\mathbb{P}(L(\hat{\theta}_{n})>\min_{\theta\in\Theta}L(\theta
)+2\delta)\to 0$, if $\limsup_{k\to\infty}\mathcal{E}_{k,\delta}<1$, $\mathbb{P}\left(\epsilon>\delta/2\right)\to 0$, $k_{1}\text{ and }k_{2}\to\infty$, $n/k_{1}\text{ and }n/k_{2}\to\infty$, and $B_{1},B_{2}\to\infty$ as $n\to\infty$.

Theorem 2 provides an exponential excess risk tail, regardless of discrete or continuous space. The first line in the bound (10) is inherited from the bound (9) for $\mathsf{MoVE}$ from majority to $\epsilon$-optimality vote. In particular, the multiplier $\left|\Theta\right|$ in (9) is now replaced by $B_{1}$, the number of retrieved models. The second line in (10) bounds the performance sacrifice due to the restriction to Phase I model candidates.

$\mathsf{ROVE}$ may be carried out with the data split between the two phases, in which case it’s referred to as $\mathsf{ROVEs}$. Data splitting makes the procedure theoretically more tractable by avoiding inter-dependency between the phases but sacrifices some statistical power from halving the data size. Empirically we find $\mathsf{ROVE}$ to be overall more effective.

The optimality threshold $\epsilon$ is allowed to be chosen in a data-driven way and the main goal guiding this choice is to be able to distinguish models of different qualities. In other words, $\epsilon$ should be chosen to create enough variability in the likelihood of being $\epsilon$-optimal across models. In our experiments, we find it a good strategy to choose an $\epsilon$ that leads to a maximum likelihood around $1/2$.

Lastly, our main theoretical results, Theorems 1 and 2, are derived using several novel techniques. First, we develop a sharper concentration result for U-statistics with binary kernels, improving upon standard Bernstein-type inequalities (e.g., Arcones (); Peel et al. ()). This refinement ensures the correct order of the bound, particularly (9), which captures the convergence of both the bootstrap approximation and the base learner, offering insights into the robustness of our methods for fast-converging base learners. Second, we perform a sensitivity analysis on the regret for the original problem (1) relative to the surrogate optimization (4), translating the superior generalization in the surrogate problem into accelerated convergence for the original. Finally, to establish asymptotic consistency for Algorithm 2 without data splitting, we develop a uniform law of large numbers (LLN) for the class of events of being $\epsilon$-optimal, using direct analysis of the second moment of the maximum deviation. Uniform LLNs are particularly challenging here because, unlike fixed classes in standard settings, this class dynamically depends on subsample size $k_{2}$ as $n\to\infty$.

###### Theorem 2 (Finite-sample bound for Algorithm 2 )

## 3 Numerical experiments

In this section, we numerically test Algorithm 1 ($\mathsf{MoVE}$), Algorithm 2 with ($\mathsf{ROVEs}$) and without ($\mathsf{ROVE}$) data splitting in training neural networks for regression problems and solving stochastic programs. Additional experimental results are provided in Appendix D due to space constraints. The code is available at: [https://anonymous.4open.science/r/vote_ensemble](https://anonymous.4open.science/r/vote_ensemble).

To empirically determine well-performing configurations for general use, we performed a comprehensive hyperparameter profiling of our algorithms in Appendix D.3. Below, we summarize the recommended configurations used in all experiments presented in this section (except Figure 4): 1) For discrete space $\Theta$, use $k=\max(10,n/200),B=200$ for $\mathsf{MoVE}$, and $k_{1}=k_{2}=\max(10,n/200),B_{1}=20,B_{2}=200$ for $\mathsf{ROVE}$ and $\mathsf{ROVEs}$; 2) For continuous space $\Theta$, use $k_{1}=\max(30,n/2),k_{2}=\max(30,n/200),B_{1}=50,B_{2}=200$ for $\mathsf{ROVE}$ and $\mathsf{ROVEs}$; 3) The $\epsilon$ in $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ is selected such that $\max_{\theta\in\mathcal{S}}(1/B_{2})\sum_{b=1}^{B_{2}}\mathbbm{1}(\theta\in
\widehat{\Theta}^{\epsilon,b}_{k_{2}})\approx 1/2$.

### 3.1 Neural networks for regression

We consider regression problems with multilayer perceptrons (MLPs) on both synthetic and real data. The base learning algorithm splits the data into training (70%) and validation (30%), and uses Adam to minimize mean squared error (MSE), with early stopping triggered when the validation improvement falls below 3% between epochs. The architecture details of the MLPs are provided in Appendix D.1. Note that $\mathsf{MoVE}$ is not included in this comparison as it’s applicable to discrete problems only.

#### Setup for Synthetic Data

Input-output pairs $(X,Y)$ are generated as $Y=(1/50)\cdot\sum_{j=1}^{50}\log(X_{j}+1)+\varepsilon$, where each $X_{j}$ is drawn independently from $\mathrm{Unif}(0,2+198(j-1)/49)$, and the noise $\varepsilon$ is independent of $X$ with zero mean. We consider both standard Gaussian noise and Pareto noise $\varepsilon=\varepsilon_{1}-\varepsilon_{2}$, where each $\varepsilon_{i}\sim\mathrm{Pareto}(2.1)$. The out-of-sample performance is estimated on a common test set of one million samples. Each algorithm is repeatedly applied to $200$ independently generated datasets to assess the average and tail performance.

#### Setup for Real Data

We use six datasets from the UCI Machine Learning Repository (Blake, ): Wine Quality (Cortez et al., ), Bike Sharing (Fanaee-T, ), Online News (Fernandes et al., ), Appliances Energy (Candanedo, ), Superconductivity (Hamidieh, ), and Gas Turbine Emission (gas, ). Each dataset is standardized (zero mean, unit variance). To evaluate the average and tail performance, we permute each dataset $100$ times, and each time use the first half for training and the second for testing.

Figure: (a) Pareto noise, $H=4$.
Refer to caption: extracted/5899831/img/NN_d50_l4_ES_plots/NN_d50_l4_ES_avg.png

Figure: (a) Appliances Energy.
Refer to caption: extracted/5899831/img/appliances_energy_plain_l4_plots/appliances_energy_plain_l4_cdf.png

#### Result.

As shown in Figure 1, in heavy-tailed noise settings (Figures 1a–1c), both $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ significantly outperform the base algorithm in terms of both expected out-of-sample MSE and tail performance under all sample sizes $n$. Notably, the performance improvement becomes more pronounced with deeper networks ($H=8$), indicating that the benefits of $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ are more apparent in models with higher expressiveness and lower bias.

In light-tailed settings (Figures 1d–1f), $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ show comparable expected out-of-sample performance to the base when $H=4$, but outperform the base as $H$ increases. Additionally, $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ outperform the base in tail probabilities even when $H=4$. This indicates that $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ provide better generalization as the model complexity grows even for light-tailed problems. Similar results for MLPs with $2$ and $6$ hidden layers can be found in Appendix D.4, where results on least squares regression and Ridge regression are also provided.

On real datasets (Figure 2), $\mathsf{ROVE}$ exhibits much lighter tails compared to the base on three out of six datasets, and similar tail behavior on the other three. $\mathsf{ROVEs}$, however, underperforms the base in these real-world scenarios, potentially due to the data split that compromises its statistical power.

### 3.2 Stochastic programs

#### Setup.

We consider four discrete stochastic programs: resource allocation, supply chain network design, maximum weight matching, and stochastic linear programming, alongside one continuous mean-variance portfolio optimization. All problems are designed to possess heavy-tailed uncertainties. For the stochastic linear program, instances with varying tail heaviness are explored to study its impact on algorithm performance. The base learning algorithm for all the problems is the SAA. Detailed descriptions of the problems are deferred to Appendix D.2 and results using DRO as the base algorithm are provided in Appendix D.4.

Figure: (a) Resource allocation.
Refer to caption: extracted/5899831/img/OptCIPlots/SSKP_SAA_error.png

Figure: (a) Influence of tail heaviness.
Refer to caption: extracted/5899831/img/LP_tail_influence.png

#### Result.

Figure 3 shows that our ensembling methods generally outperform the base algorithm in all cases, except for the linear program case (Figure 3(e)). Notably, $\mathsf{ROVE}$ still outperforms the base in the linear program case, demonstrating its robustness, while $\mathsf{MoVE}$ performs slightly worse than the base under small sample sizes. Comparing $\mathsf{ROVE}$ and $\mathsf{ROVEs}$, $\mathsf{ROVE}$ consistently exhibits superior performance than $\mathsf{ROVEs}$ in all cases.

When there is a unique optimal solution, $\mathsf{MoVE}$ and $\mathsf{ROVE}$ perform similarly, both generally better than $\mathsf{ROVEs}$, as seen in Figures 3(a)-3(d). However, in cases with multiple optima (Figures 3(e) and 4(a)), the performance of $\mathsf{MoVE}$ deteriorates while $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ stay strong. This is in accordance with our discussion on the advantage of $\epsilon$-optimality vote in Section 2.2. Additional results in Appendix D.4 shall further explain that optima multiplicity weakens the base learner for $\mathsf{MoVE}$ in the sense of decreasing the $\eta_{k,\delta}$ and hence inflating the tail bound in Theorem 1.

As shown in Figure 4(a), the performance gap between $\mathsf{ROVE}$, $\mathsf{ROVEs}$, and the base algorithm becomes increasingly significant as the tail of the uncertainty becomes heavier. This supports the effectiveness of $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ in handling heavy-tailed uncertainty, where the base algorithm’s performance suffers. Note that here $\mathsf{MoVE}$ behaves similarly as the base due to optima multiplicity.

The running time comparison in Figure 4(b) shows that, despite requiring multiple runs on subsamples, our ensembling methods do not introduce a significantly higher computational burden compared to running the base algorithm on the full sample, and can even be advantageous under large sample sizes. This is because, in problems like DRO (Ben-Tal et al., ; Mohajerin Esfahani & Kuhn, ) and two-stage stochastic programming, solving the optimization on the full sample often leads to a substantial increase in problem size, as the decision space and constraints grow at least linearly with the sample size. Subsampled optimizations, as performed in our approach, result in smaller, more manageable problems that can be solved more efficiently. Moreover, our theory indicates that solving more than
$\mathcal{O}(n/k)$ subsamples does not further improve generalization performance, ensuring that computational efficiency is maintained. Additionally, parallel processing of subsamples further reduces computational time.

Finally, among the three proposed ensemble methods, $\mathsf{ROVE}$ is the preferred choice over $\mathsf{MoVE}$ and $\mathsf{ROVEs}$ for general use as it’s applicable to both discrete and continuous problems and consistently delivers superior and stable performance across all scenarios.

## 4 Related work

This work is closely connected to various topics in optimization and machine learning, and we only review the most relevant ones. See Appendix A for additional literature review.

#### Ensemble learning.

Ensemble learning (Dietterich, ; Zhou, ; Sagi & Rokach, ) has been widely studied for improving model performance by combining multiple weak learners into strong ones. Popular ensemble methods include bagging (Breiman, ), boosting (Freund et al., ) and stacking (Wolpert, ; Džeroski & Ženko, ). Bagging enhances model stability by training models on different bootstrap samples and combining their predictions through majority voting or averaging, effectively reducing variance, especially for unstable learners like decision trees that underpin random forests (Breiman, ). Subagging (Bühlmann & Yu, ) is a variant of bagging that constructs the ensemble from subsamples in place of bootstrap samples. Boosting is a sequential process where each subsequent model corrects its predecessors’ errors, reducing both bias and variance (Ibragimov & Gusev, ; Ghosal & Hooker, ). Prominent boosting methods include AdaBoost (Freund et al., ), Stochastic Gradient Boosting (SGB) (Friedman, ; ), and Extreme Gradient Boosting (XGB) (Friedman et al., ) which differ in their approaches to weighting training data and hypotheses. Boosting is commonly used with decision trees as Gradient Boosted Decision Trees (GBDT), including XGBoost (Chen & Guestrin, ), LightGBM (Ke et al., ), and CatBoost (Hancock & Khoshgoftaar, ). Instead of using simple aggregation like weighted averaging or majority voting, stacking trains a model to combine base predictions in a more sophisticated way, further improving performance. A key procedural difference of our approach from these ensemble methods is that we perform majority voting at the model level, rather than at the prediction level, to select a single best model from the ensemble. As a result, our method consistently outputs models within the same space as the base learner, making it applicable to general stochastic optimization problems. In contrast, most existing ensemble methods yield aggregated models outside the base space. Additionally, compared to the bias/variance reduction of typical ensembles, our approach guarantees exponentially decaying excess risk tails and hence is particularly effective in settings with heavy-tailed noise.

#### Optimization and learning with heavy tails.

Optimization with heavy-tailed noises has garnered significant attention due to its relevance in traditional fields such as portfolio management (Mainik et al., ) and scheduling (Im et al., ), as well as emerging domains like large language models (Brown et al., ; Achiam et al., ). Tail bounds of most existing algorithms are guaranteed to decay exponentially under sub-Gaussian or uniformly bounded costs but deteriorate to a slow polynomial decay under heavy-tailedness (Kaňková & Houda, ; Jiang et al., ; Jiang & Li, ; Oliveira & Thompson, ).
For SAA or ERM, faster rates are possible under the small-ball (Mendelson, ; ; Roy et al., ) or Bernstein’s condition (Dinh et al., ) on the function class, while our approach is free from such conditions.
Considerable effort has been made to mitigate the adverse effects of heavy-tailedness with robust procedures among which the geometric median (Minsker, ), or more generally, median-of-means (MOM) (Lugosi & Mendelson, ; ) approach is most similar to ours.
The basic idea there is to estimate a true mean by dividing the data into disjoint subsamples, computing an estimate on each, and then taking the median. Lecué & Lerasle (); Lugosi & Mendelson (); Lecué & Lerasle () use MOM in estimating the expected cost and establish exponential tail bounds for the mean squared loss and convex function classes. Hsu & Sabato (; ) apply MOM directly on the solution level for continuous problems and require strong convexity from the cost to establish generalization bounds. Besides MOM, another approach estimates the expected cost via truncation (Catoni, ) and allows heavy tails for linear regression (Audibert & Catoni, ; Zhang & Zhou, ) or problems with uniformly bounded function classes (Brownlees et al., ), but is computationally intractable due to the truncation and thus more of theoretical interest. In contrast, our ensemble approach is a meta algorithm that acts on any learning algorithm to provide exponential tail bounds regardless of the underlying problem characteristics. Relatedly, various techniques such as gradient clipping (Cutkosky & Mehta, ; Gorbunov et al., ) and MOM (Puchkin et al., ) have been adopted in stochastic gradient descent (SGD) algorithms for handling heavy-tailed gradient noises, but their focus is the faster convergence of SGD rather than generalization.

#### Machine learning for optimization.

Learning to optimize (L2O) studies the use of machine learning in accelerating existing or discovering novel optimization algorithms. Much effort has been in training models via supervised or reinforcement learning to make critical algorithmic decisions such as cut selection (e.g., Deza & Khalil (); Tang et al. ()), search strategies (e.g., Khalil et al. (); He et al. (); Scavuzzo et al. ()), scaling (Berthold & Hendel, ), and primal heuristics (Shen et al., ) in mixed-integer optimization, or even directly generate high-quality solutions (e.g., neural combinatorial optimization pioneered by Bello et al. ()). See Chen et al. (; ); Bengio et al. (); Zhang et al. () for comprehensive surveys on L2O. This line of research is orthogonal to our goal, and L2O techniques can work as part of or directly serve as the base learning algorithm within our framework.

## 5 Conclusion and limitation

This paper introduces a novel ensemble technique that significantly improves generalization by aggregating base learners via majority voting. In particular, our approach converts polynomially decaying generalization tails into exponential decay, thus providing order-of-magnitude improvements as opposed to constant factor improvements exhibited by variance reduction. Extensive numerical experiments in both machine learning and stochastic programming validate its effectiveness, especially for scenarios with heavy-tailed data and slow convergence rates. This work underscores the powerful potential of our new ensemble approach across a broad range of machine learning applications.

While our method accelerates tail convergence, it may increase model bias, similar to other subsampling-based techniques like subagging (Bühlmann & Yu, ). This makes it best suited for applications with relatively low bias, e.g., when the model is sufficiently expressive.

## Appendix A Additional related work

#### Bagging for stochastic optimization.

Bagging has been adopted in stochastic optimization for various purposes. The most relevant line of works (Biggs et al., ; Perakis & Thayaparan, ; Wang et al., ; Biggs & Perakis, ) study mixed integer reformulations for stochastic optimization with bagging approximated objectives such as random forests and ensembles of neural networks with the ReLU activation. These works focus on computational tractability instead of generalization performance. Anderson & Nguyen () empirically evaluates several statistical techniques including bagging against the plain SAA and finds bagging advantageous for portfolio optimization problems. Birge () investigates a batch mean approach for continuous optimization that creates subsamples by dividing the data set into non-overlapping batches instead of resampling and aggregates SAA solutions on the subsamples via averaging, which is empirically demonstrated to reduce solution errors for constrained and high-dimensional problems. Another related batch of works (Lam & Qian, ; ; Chen & Woodruff, ; ; Eichhorn & Römisch, ) concern the use of bagging for constructing confidence bounds for generalization errors of data-driven solutions, but they do not attempt to improve generalization. Related to bagging, bootstrap has been utilized to quantify algorithmic uncertainties for randomized algorithms such as randomized least-squares algorithms (Lopes et al., ), randomized Newton methods (Chen & Lopes, ), and stochastic gradient descent (Fang et al., ; Zhong et al., ), which is orthogonal to our focus on generalization performance.

## Appendix B Implications of Theorem 1 for strong base learners

We provide a brief discussion of Theorem 1 applied to fast convergent base learners. Based on Theorem 1, the way $p_{k}^{\max}$ and $\mathcal{E}_{k,\delta}$ enter into (9) reflects how the generalization performance of the base learning algorithm is inherited by our framework. To explain, large $p_{k}^{\max}$ and small $\mathcal{E}_{k,\delta}$ correspond to better generalization of the base learning algorithm. This can be exploited by the bound (9) with the presence of $\max(1-p_{k}^{\max},\;\mathcal{E}_{k,\delta})$, which is captured with our sharper concentration of U-statistics with binary kernels. In particular, for base learning algorithms with fast generalization convergence, say $1-p_{k}^{\max}=\mathcal{O}(e^{-k})$ and $\mathcal{E}_{k,\delta}=\mathcal{O}(e^{-k})$ for simplicity, we have $C_{1}\max(1-p_{k}^{\max},\;\mathcal{E}_{k,\delta})=\mathcal{O}(e^{-k})$ and hence the first term in (9) becomes $\mathcal{O}(e^{-n})$ which matches the error of the base learning algorithm applied directly to the full data set.

## Appendix C Technical proofs

### C.1 Preliminaries

An important tool in the development of our theories is the U-statistic that naturally arises in subsampling without replacement. We first present the definition of U-statistic and its concentration properties.

###### Definition 1

Given the i.i.d. data set $\{z_{1},\ldots,z_{n}\}\subset\mathcal{Z}$ and a (not necessarily symmetric) kernel of order $k\leq n$ is a function $\kappa:\mathcal{Z}^{k}\to\mathbb{R}$ such that $\mathbb{E}\left[\left|\kappa(z_{1},\ldots,z_{k})\right|\right]<\infty$, the U-statistic associated with the kernel $\kappa$ is

$$ $U(z_{1},\ldots,z_{n})=\frac{1}{n(n-1)\cdots(n-k+1)}\sum_{1\leq i_{1},i_{2}, \cdots,i_{k}\leq n\ \mathrm{s.t.}\ i_{s}\neq i_{t}\ \forall 1\leq s<t\leq k} \kappa(z_{i_{1}},\ldots,z_{i_{k}}).$ $$

###### Lemma 1 (MGF dominance of U-statistics from Hoeffding ( 1963 ) )

For any integer $0<k\leq n$ and any kernel $\kappa(z_{1},\ldots,z_{k})$, let $U(z_{1},\ldots,z_{n})$ be the corresponding U-statistic defined in Definition 1, and

$$ $\bar{\kappa}(z_{1},\ldots,z_{n})=\frac{1}{\lfloor n/k\rfloor}\sum_{i=1}^{ \lfloor n/k\rfloor}\kappa(z_{k(i-1)+1},\ldots,z_{ki})$ (11) $$

be the average of the kernel across the first $\lfloor n/k\rfloor k$ data.
Then, for every $t\in\mathbb{R}$, it holds that

$$ $\mathbb{E}\left[\exp(tU)\right]\leq\mathbb{E}\left[\exp(t\bar{\kappa})\right].$ $$

Proof of Lemma 1. By symmetry, we have that

$$ $U(z_{1},\ldots,z_{n})=\frac{1}{n!}\sum_{\mathrm{bijection}\ \pi:[n]\rightarrow [n]}\bar{\kappa}(z_{\pi(1)},\ldots,z_{\pi(n)}),$ $$

where we denote $[n]:=\{1,\ldots,n\}$.
Then, by the convexity of the exponential function and Jensen’s inequality, we have that

$$ $\displaystyle\mathbb{E}\left[\exp(tU)\right]$ $\displaystyle=$ $\displaystyle\mathbb{E}\left[\exp\left(t\cdot\frac{1}{n!}\sum_{\mathrm{ bijection}\ \pi:[n]\rightarrow[n]}\bar{\kappa}(z_{\pi(1)},\ldots,z_{\pi(n)}) \right)\right]$ $\displaystyle\leq$ $\displaystyle\mathbb{E}\left[\frac{1}{n!}\sum_{\mathrm{bijection}\ \pi:[n] \rightarrow[n]}\exp\left(t\cdot\bar{\kappa}(z_{\pi(1)},\ldots,z_{\pi(n)}) \right)\right]$ $\displaystyle=$ $\displaystyle\mathbb{E}\left[\exp\left(t\cdot\bar{\kappa}(z_{1},\ldots,z_{n}) \right)\right].$ $$

This completes the proof.
$\Box$

Next, we present our sharper concentration bound for U-statistics with binary kernels:

###### Lemma 2 (Concentration bound for U-statistics with binary kernels)

Let $\kappa(z_{1},\dots,z_{k};\omega)$ be a $\{0,1\}$-valued kernel of order $k\leq n$ that possibly depends on additional randomness $\omega$ that is independent of the data $\left\{z_{1},\ldots,z_{n}\right\}$, $\kappa^{*}(z_{1},\dots,z_{k}):=\mathbb{E}\left[\kappa(z_{1},\dots,z_{k};\omega
)|z_{1},\dots,z_{k}\right]$, and $U(z_{1},\dots,z_{n})$ be the U-statistic associated with $\kappa^{*}$. Then, it holds that

$$ $\displaystyle\mathbb{P}\left(U-\mathbb{E}\left[\kappa\right]\geq\epsilon\right )\leq\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(\mathbb{E}\left[ \kappa\right]+\epsilon\|\mathbb{E}\left[\kappa\right]\right)\right),$ $\displaystyle\mathbb{P}\left(U-\mathbb{E}\left[\kappa\right]\leq-\epsilon \right)\leq\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(\mathbb{E} \left[\kappa\right]-\epsilon\|\mathbb{E}\left[\kappa\right]\right)\right),$ $$

where $D_{\operatorname{KL}}(p\|q):=p\ln\frac{p}{q}+(1-p)\ln\frac{1-p}{1-q}$ is the KL-divergence between two Bernoulli random variables with parameters $p$ and $q$, respectively.

Proof of Lemma 2. We first consider the direction $U-\mathbb{E}\left[\kappa\right]\geq\epsilon$. Let

$$ $\tilde{\kappa}^{*}:=\dfrac{1}{\hat{n}}\sum_{i=1}^{\hat{n}}\kappa^{*}(z_{k(i-1) +1},\ldots,z_{ki}),$ $$

and

$$ $\tilde{\kappa}:=\dfrac{1}{\hat{n}}\sum_{i=1}^{\hat{n}}\kappa(z_{k(i-1)+1}, \ldots,z_{ki};\omega_{i}),$ $$

where we use the shorthand notation $\hat{n}:=\lfloor\frac{n}{k}\rfloor$, and $\omega_{i}$’s are mutually independent and also independent from $\left\{z_{1},\ldots,z_{n}\right\}$.
Then, since $\mathbb{E}\left[\kappa\right]=\mathbb{E}\left[\kappa^{*}\right]$, for all $t>0$ it holds that

$$ $\displaystyle\mathbb{P}\left(U-\mathbb{E}\left[\kappa\right]\geq\epsilon\right)$ $\displaystyle=\mathbb{P}\left(\exp\left(tU\right)\geq\exp\left(t\left(\mathbb{ E}\left[\kappa\right]+\epsilon\right)\right)\right)$ (12) $\displaystyle\overset{(i)}{\leq}\exp\left(-t\left(\mathbb{E}\left[\kappa\right ]+\epsilon\right)\right)\cdot\mathbb{E}\left[\exp\left(tU\right)\right]$ $\displaystyle\overset{(ii)}{\leq}\exp\left(-t\left(\mathbb{E}\left[\kappa \right]+\epsilon\right)\right)\cdot\mathbb{E}\left[\exp\left(t\tilde{\kappa}^{ *}\right)\right]$ $\displaystyle\overset{(iii)}{\leq}\exp\left(-t\left(\mathbb{E}\left[\kappa \right]+\epsilon\right)\right)\cdot\mathbb{E}\left[\exp\left(t\tilde{\kappa} \right)\right],$ $$

where we apply the Markov inequality in $(i)$, step $(ii)$ is due to Lemma 1, and step $(iii)$ uses Jensen’s inequality and the convexity of the exponential function.
Due to independence, $\tilde{\kappa}$ can be viewed as the sample average of $\hat{n}$ i.i.d. Bernoulli random variables, i.e., $\tilde{\kappa}\sim\frac{1}{\hat{n}}\sum_{i=1}^{\hat{n}}\operatorname{Bernoulli
}\left(\mathbb{E}\left[\kappa\right]\right)$.
Hence, we have that

$$ $\displaystyle\mathbb{E}\left[\exp\left(t\tilde{\kappa}\right)\right]$ $\displaystyle=\mathbb{E}\left[\exp\left(\frac{t}{\hat{n}}\sum_{i=1}^{\hat{n}} \operatorname{Bernoulli}\left(\mathbb{E}\left[\kappa\right]\right)\right)\right]$ (13) $\displaystyle=\left(\mathbb{E}\left[\exp\left(\frac{t}{\hat{n}}\operatorname{ Bernoulli}\left(\mathbb{E}\left[\kappa\right]\right)\right)\right]\right)^{ \hat{n}}$ $\displaystyle=\left[(1-\mathbb{E}\left[\kappa\right])+\mathbb{E}\left[\kappa \right]\cdot\exp\left(\dfrac{t}{\hat{n}}\right)\right]^{\hat{n}},$ $$

where we use the moment-generating function of Bernoulli random variables in the last line.
Substituting (13) into (12), we have that

$$ $\mathbb{P}\left(U-\mathbb{E}\left[\kappa\right]\geq\epsilon\right)\leq\exp \left(-t\left(\mathbb{E}\left[\kappa\right]+\epsilon\right)\right)\cdot\left[( 1-\mathbb{E}\left[\kappa\right])+\mathbb{E}\left[\kappa\right]\cdot\exp\left( \dfrac{t}{\hat{n}}\right)\right]^{\hat{n}}=:f(t).$ (14) $$

Now, we consider minimizing $f(t)$ for $t>0$. Let $g(t)=\log f(t)$, then it holds that

$$ $g^{\prime}(t)=-(\mathbb{E}\left[\kappa\right]+\epsilon)+\dfrac{\mathbb{E}\left [\kappa\right]\cdot\exp\left(\frac{t}{\hat{n}}\right)}{(1-\mathbb{E}\left[ \kappa\right])+\mathbb{E}\left[\kappa\right]\cdot\exp\left(\frac{t}{\hat{n}} \right)}.$ $$

By setting $g^{\prime}(t)=0$, it is easy to verify that the minimum point of $f(t)$, denoted by $t^{\star}$, satisfies that

$$ $\displaystyle\mathbb{E}\left[\kappa\right]\cdot\exp\left(\frac{t}{\hat{n}} \right)\cdot(1-\mathbb{E}\left[\kappa\right]-\epsilon)=(1-\mathbb{E}\left[ \kappa\right])\cdot(\mathbb{E}\left[\kappa\right]+\epsilon)$ (15) $\displaystyle\Leftrightarrow$ $\displaystyle\exp(t)=\left[\dfrac{(1-\mathbb{E}\left[\kappa\right])\cdot( \mathbb{E}\left[\kappa\right]+\epsilon)}{\mathbb{E}\left[\kappa\right]\cdot(1- \mathbb{E}\left[\kappa\right]-\epsilon)}\right]^{\hat{n}}.$ $$

Substituting (15) into (14) gives

$$ $\displaystyle\mathbb{P}\left(U-\mathbb{E}\left[\kappa\right]\geq\epsilon\right)$ $\displaystyle\leq$ $\displaystyle\left(\dfrac{1-\mathbb{E}\left[\kappa\right]}{1-\mathbb{E}\left[ \kappa\right]-\epsilon}\right)^{\hat{n}}\cdot\left[\dfrac{\mathbb{E}\left[ \kappa\right]\cdot\left(1-\mathbb{E}\left[\kappa\right]-\epsilon\right)}{\left (1-\mathbb{E}\left[\kappa\right]\right)\left(\mathbb{E}\left[\kappa\right]+ \epsilon\right)}\right]^{\hat{n}\left(\mathbb{E}\left[\kappa\right]+\epsilon \right)}$ (16) $\displaystyle=$ $\displaystyle\left[\left(\dfrac{1-\mathbb{E}\left[\kappa\right]}{1-\mathbb{E} \left[\kappa\right]-\epsilon}\right)^{1-\mathbb{E}\left[\kappa\right]-\epsilon }\cdot\left(\dfrac{\mathbb{E}\left[\kappa\right]}{\mathbb{E}\left[\kappa\right ]+\epsilon}\right)^{\mathbb{E}\left[\kappa\right]+\epsilon}\right]^{\hat{n}}$ $\displaystyle=$ $\displaystyle\exp\left(-\hat{n}\cdot D_{\operatorname{KL}}\left(\mathbb{E} \left[\kappa\right]+\epsilon\|\mathbb{E}\left[\kappa\right]\right)\right).$ $$

Since $n/k\leq 2\hat{n}$, the first bound immediately follows from (16).

Since $D_{\operatorname{KL}}(p\|q)=D_{\operatorname{KL}}(1-p\|1-q)$, the bound for the reverse side $U-\mathbb{E}\left[\kappa\right]\leq-\epsilon$ then follows by applying the first bound to the flipped binary kernel $1-\kappa$ and $1-U$. This completes the proof of Lemma 2.
$\Box$

Next lemma gives lower bounds for KL divergences which help analyze the bounds in Lemma 2:

###### Lemma 3

Let $D_{\operatorname{KL}}(p\|q):=p\ln\frac{p}{q}+(1-p)\ln\frac{1-p}{1-q}$ be the KL-divergence between two Bernoulli random variables with parameters $p$ and $q$, respectively.
Then, it holds that

$$ $D_{\operatorname{KL}}(p\|q)\geq p\ln\frac{p}{q}+q-p.$ (17) $$

If $p\in[\gamma,1-\gamma]$ for some $\gamma\in(0,\frac{1}{2}]$, it also holds that

$$ $D_{\operatorname{KL}}(p\|q)\geq-\ln\left(2(q(1-q))^{\gamma}\right).$ (18) $$

Proof of Lemma 3. To show (17), some basic calculus shows that for any fixed $q$, the function $g(p):=(1-p)\ln\frac{1-p}{1-q}$ is convex in $p$, and we have that

$$ $g(q)=0,g^{\prime}(q)=-1.$ $$

Therefore $g(p)\geq g(q)+g^{\prime}(q)(p-q)=q-p$, which implies (17) immediately.

The lower bound (18) follows from

$$ $\displaystyle D_{\operatorname{KL}}(p\|q)$ $\displaystyle\geq$ $\displaystyle-p\ln q-(1-p)\ln(1-q)+\min_{p\in[\gamma,1-\gamma]}\{p\ln p+(1-p) \ln(1-p)\}$ $\displaystyle\geq$ $\displaystyle-\gamma\ln q-\gamma\ln(1-q)-\ln 2=-\ln(2(q(1-q))^{\gamma}).$ $$

This completes the proof of Lemma 3.
$\Box$

To incorporate all the proposed algorithms in a unified theoretical framework, we consider a set-valued mapping

$$ $\mathbb{A}(z_{1},\ldots,z_{k};\omega):\mathcal{Z}^{k}\times\Omega\to 2^{\Theta}$ (19) $$

where $\omega$ denotes algorithmic randomness that is independent of the data $\left\{z_{1},\ldots,z_{k}\right\}$. Each of our proposed algorithms attempts to solve the probability-maximization problem

$$ $\max_{\theta\in\Theta}\ \hat{p}_{k}(\theta):=\mathbb{P}_{*}\left(\theta\in \mathbb{A}(z^{*}_{1},\ldots,z^{*}_{k};\omega)\right),$ (20) $$

for a certain choice of $\mathbb{A}$, where $\left\{z^{*}_{1},\ldots,z^{*}_{k}\right\}$ is subsampled from the i.i.d. data $\left\{z_{1},\ldots,z_{n}\right\}$ uniformly without replacement, and $\mathbb{P}_{*}$ denotes the probability with respect to the algorithmic randomness $\omega$ and the subsampling randomness conditioned on the data. Note that this problem is an empirical approximation of the problem

$$ $\max_{\theta\in\Theta}\ p_{k}(\theta):=\mathbb{P}\left(\theta\in\mathbb{A}(z_{ 1},\ldots,z_{k};\omega)\right).$ (21) $$

The problem actually solved with a finite number of subsamples is

$$ $\max_{\theta\in\Theta}\ \bar{p}_{k}(\theta):=\frac{1}{B}\sum_{b=1}^{B}\mathbbm {1}(\theta\in\mathbb{A}(z_{1}^{b},\ldots,z_{k}^{b};\omega_{b})).$ (22) $$

Specifically, Algorithm 1 uses

$$ $\mathbb{A}(z^{*}_{1},\ldots,z^{*}_{k};\omega)=\left\{\mathcal{A}(z^{*}_{1}, \ldots,z^{*}_{k};\omega)\right\}$ (23) $$

where $\mathcal{A}$ denotes the base learning algorithm, and Algorithm 2 uses

$$ $\mathbb{A}(z^{*}_{1},\ldots,z^{*}_{k_{2}};\omega)=\left\{\theta\in\mathcal{S}: \frac{1}{k_{2}}\sum_{i=1}^{k_{2}}l(\theta,z^{*}_{i})\leq\min_{\theta^{\prime} \in\mathcal{S}}\frac{1}{k_{2}}\sum_{i=1}^{k_{2}}l(\theta^{\prime},z^{*}_{i})+ \epsilon\right\}$ (24) $$

conditioned on the solution set $\mathcal{S}$ retrieved in Phase I. Note that no algorithmic randomness is involved in (24) once the set $\mathcal{S}$ is given. We define:

###### Definition 2

For any $\delta\in[0,1]$, let

$$ $\mathcal{P}^{\delta}_{k}:=\{\theta\in\Theta:p_{k}(\theta)\geq\max_{\theta^{ \prime}\in\Theta}p_{k}(\theta^{\prime})-\delta\}$ (25) $$

be the set of $\delta$-optimal solutions of problem (21). Let

$$ $\theta_{k}^{\max}\in\operatorname*{arg\,max}_{\theta\in\Theta}p_{k}(\theta)$ $$

be a solution with maximum probability that is chosen in a unique manner if there are multiple such solutions. Let

$$ $\widehat{\mathcal{P}}^{\delta}_{k}:=\{\theta\in\Theta:\hat{p}_{k}(\theta)\geq \hat{p}_{k}(\theta_{k}^{\max})-\delta\}$ (26) $$

be the set of $\delta$-optimal solutions relative to $\theta_{k}^{\max}$ for problem (20).

and

###### Definition 3

Let

$$ $\Theta^{\delta}:=\left\{\theta\in\Theta:L(\theta)\leq\min_{\theta^{\prime}\in \Theta}L(\theta^{\prime})+\delta\right\}$ (27) $$

be the set of $\delta$-optimal solutions of problem (1). In particular, $\Theta^{0}$ represents the set of optimal solutions. Let

$$ $\widehat{\Theta}_{k}^{\delta}:=\left\{\theta\in\Theta:\frac{1}{k}\sum_{i=1}^{k }l(\theta,z_{i})\leq\min_{\theta^{\prime}\in\Theta}\frac{1}{k}\sum_{i=1}^{k}l( \theta^{\prime},z_{i})+\delta\right\}$ (28) $$

be the set of $\delta$-optimal solutions of the SAA with i.i.d. data $(z_{1},\ldots,z_{k})$.

###### Definition 1

###### Lemma 1 (MGF dominance of U-statistics from Hoeffding ( 1963 ) )

###### Lemma 2 (Concentration bound for U-statistics with binary kernels)

###### Lemma 3

###### Definition 2

###### Definition 3

### C.2 Proof for Theorem 1

We consider Algorithm 3, a more general version of Algorithm 1 that operates on the set-valued learning algorithm $\mathbb{A}$ in (19) and reduces to exactly Algorithm 1 in the special case (23). Again we omit the algorithmic randomness $\omega$ in $\mathbb{A}$ for convenience.

Figure: Algorithm 3 Majority Vote Ensembling for Set-Valued Learning Algorithms

We have the following finite-sample result for Algorithm 3:

###### Theorem 3 (Finite-sample bound for Algorithm 3 )

Consider discrete decision space $\Theta$. Recall $p_{k}(\theta)$ defined in (21). Let $p_{k}^{\max}:=\max_{\theta\in\Theta}p_{k}(\theta)$ and

$$ $\bar{\eta}_{k,\delta}:=p_{k}^{\max}-\max_{\theta\in\Theta\backslash\Theta^{ \delta}}p_{k}(\theta),$ (29) $$

where $\max_{\theta\in\Theta\backslash\Theta^{\delta}}p_{k}(\theta)$ evaluates to 0 0 if $\Theta\backslash\Theta^{\delta}$ is empty. For every $k\leq n$ and $\delta\geq 0$ such that $\bar{\eta}_{k,\delta}>0$, the solution output by Algorithm 3 satisfies that

$$ $\displaystyle\mathbb{P}\left(L(\hat{\theta}_{n})>\min_{\theta\in\Theta}L( \theta)+\delta\right)$ (30) $\displaystyle\leq$ $\displaystyle\left|\Theta\right|\Bigg{[}\exp\left(-\dfrac{n}{2k}\cdot D_{ \operatorname{KL}}\left(p_{k}^{\max}-\dfrac{3\eta}{4}\Big{\|}p_{k}^{\max}-\eta \right)\right)+2\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k} ^{\max}-\dfrac{\eta}{4}\Big{\|}p_{k}^{\max}\right)\right)+$ $\displaystyle\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{\min\left(p_{k}^{ \max},1-p_{k}^{\max}\right)+3\eta/4}\right)+$ $\displaystyle\mathbbm{1}\left(p_{k}^{\max}+\dfrac{\eta}{4}\leq 1\right)\cdot \exp\Big{(}-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{\max}+\dfrac{ \eta}{4}\Big{\|}p_{k}^{\max}\right)-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{1-p_{k} ^{\max}+\eta/4}\Big{)}\Bigg{]}$ $$

for every $\eta\in(0,\bar{\eta}_{k,\delta}]$. In particular, if $\bar{\eta}_{k,\delta}>4/5$, (30) is further bounded by

$$ $\left|\Theta\right|\left(3\min\left(e^{-2/5},C_{1}\max(1-p_{k}^{\max},\;\max_{ \theta\in\Theta\backslash\Theta^{\delta}}p_{k}(\theta))\right)^{\frac{n}{C_{2} k}}+\exp\left(-\frac{B}{C_{3}}\right)\right),$ (31) $$

where $C_{1},C_{2},C_{3}>0$ are universal constants, and $D_{\operatorname{KL}}(p\|q):=p\ln\frac{p}{q}+(1-p)\ln\frac{1-p}{1-q}$ is the Kullback–Leibler divergence between two Bernoulli distributions with means $p$ and $q$.

Proof of Theorem 3. We first prove excess risk tail bounds for the problem (21), split into two lemmas, Lemmas 4 and 5 below.

###### Lemma 4

Consider discrete decision space $\Theta$. Recall from Definition 2 that $p_{k}^{\max}=p_{k}(\theta_{k}^{\max})$ holds for $\theta_{k}^{\max}$. For every $0\leq\epsilon\leq\delta\leq p_{k}^{\max}$, it holds that

$$ $\displaystyle\mathbb{P}\left(\widehat{\mathcal{P}}^{\epsilon}_{k}\not\subseteq \mathcal{P}^{\delta}_{k}\right)\leq\left|\Theta\right|$ $\displaystyle\left[\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_ {k}^{\max}-\dfrac{\delta+\epsilon}{2}\Big{\|}p_{k}^{\max}-\delta\right)\right)\right.$ $\displaystyle+$ $\displaystyle\left.\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_ {k}^{\max}-\dfrac{\delta-\epsilon}{2}\Big{\|}p_{k}^{\max}\right)\right)\right].$ $$

Proof of Lemma 4. By Definition 2, we observe the following equivalence

$$ $\left\{\widehat{\mathcal{P}}^{\epsilon}_{k}\not\subseteq\mathcal{P}^{\delta}_{ k}\right\}=\bigcup_{\theta\in\Theta\backslash\mathcal{P}^{\delta}_{k}}\left\{ \theta\in\widehat{\mathcal{P}}^{\epsilon}_{k}\right\}=\bigcup_{\theta\in\Theta \backslash\mathcal{P}^{\delta}_{k}}\left\{{\hat{p}}_{k}(\theta)\geq{\hat{p}}_{ k}\left(\theta_{k}^{\max}\right)-\epsilon\right\}.$ $$

Hence, by the union bound, it holds that

$$ $\mathbb{P}\left(\widehat{\mathcal{P}}^{\epsilon}_{k}\not\subseteq\mathcal{P}^{ \delta}_{k}\right)\leq\sum_{\theta\in\Theta\backslash\mathcal{P}^{\delta}_{k}} \mathbb{P}\left({\hat{p}}_{k}(\theta)\geq{\hat{p}}_{k}\left(\theta_{k}^{\max} \right)-\epsilon\right).$ $$

We further bound the probability $\mathbb{P}\left(\left\{{\hat{p}}_{k}(\theta)\geq{\hat{p}}_{k}\left(\theta_{k}^
{\max}\right)-\epsilon\right\}\right)$ as follows

$$ $\displaystyle\mathbb{P}\left({\hat{p}}_{k}(\theta)\geq{\hat{p}}_{k}\left( \theta_{k}^{\max}\right)-\epsilon\right)$ (32) $\displaystyle\leq$ $\displaystyle\mathbb{P}\left(\left\{{\hat{p}}_{k}(\theta)\geq p_{k}(\theta_{k} ^{\max})-\dfrac{\delta+\epsilon}{2}\right\}\cap\left\{{\hat{p}}_{k}\left( \theta_{k}^{\max}\right)\leq p_{k}(\theta_{k}^{\max})-\dfrac{\delta-\epsilon}{ 2}\right\}\right)$ $\displaystyle\leq$ $\displaystyle\mathbb{P}\left({\hat{p}}_{k}(\theta)\geq p_{k}(\theta_{k}^{\max} )-\dfrac{\delta+\epsilon}{2}\right)+\mathbb{P}\left({\hat{p}}_{k}\left(\theta_ {k}^{\max}\right)\leq p_{k}(\theta_{k}^{\max})-\dfrac{\delta-\epsilon}{2} \right).$ $$

On one hand, the first probability in (32) is solely determined by and increasing in $p_{k}(\theta)=\mathbb{E}\left[{\hat{p}}_{k}(\theta)\right]$. On the other hand, we have $p_{k}(\theta)<p_{k}\left(\theta_{k}^{\max}\right)-\delta$ for every $\theta\in\Theta\backslash\mathcal{P}^{\delta}_{k}$ by the definition of $\mathcal{P}^{\delta}_{k}$. Therefore we can slightly abuse the notation to write

$$ $\displaystyle\mathbb{P}\left({\hat{p}}_{k}(\theta)\geq{\hat{p}}_{k}\left( \theta_{k}^{\max}\right)-\epsilon\right)$ $\displaystyle\leq$ $\displaystyle\mathbb{P}\left({\hat{p}}_{k}(\theta)\geq p_{k}(\theta_{k}^{\max} )-\dfrac{\delta+\epsilon}{2}\Big{|}p_{k}(\theta)=p_{k}(\theta_{k}^{\max})- \delta\right)$ $\displaystyle\hskip 8.61108pt+\mathbb{P}\left({\hat{p}}_{k}\left(\theta_{k}^{ \max}\right)\leq p_{k}(\theta_{k}^{\max})-\dfrac{\delta-\epsilon}{2}\right)$ $\displaystyle\leq$ $\displaystyle\mathbb{P}\left({\hat{p}}_{k}(\theta)-p_{k}(\theta)\geq\dfrac{ \delta-\epsilon}{2}\Big{|}p_{k}(\theta)=p_{k}(\theta_{k}^{\max})-\delta\right)$ $\displaystyle\hskip 8.61108pt+\mathbb{P}\left({\hat{p}}_{k}\left(\theta_{k}^{ \max}\right)-p_{k}(\theta_{k}^{\max})\leq-\dfrac{\delta-\epsilon}{2}\right).$ $$

Note that, with $\kappa(z_{1},\ldots,z_{k};\omega):=\mathbf{1}\left(\theta\in\mathbb{A}(z_{1},
\ldots,z_{k};\omega)\right)$, the probability ${\hat{p}}_{k}(\theta)$ can be viewed as a U-statistic with the kernel $\kappa^{*}(z_{1},\ldots,z_{k}):=\mathbb{E}\left[\kappa(z_{1},\ldots,z_{k};
\omega)|z_{1},\ldots,z_{k}\right]$.
A similar representation holds for ${\hat{p}}_{k}\left(\theta_{k}^{\max}\right)$ as well. Therefore, we can apply Lemma 2 to conclude that

$$ $\displaystyle\mathbb{P}\left(\widehat{\mathcal{P}}^{\epsilon}_{k}\not\subseteq \mathcal{P}^{\delta}_{k}\right)$ $\displaystyle\leq\sum_{\theta\in\Theta\backslash\mathcal{P}^{\delta}_{k}} \mathbb{P}\left({\hat{p}}_{k}(\theta)\geq{\hat{p}}_{k}\left(\theta_{k}^{\max} \right)-\epsilon\right)$ $\displaystyle\leq\left|\Theta\backslash\mathcal{P}^{\delta}_{k}\right|\left[ \mathbb{P}\left({\hat{p}}_{k}(\theta)-p_{k}(\theta)\geq\dfrac{\delta-\epsilon} {2}\Big{|}p_{k}(\theta)=p_{k}(\theta_{k}^{\max})-\delta\right)\right.$ $\displaystyle\hskip 51.6665pt+\left.\mathbb{P}\left(p_{k}\left(\theta_{k}^{ \max}\right)-{\hat{p}}_{k}\left(\theta_{k}^{\max}\right)\leq-\dfrac{\delta- \epsilon}{2}\right)\right]$ $\displaystyle\leq\left|\Theta\right|\left[\exp\left(-\dfrac{n}{2k}\cdot D_{ \operatorname{KL}}\left(p_{k}\left(\theta_{k}^{\max}\right)-\delta+\dfrac{ \delta-\epsilon}{2}\Big{\|}p_{k}\left(\theta_{k}^{\max}\right)-\delta\right) \right)\right.$ $\displaystyle\hskip 54.06006pt+\left.\exp\left(-\dfrac{n}{2k}\cdot D_{ \operatorname{KL}}\left(p_{k}\left(\theta_{k}^{\max}\right)-\dfrac{\delta- \epsilon}{2}\Big{\|}p_{k}\left(\theta_{k}^{\max}\right)\right)\right)\right],$ $$

which completes the proof of Lemma 4.
$\Box$

###### Lemma 5

Consider discrete decision space $\Theta$. For every $\epsilon\in[0,1]$ it holds for the solution output by Algorithm 3 that

$$ $\mathbb{P}_{*}\left(\hat{\theta}_{n}\notin\widehat{\mathcal{P}}^{\epsilon}_{k} \right)\leq\left|\Theta\right|\cdot\exp\left(-\dfrac{B}{6}\cdot\dfrac{\epsilon ^{2}}{\min\left(\hat{p}_{k}(\theta_{k}^{\max}),1-\hat{p}_{k}(\theta_{k}^{\max} )\right)+\epsilon}\right),$ $$

where $\left|\cdot\right|$ denotes the cardinality of a set and $\mathbb{P}_{*}$ denotes the probability with respect to both the resampling randomness conditioned on the observations and the algorithmic randomness.

Proof of Lemma 5. We observe that $\bar{p}_{k}(\theta)$ is an conditionally unbiased estimator for $\hat{p}_{k}(\theta)$, i.e., $\mathbb{E}_{*}\left[\bar{p}_{k}(\theta)\right]={\hat{p}}_{k}(\theta)$. We can express the difference between $\bar{p}_{k}(\theta)$ and $\bar{p}_{k}(\theta_{k}^{\max})$ as the sample average

$$ $\bar{p}_{k}(\theta)-\bar{p}_{k}(\theta_{k}^{\max})=\dfrac{1}{B}\sum_{b=1}^{B} \left[\mathbbm{1}(\theta\in\mathbb{A}(z_{1}^{b},\ldots,z_{k}^{b}))-\mathbbm{1} (\theta_{k}^{\max}\in\mathbb{A}(z_{1}^{b},\ldots,z_{k}^{b}))\right],$ $$

whose expectation is equal to ${\hat{p}}_{k}(\theta)-{\hat{p}}_{k}(\theta^{\max}_{k})$. We denote by

$$ $\mathbbm{1}_{\theta}^{*}:=\mathbbm{1}(\theta\in\mathbb{A}(z_{1}^{*},\ldots,z_{ k}^{*}))\text{\ for }\theta\in\Theta$ $$

for convenience, where $(z_{1}^{*},\ldots,z_{k}^{*})$ represents a random subsample. Then by Bernstein’s inequality, we have every $t\geq 0$ that

$$ $\displaystyle\mathbb{P}_{*}\left(\bar{p}_{k}(\theta)-\bar{p}_{k}(\hat{\theta}_ {k}^{\max})-(\hat{p}_{k}(\theta)-\hat{p}_{k}(\theta_{k}^{\max}))\geq t\right)$ (33) $\displaystyle\leq$ $\displaystyle\exp\left(-B\cdot\dfrac{t^{2}}{2\mathrm{Var}_{*}(\mathbbm{1}_{ \theta}^{*}-\mathbbm{1}_{\theta_{k}^{\max}}^{*})+4/3\cdot t}\right).$ $$

Since

$$ $\displaystyle\mathrm{Var}_{*}(\mathbbm{1}_{\theta}^{*}-\mathbbm{1}_{\theta_{k} ^{\max}}^{*})$ $\displaystyle\leq$ $\displaystyle\mathbb{E}_{*}\left[(\mathbbm{1}_{\theta}^{*}-\mathbbm{1}_{\theta _{k}^{\max}}^{*})^{2}\right]$ $\displaystyle\leq$ $\displaystyle\hat{p}_{k}(\theta)+\hat{p}_{k}(\theta_{k}^{\max})\leq 2\hat{p}_{ k}(\theta_{k}^{\max}),$ $$

and

$$ $\displaystyle\mathrm{Var}_{*}(\mathbbm{1}_{\theta}^{*}-\mathbbm{1}_{\theta_{k} ^{\max}}^{*})$ $\displaystyle\leq$ $\displaystyle\mathrm{Var}_{*}(1-\mathbbm{1}_{\theta}^{*}-1+\mathbbm{1}_{\theta _{k}^{\max}}^{*})$ $\displaystyle\leq$ $\displaystyle\mathbb{E}_{*}\left[(1-\mathbbm{1}_{\theta}^{*}-1+\mathbbm{1}_{ \theta_{k}^{\max}}^{*})^{2}\right]$ $\displaystyle\leq$ $\displaystyle 1-\hat{p}_{k}(\theta)+1-\hat{p}_{k}(\theta_{k}^{\max})\leq 2(1- \hat{p}_{k}(\theta)),$ $$

we have $\mathrm{Var}_{*}(\mathbbm{1}_{\theta}^{*}-\mathbbm{1}_{\theta_{k}^{\max}}^{*})
\leq 2\min(\hat{p}_{k}(\theta_{k}^{\max}),1-\hat{p}_{k}(\theta))$. Substituting this bound to (33) and taking $t=\hat{p}_{k}(\theta_{k}^{\max})-\hat{p}_{k}(\theta)$ lead to

$$ $\displaystyle\mathbb{P}_{*}\left(\bar{p}_{k}(\theta)-\bar{p}_{k}(\hat{\theta}_ {k}^{\max})\geq 0\right)$ $\displaystyle\leq$ $\displaystyle\exp\left(-B\cdot\dfrac{(\hat{p}_{k}(\theta_{k}^{\max})-\hat{p}_{ k}(\theta))^{2}}{4\min(\hat{p}_{k}(\theta_{k}^{\max}),1-\hat{p}_{k}(\theta))+4 /3\cdot(\hat{p}_{k}(\theta_{k}^{\max})-\hat{p}_{k}(\theta))}\right)$ $\displaystyle\leq$ $\displaystyle\exp\left(-B\cdot\dfrac{(\hat{p}_{k}(\theta_{k}^{\max})-\hat{p}_{ k}(\theta))^{2}}{4\min(\hat{p}_{k}(\theta_{k}^{\max}),1-\hat{p}_{k}(\theta_{k} ^{\max}))+16/3\cdot(\hat{p}_{k}(\theta_{k}^{\max})-\hat{p}_{k}(\theta))}\right)$ $\displaystyle\leq$ $\displaystyle\exp\left(-\dfrac{B}{6}\cdot\dfrac{(\hat{p}_{k}(\theta_{k}^{\max} )-\hat{p}_{k}(\theta))^{2}}{\min(\hat{p}_{k}(\theta_{k}^{\max}),1-\hat{p}_{k}( \theta_{k}^{\max}))+\hat{p}_{k}(\theta_{k}^{\max})-\hat{p}_{k}(\theta)}\right).$ $$

Therefore, we have that

$$ $\displaystyle\mathbb{P}_{*}\left(\hat{\theta}_{n}\notin\widehat{\mathcal{P}}^{ \epsilon}_{k}\right)$ $\displaystyle=\mathbb{P}_{*}\left(\bigcup_{\theta\in\Theta\backslash\widehat{ \mathcal{P}}^{\epsilon}_{k}}\left\{\bar{p}_{k}(\theta)=\max_{\theta^{\prime} \in\Theta}\bar{p}_{k}(\theta^{\prime})\right\}\right)$ $\displaystyle\leq\sum_{\theta\in\Theta\backslash\widehat{\mathcal{P}}^{ \epsilon}_{k}}\mathbb{P}_{*}\left(\bar{p}_{k}(\theta)=\max_{\theta^{\prime}\in \Theta}\bar{p}_{k}(\theta^{\prime})\right)$ $\displaystyle\leq\sum_{\theta\in\Theta\backslash\widehat{\mathcal{P}}^{ \epsilon}_{k}}\mathbb{P}_{*}\left(\bar{p}_{k}(\theta)\geq\bar{p}_{k}(\theta_{k }^{\max})\right)$ $\displaystyle\leq\sum_{\theta\in\Theta\backslash\widehat{\mathcal{P}}^{ \epsilon}_{k}}\exp\left(-\dfrac{B}{6}\cdot\dfrac{(\hat{p}_{k}(\theta_{k}^{\max })-\hat{p}_{k}(\theta))^{2}}{\min(\hat{p}_{k}(\theta_{k}^{\max}),1-\hat{p}_{k} (\theta_{k}^{\max}))+\hat{p}_{k}(\theta_{k}^{\max})-\hat{p}_{k}(\theta)}\right).$ $$

Note that the function $x^{2}/(\min(\hat{p}_{k}(\theta_{k}^{\max}),1-\hat{p}_{k}(\theta_{k}^{\max}))+x)$ in $x\in[0,1]$ is monotonically increasing and that $\hat{p}_{k}(\theta_{k}^{\max})-\hat{p}_{k}(\theta)>\epsilon$ for all $\theta\in\Theta\backslash\widehat{\mathcal{P}}^{\epsilon}_{k}$. Therefore, we can further bound the probability as

$$ $\displaystyle\mathbb{P}_{*}\left(\hat{\theta}_{n}\notin\widehat{\mathcal{P}}^{ \epsilon}_{k}\right)$ $\displaystyle\leq$ $\displaystyle\left|\Theta\backslash\widehat{\mathcal{P}}^{\epsilon}_{k}\right| \cdot\exp\left(-\dfrac{B}{6}\cdot\dfrac{\epsilon^{2}}{\min\left(\hat{p}_{k}^{ \max},1-\hat{p}_{k}^{\max}\right)+\epsilon}\right).$ $$

Noting that $\left|\Theta\backslash\widehat{\mathcal{P}}^{\epsilon}_{k}\right|\leq\left|
\Theta\right|$ completes the proof of Lemma 5.
$\Box$

We are now ready for the proof of Theorem 3. We first note that, if $\bar{\eta}_{k,\delta}>0$,
it follows from Definition 2 that

$$ $\mathcal{P}_{k}^{\eta}\subseteq\Theta^{\delta}\text{\ for any\ }\eta\in(0,\; \bar{\eta}_{k,\delta}).$ $$

Therefore, for any $\eta\in(0,\;\bar{\eta}_{k,\delta})$, we can write that

$$ $\displaystyle\mathbb{P}\left(\hat{\theta}_{n}\notin\Theta^{\delta}\right)\leq \mathbb{P}\left(\hat{\theta}_{n}\notin\mathcal{P}^{\eta}_{k}\right)$ $\displaystyle\leq$ $\displaystyle\mathbb{P}\left(\left\{\hat{\theta}_{n}\notin\widehat{\mathcal{P} }^{\eta/2}_{k}\right\}\cup\left\{\widehat{\mathcal{P}}^{\eta/2}_{k}\not \subseteq\mathcal{P}^{\eta}_{k}\right\}\right)$ (34) $\displaystyle\leq$ $\displaystyle\mathbb{P}\left(\hat{\theta}_{n}\notin\widehat{\mathcal{P}}^{\eta /2}_{k}\right)+\mathbb{P}\left(\widehat{\mathcal{P}}^{\eta/2}_{k}\not\subseteq \mathcal{P}^{\eta}_{k}\right).$ $$

We first evaluate the second probability on the right-hand side of (34).
Lemma 4 gives that

$$ $\displaystyle\mathbb{P}\left(\widehat{\mathcal{P}}^{\eta/2}_{k}\not\subseteq \mathcal{P}^{\eta}_{k}\right)\leq\left|\Theta\right|\Bigg{[}$ $\displaystyle\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{ \max}-\dfrac{3\eta}{4}\Big{\|}p_{k}^{\max}-\eta\right)\right)$ (35) $\displaystyle+$ $\displaystyle\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{ \max}-\dfrac{\eta}{4}\Big{\|}p_{k}^{\max}\right)\right)\Bigg{]}.$ $$

Next, by applying Lemma 5 with $\epsilon=\eta/2$, we can bound the first probability on the right-hand side of (34) as

$$ $\mathbb{P}\left(\hat{\theta}_{n}\notin\widehat{\mathcal{P}}^{\eta/2}_{k}\right )\leq\left|\Theta\right|\cdot\mathbb{E}\left[\exp\left(-\dfrac{B}{24}\cdot \dfrac{\eta^{2}}{\min\left(\hat{p}_{k}(\theta_{k}^{\max}),1-\hat{p}_{k}(\theta _{k}^{\max})\right)+\eta/2}\right)\right].$ (36) $$

Conditioned on the value of $\hat{p}_{k}(\theta_{k}^{\max})$, we can further upper-bound the right-hand side of (36) as follows

$$ $\displaystyle\mathbb{E}\left[\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{ \min\left(\hat{p}_{k}(\theta_{k}^{\max}),1-\hat{p}_{k}(\theta_{k}^{\max}) \right)+\eta/2}\right)\right]$ $\displaystyle\leq$ $\displaystyle\mathbb{P}\left(\hat{p}_{k}(\theta_{k}^{\max})\leq p_{k}^{\max}- \frac{\eta}{4}\right)\cdot\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{p_{k}^ {\max}+\eta/4}\right)+$ $\displaystyle\mathbb{P}\left(\left|\hat{p}_{k}(\theta_{k}^{\max})-p_{k}^{\max} \right|<\frac{\eta}{4}\right)\cdot\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2} }{\min\left(p_{k}^{\max},1-p_{k}^{\max}\right)+3\eta/4}\right)+$ $\displaystyle\mathbb{P}\left(\hat{p}_{k}(\theta_{k}^{\max})\geq p_{k}^{\max}+ \frac{\eta}{4}\right)\cdot\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{1-p_{k }^{\max}+\eta/4}\right)$ $\displaystyle\leq$ $\displaystyle\mathbb{P}\left(\hat{p}_{k}(\theta_{k}^{\max})\leq p_{k}^{\max}- \frac{\eta}{4}\right)+\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{\min\left( p_{k}^{\max},1-p_{k}^{\max}\right)+3\eta/4}\right)+$ $\displaystyle\mathbb{P}\left(\hat{p}_{k}(\theta_{k}^{\max})\geq p_{k}^{\max}+ \frac{\eta}{4}\right)\cdot\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{1-p_{k }^{\max}+\eta/4}\right)$ $\displaystyle\overset{(i)}{\leq}$ $\displaystyle\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{ \max}-\dfrac{\eta}{4}\Big{\|}p_{k}^{\max}\right)\right)+$ $\displaystyle\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{\min\left(p_{k}^{ \max},1-p_{k}^{\max}\right)+3\eta/4}\right)+$ $\displaystyle\mathbbm{1}\left(p_{k}^{\max}+\dfrac{\eta}{4}\leq 1\right)\cdot \exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{\max}+\dfrac{ \eta}{4}\Big{\|}p_{k}^{\max}\right)\right)\cdot\exp\left(-\dfrac{B}{24}\cdot \dfrac{\eta^{2}}{1-p_{k}^{\max}+\eta/4}\right)$ $$

where inequality $(i)$ results from applying Lemma 2 with $\hat{p}_{k}(\theta_{k}^{\max})$, the U-statistic estimate for $p_{k}^{\max}$.
Together, the above equations imply that

$$ $\displaystyle\mathbb{P}\left(\hat{\theta}_{n}\notin\Theta^{\delta}\right)$ $\displaystyle\leq$ $\displaystyle\left|\Theta\right|\Bigg{[}\exp\left(-\dfrac{n}{2k}\cdot D_{ \operatorname{KL}}\left(p_{k}^{\max}-\dfrac{3\eta}{4}\Big{\|}p_{k}^{\max}-\eta \right)\right)+$ $\displaystyle 2\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^ {\max}-\dfrac{\eta}{4}\Big{\|}p_{k}^{\max}\right)\right)+$ $\displaystyle\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{\min\left(p_{k}^{ \max},1-p_{k}^{\max}\right)+3\eta/4}\right)+$ $\displaystyle\mathbbm{1}\left(p_{k}^{\max}+\dfrac{\eta}{4}\leq 1\right)\cdot \exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{\max}+\dfrac{ \eta}{4}\Big{\|}p_{k}^{\max}\right)-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{1-p_{k} ^{\max}+\eta/4}\right)\Bigg{]}.$ $$

Since the above probability bound is left-continuous in $\eta$ and $\eta$ can be arbitrarily chosen from $(0,\;\bar{\eta}_{k,\delta})$, the validity of the case $\eta=\bar{\eta}_{k,\delta}$ follows from pushing $\eta$ to the limit $\bar{\eta}_{k,\delta}$. This gives (30).

To simplify the bound in the case $\bar{\eta}_{k,\delta}>4/5$. Consider the bound (30) with $\eta=\bar{\eta}_{k,\delta}$. Since $p_{k}^{\max}\geq\bar{\eta}_{k,\delta}$ by the definition of $\bar{\eta}_{k,\delta}$, it must hold that $p_{k}^{\max}+\bar{\eta}_{k,\delta}/4>4/5+1/5=1$, therefore the last term in the finite-sample bound (30) vanishes. To simplify the first two terms in the finite-sample bound, we note that

$$ $\displaystyle p_{k}^{\max}-\dfrac{3\bar{\eta}_{k,\delta}}{4}$ $\displaystyle\leq 1-\dfrac{3}{4}\cdot\dfrac{4}{5}=\dfrac{2}{5},$ $\displaystyle p_{k}^{\max}-\dfrac{3\bar{\eta}_{k,\delta}}{4}$ $\displaystyle\geq\bar{\eta}_{k,\delta}-\dfrac{3\bar{\eta}_{k,\delta}}{4}\geq \dfrac{1}{5},$ $\displaystyle p_{k}^{\max}-\dfrac{\bar{\eta}_{k,\delta}}{4}$ $\displaystyle\leq 1-\dfrac{1}{4}\cdot\dfrac{4}{5}=\dfrac{4}{5},$ $\displaystyle p_{k}^{\max}-\dfrac{\bar{\eta}_{k,\delta}}{4}$ $\displaystyle\geq\bar{\eta}_{k,\delta}-\dfrac{\bar{\eta}_{k,\delta}}{4}\geq \dfrac{3}{5},$ $$

and that $p_{k}^{\max}-\bar{\eta}_{k,\delta}\leq 1-\bar{\eta}_{k,\delta}\leq 1/5$, therefore by the bound (18) from Lemma 3, we can bound the first two terms as

$$ $\displaystyle\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{ \max}-\dfrac{3\bar{\eta}_{k,\delta}}{4}\Big{\|}p_{k}^{\max}-\bar{\eta}_{k, \delta}\right)\right)$ $\displaystyle\leq$ $\displaystyle\exp\left(\dfrac{n}{2k}\ln\left(2((p_{k}^{\max}-\bar{\eta}_{k, \delta})(1-p_{k}^{\max}+\bar{\eta}_{k,\delta}))^{1/5}\right)\right)$ $\displaystyle=$ $\displaystyle\left(2((p_{k}^{\max}-\bar{\eta}_{k,\delta})(1-p_{k}^{\max}+\bar{ \eta}_{k,\delta}))^{1/5}\right)^{n/(2k)}$ $\displaystyle\leq$ $\displaystyle\left(2(p_{k}^{\max}-\bar{\eta}_{k,\delta})^{1/5}\right)^{n/(2k)}$ $\displaystyle=$ $\displaystyle\left(2^{5}(p_{k}^{\max}-\bar{\eta}_{k,\delta})\right)^{n/(10k)},$ $$

and similarly

$$ $\displaystyle\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{ \max}-\dfrac{\bar{\eta}_{k,\delta}}{4}\Big{\|}p_{k}^{\max}\right)\right)$ $\displaystyle\leq$ $\displaystyle\exp\left(\dfrac{n}{2k}\ln\left(2(p_{k}^{\max}(1-p_{k}^{\max}))^{ 1/5}\right)\right)$ $\displaystyle=$ $\displaystyle\left(2(p_{k}^{\max}(1-p_{k}^{\max}))^{1/5}\right)^{n/(2k)}$ $\displaystyle\leq$ $\displaystyle\left(2(1-p_{k}^{\max})^{1/5}\right)^{n/(2k)}$ $\displaystyle=$ $\displaystyle\left(2^{5}(1-p_{k}^{\max})\right)^{n/(10k)}.$ $$

On the other hand, by Lemma 3 both $D_{\operatorname{KL}}\left(p_{k}^{\max}-3\bar{\eta}_{k,\delta}/4\|p_{k}^{\max}
-\bar{\eta}_{k,\delta}\right)$ and $D_{\operatorname{KL}}\left(p_{k}^{\max}-\bar{\eta}_{k,\delta}/4\|p_{k}^{\max}\right)$ are bounded below by $\bar{\eta}_{k,\delta}^{2}/8$, therefore

$$ $\exp\left(-\dfrac{n}{2k}\cdot D_{\operatorname{KL}}\left(p_{k}^{\max}-\dfrac{3 \bar{\eta}_{k,\delta}}{4}\Big{\|}p_{k}^{\max}-\bar{\eta}_{k,\delta}\right) \right)\leq\exp\left(-\dfrac{n}{2k}\cdot\dfrac{\bar{\eta}_{k,\delta}^{2}}{8} \right)\leq\exp\left(-\dfrac{n}{25k}\right),$ $$

and the same holds for $\exp\left(-n/(2k)\cdot D_{\operatorname{KL}}\left(p_{k}^{\max}-\bar{\eta}_{k,
\delta}/4\|p_{k}^{\max}\right)\right)$. For the third term in the bound (30) we have

$$ $\dfrac{\bar{\eta}_{k,\delta}^{2}}{\min(p_{k}^{\max},1-p_{k}^{\max})+3\bar{\eta }_{k,\delta}/4}\geq\dfrac{(4/5)^{2}}{\min(1,1/5)+3/4}\geq\dfrac{16}{25},$ $$

and hence

$$ $\exp\left(-\dfrac{B}{24}\cdot\dfrac{\bar{\eta}_{k,\delta}^{2}}{\min\left(p_{k} ^{\max},1-p_{k}^{\max}\right)+3\bar{\eta}_{k,\delta}/4}\right)\leq\exp\left(- \dfrac{B}{75/2}\right).$ $$

The first desired bound then follows by setting $C_{1},C_{2},C_{3}$ to be the appropriate constants. This completes the proof of Theorem 3.
$\Box$

Proof of Theorem 1. Algorithm 1 is a special case of Algorithm 3 with the learning algorithm (23) that outputs a singleton, therefore the results of Theorem 3 automatically apply. Since $\mathcal{E}_{k,\delta}=\sum_{\theta\in\Theta\backslash\Theta^{\delta}}p_{k}(
\theta)\geq\max_{\theta\in\Theta\backslash\Theta^{\delta}}p_{k}(\theta)$, it holds that $\eta_{k,\delta}\leq\bar{\eta}_{k,\delta}$. When $\eta_{k,\delta}>0$ we also have $\bar{\eta}_{k,\delta}>0$, therefore (8) follows from setting $\eta$ to be $\eta_{k,\delta}$ in (30), and (9) follows from upper bounding $\max_{\theta\in\Theta\backslash\Theta^{\delta}}p_{k}(\theta)$ with $\mathcal{E}_{k,\delta}$ in (31).
$\Box$

###### Theorem 3 (Finite-sample bound for Algorithm 3 )

###### Lemma 4

###### Lemma 5

### C.3 Proof for Theorem 2

We first present two lemmas to be used in the main proof. The first lemma characterizes the exponentially improving quality of the solution set retrieved in Phase I:

###### Lemma 6 (Quality of retrieved solutions in Algorithm 2 )

For every $k$ and $\delta\geq 0$, the set of retrieved solutions $\mathcal{S}$ from Phase I of Algorithm 2 with $k_{1}=k$ and without data splitting satisfies that

$$ $\mathbb{P}\left(\mathcal{S}\cap\Theta^{\delta}=\emptyset\right)\leq\min\left(e ^{-(1-\mathcal{E}_{k,\delta})/C_{4}},C_{5}\mathcal{E}_{k,\delta}\right)^{\frac {n}{C_{6}k}}+\exp\left(-\dfrac{B_{1}}{C_{7}}(1-\mathcal{E}_{k,\delta})\right),$ (37) $$

where $C_{4},C_{5},C_{6},C_{7}>0$ are universal constants. The same bound with $n$ replaced by $n/2$ holds true for Algorithm 2 with data splitting.

Proof of Lemma 6. Let $(z_{1}^{*},\ldots,z_{k}^{*})$ be a random subsample and $\mathbb{P}_{*}$ be the probability with respect to the subsampling randomness conditioned on the data and the algorithmic randomness. Consider the two probabilities

$$ $\mathbb{P}\left(\mathcal{A}(z_{1},\ldots,z_{k})\in\Theta^{\delta}\right),\ \mathbb{P}_{*}\left(\mathcal{A}(z_{1}^{*},\ldots,z_{k}^{*})\in\Theta^{\delta} \right).$ $$

We have $1-\mathcal{E}_{k,\delta}=\mathbb{P}\left(\mathcal{A}(z_{1},\ldots,z_{k})\in
\Theta^{\delta}\right)$,
and the conditional probability

$$ $\mathbb{P}\left(\mathcal{S}\cap\Theta^{\delta}=\emptyset\Big{\rvert}\mathbb{P} _{*}\left(\mathcal{A}(z_{1}^{*},\ldots,z_{k}^{*})\in\Theta^{\delta}\right) \right)=\left(1-\mathbb{P}_{*}\left(\mathcal{A}(z_{1}^{*},\ldots,z_{k}^{*})\in \Theta^{\delta}\right)\right)^{B_{1}}.$ $$

Therefore we can write

$$ $\displaystyle\mathbb{P}\left(\mathcal{S}\cap\Theta^{\delta}=\emptyset\right)$ $\displaystyle=$ $\displaystyle\mathbb{E}\left[\left(1-\mathbb{P}_{*}\left(\mathcal{A}(z_{1}^{*} ,\ldots,z_{k}^{*})\in\Theta^{\delta}\right)\right)^{B_{1}}\right]$ (38) $\displaystyle\leq$ $\displaystyle\mathbb{P}\left(\mathbb{P}_{*}\left(\mathcal{A}(z_{1}^{*},\ldots, z_{k}^{*})\in\Theta^{\delta}\right)<\frac{1-\mathcal{E}_{k,\delta}}{e}\right)+ \left(1-\frac{1-\mathcal{E}_{k,\delta}}{e}\right)^{B_{1}}$ $$

where $e$ is the base of the natural logarithm. Applying Lemma 2 with $\kappa(z_{1},\ldots,z_{k};\omega):=\mathbbm{1}\left(\mathcal{A}(z_{1},\ldots,z
_{k};\omega)\in\Theta^{\delta}\right)$ gives

$$ $\mathbb{P}\left(\mathbb{P}_{*}\left(\mathcal{A}(z_{1}^{*},\ldots,z_{k}^{*})\in \Theta^{\delta}\right)<\frac{1-\mathcal{E}_{k,\delta}}{e}\right)\leq\exp\left( -\frac{n}{2k}\cdot D_{\operatorname{KL}}\left(\frac{1-\mathcal{E}_{k,\delta}}{ e}\Big{\|}1-\mathcal{E}_{k,\delta}\right)\right).$ $$

Further applying the bound (17) from Lemma 3 to the KL divergence on the right-hand side leads to

$$ $D_{\operatorname{KL}}\left(\frac{1-\mathcal{E}_{k,\delta}}{e}\Big{\|}1- \mathcal{E}_{k,\delta}\right)\geq\frac{1-\mathcal{E}_{k,\delta}}{e}\ln\frac{1} {e}+1-\mathcal{E}_{k,\delta}-\frac{1-\mathcal{E}_{k,\delta}}{e}=\left(1-\frac{ 2}{e}\right)(1-\mathcal{E}_{k,\delta})$ $$

and

$$ $\displaystyle D_{\operatorname{KL}}\left(\frac{1-\mathcal{E}_{k,\delta}}{e} \Big{\|}1-\mathcal{E}_{k,\delta}\right)$ $\displaystyle=$ $\displaystyle D_{\operatorname{KL}}\left(1-\frac{1-\mathcal{E}_{k,\delta}}{e} \Big{\|}\mathcal{E}_{k,\delta}\right)$ $\displaystyle\geq$ $\displaystyle\left(1-\frac{1-\mathcal{E}_{k,\delta}}{e}\right)\ln\frac{1-(1- \mathcal{E}_{k,\delta})/e}{\mathcal{E}_{k,\delta}}-(1-\mathcal{E}_{k,\delta})+ \frac{1-\mathcal{E}_{k,\delta}}{e}\text{\ \ by bound (\ref{KL lower bound: ratio})}$ $\displaystyle\geq$ $\displaystyle\left(1-\frac{1-\mathcal{E}_{k,\delta}}{e}\right)\ln\left(1-\frac {1-\mathcal{E}_{k,\delta}}{e}\right)-\left(1-\frac{1}{e}\right)\ln\mathcal{E}_ {k,\delta}-1+\frac{1}{e}$ $\displaystyle\geq$ $\displaystyle\left(1-\frac{1}{e}\right)\ln\left(1-\frac{1}{e}\right)-\left(1- \frac{1}{e}\right)\ln\mathcal{E}_{k,\delta}-1+\frac{1}{e}$ $\displaystyle=$ $\displaystyle\left(1-\frac{1}{e}\right)\ln\frac{e-1}{e^{2}\mathcal{E}_{k, \delta}}.$ $$

Combining the two bounds for the KL divergence we have

$$ $\displaystyle\mathbb{P}\left(\mathbb{P}_{*}\left(\mathcal{A}(z_{1}^{*},\ldots, z_{k}^{*})\in\Theta^{\delta}\right)<\frac{1-\mathcal{E}_{k,\delta}}{e}\right)$ $\displaystyle\leq$ $\displaystyle\min\left(\exp\left(-\frac{n}{2k}\cdot\left(1-\frac{2}{e}\right)( 1-\mathcal{E}_{k,\delta})\right),\left(\frac{e^{2}\mathcal{E}_{k,\delta}}{e-1} \right)^{(1-1/e)\frac{n}{2k}}\right).$ $$

Note that the second term on the right-hand side of (38) satisfies that $\left(1-(1-\mathcal{E}_{k,\delta})/e\right)^{B_{1}}\leq\exp\left(-B_{1}(1-
\mathcal{E}_{k,\delta})/e\right)$.
Thus, we derive that

$$ $\displaystyle\mathbb{P}\left(\mathcal{S}\cap\Theta^{\delta}=\emptyset\right)$ $\displaystyle\leq$ $\displaystyle\min\left(\exp\left(-\frac{n}{2k}\cdot\left(1-\frac{2}{e}\right)( 1-\mathcal{E}_{k,\delta})\right),\left(\frac{e^{2}\mathcal{E}_{k,\delta}}{e-1} \right)^{(1-1/e)\frac{n}{2k}}\right)+\exp\left(-\frac{B_{1}(1-\mathcal{E}_{k, \delta})}{e}\right)$ $\displaystyle\leq$ $\displaystyle\min\left(\exp\left(-\dfrac{1-2/e}{1-1/e}\cdot(1-\mathcal{E}_{k, \delta})\right),\frac{e^{2}\mathcal{E}_{k,\delta}}{e-1}\right)^{(1-1/e)\frac{n }{2k}}+\exp\left(-\frac{B_{1}(1-\mathcal{E}_{k,\delta})}{e}\right).$ $$

The conclusion then follows by setting $C_{4},C_{5},C_{6},C_{7}$ to be the appropriate constants.
$\Box$

The second lemma gives bounds for the excess risk sensitivity $\bar{\eta}_{k,\delta}$ in the case of the set-valued learning algorithm (24):

###### Lemma 7 (Bounds of $\bar{\eta}_{k,\delta}$ for the set-valued learning algorithm ( 24 ))

Consider discrete decision space $\Theta$. If the set-valued learning algorithm

$$ $\mathbb{A}(z_{1},\ldots,z_{k};\omega):=\left\{\theta\in\Theta:\frac{1}{k}\sum_ {i=1}^{k}l(\theta,z_{i})\leq\min_{\theta^{\prime}\in\Theta}\frac{1}{k}\sum_{i= 1}^{k}l(\theta^{\prime},z_{i})+\epsilon\right\}$ $$

is used with $\epsilon\geq 0$, it holds that

$$ $\displaystyle p_{k}^{\max}=\max_{\theta\in\Theta}p_{k}(\theta)\geq 1-T_{k} \left(\dfrac{\epsilon}{2}\right),$ (39) $\displaystyle\max_{\theta\in\Theta\backslash\Theta^{\delta}}p_{k}(\theta)\leq T _{k}\left(\dfrac{\delta-\epsilon}{2}\right),$ (40) $$

and hence

$$ $\bar{\eta}_{k,\delta}\geq 1-T_{k}\left(\dfrac{\epsilon}{2}\right)-T_{k}\left( \dfrac{\delta-\epsilon}{2}\right),$ (41) $$

where $T_{k}$ is the tail probability defined in Theorem 2.

Proof of Lemma 7. Let $\hat{L}_{k}(\theta):=\frac{1}{k}\sum_{i=1}^{k}l(\theta,z_{i})$. Let $\theta^{*}$ be an optimal solution of (1). We have

$$ $\max_{\theta\in\Theta}p_{k}(\theta)\geq p_{k}(\theta^{*})=\mathbb{P}\left( \theta^{*}\in\widehat{\Theta}_{k}^{\epsilon}\right)\geq\mathbb{P}\left(\Theta^ {0}\subseteq\widehat{\Theta}_{k}^{\epsilon}\right).$ $$

To bound the probability on the right hand side, we write

$$ $\displaystyle\left\{\Theta^{0}\not\subseteq\widehat{\Theta}_{k}^{\epsilon}\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta\in\Theta^{0},\theta^{\prime}\in\Theta}\left\{\hat {L}_{k}(\theta)>\hat{L}_{k}(\theta^{\prime})+\epsilon\right\}$ $\displaystyle=$ $\displaystyle\bigcup_{\theta\in\Theta^{0},\theta^{\prime}\in\Theta}\left\{\hat {L}_{k}(\theta)-L(\theta)>\hat{L}_{k}(\theta^{\prime})-L(\theta^{\prime})+L( \theta^{\prime})-L(\theta)+\epsilon\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta\in\Theta^{0},\theta^{\prime}\in\Theta}\left\{\hat {L}_{k}(\theta)-L(\theta)>\hat{L}_{k}(\theta^{\prime})-L(\theta^{\prime})+ \epsilon\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta\in\Theta^{0},\theta^{\prime}\in\Theta}\left\{\hat {L}_{k}(\theta)-L(\theta)>\dfrac{\epsilon}{2}\text{ or }\hat{L}_{k}(\theta^{ \prime})-L(\theta^{\prime})<-\dfrac{\epsilon}{2}\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta\in\Theta}\left\{\left|\hat{L}_{k}(\theta)-L( \theta)\right|>\dfrac{\epsilon}{2}\right\}$ $\displaystyle=$ $\displaystyle\left\{\max_{\theta\in\Theta}\left|\hat{L}_{k}(\theta)-L(\theta) \right|>\dfrac{\epsilon}{2}\right\},$ $$

therefore

$$ $\max_{\theta\in\Theta}p_{k}(\theta)\geq\mathbb{P}\left(\max_{\theta\in\Theta} \left|\hat{L}_{k}(\theta)-L(\theta)\right|\leq\dfrac{\epsilon}{2}\right)\geq 1 -T_{k}\left(\dfrac{\epsilon}{2}\right).$ (42) $$

This proves (39). To bound the other term $\max_{\theta\in\Theta\backslash\Theta^{\delta}}p_{k}(\theta)$, for any $\theta\in\Theta\backslash\Theta^{\delta}$ it holds that

$$ $p_{k}(\theta)=\mathbb{P}\left(\theta\in\widehat{\Theta}_{k}^{\epsilon}\right) \leq\mathbb{P}\left(\widehat{\Theta}_{k}^{\epsilon}\not\subseteq\Theta^{\delta }\right),$ (43) $$

and hence $\max_{\theta\in\Theta\backslash\Theta^{\delta}}p_{k}(\theta)\leq\mathbb{P}
\left(\widehat{\Theta}_{k}^{\epsilon}\not\subseteq\Theta^{\delta}\right)$. To bound the latter, we have

$$ $\displaystyle\left\{\widehat{\Theta}_{k}^{\epsilon}\not\subseteq\Theta^{\delta }\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta,\theta^{\prime}\in\Theta\text{ s.t. }L(\theta^{ \prime})-L(\theta)>\delta}\left\{\hat{L}_{k}(\theta^{\prime})\leq\hat{L}_{k}( \theta)+\epsilon\right\}$ $\displaystyle=$ $\displaystyle\bigcup_{\theta,\theta^{\prime}\in\Theta\text{ s.t. }L(\theta^{ \prime})-L(\theta)>\delta}\left\{\hat{L}_{k}(\theta^{\prime})-L(\theta^{\prime })+L(\theta^{\prime})-L(\theta)\leq\hat{L}_{k}(\theta)-L(\theta)+\epsilon\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta,\theta^{\prime}\in\Theta\text{ s.t. }L(\theta^{ \prime})-L(\theta)>\delta}\left\{\hat{L}_{k}(\theta^{\prime})-L(\theta^{\prime })+\delta<\hat{L}_{k}(\theta)-L(\theta)+\epsilon\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta,\theta^{\prime}\in\Theta\text{ s.t. }L(\theta^{ \prime})-L(\theta)>\delta}\left\{\hat{L}_{k}(\theta^{\prime})-L(\theta^{\prime })<-\dfrac{\delta-\epsilon}{2}\text{ or }\hat{L}_{k}(\theta)-L(\theta)>\dfrac{ \delta-\epsilon}{2}\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta\in\Theta}\left\{\left|\hat{L}_{k}(\theta)-L( \theta)\right|>\dfrac{\delta-\epsilon}{2}\right\}$ $\displaystyle=$ $\displaystyle\left\{\max_{\theta\in\Theta}\left|\hat{L}_{k}(\theta)-L(\theta) \right|>\dfrac{\delta-\epsilon}{2}\right\},$ $$

therefore

$$ $\max_{\theta\in\Theta\backslash\Theta^{\delta}}p_{k}(\theta)\leq\mathbb{P} \left(\max_{\theta\in\Theta}\left|\hat{L}_{k}(\theta)-L(\theta)\right|>\dfrac{ \delta-\epsilon}{2}\right)\leq T_{k}\left(\dfrac{\delta-\epsilon}{2}\right).$ (44) $$

This immediately gives (40). (41) is obvious given (39) and (40).
$\Box$

To prove Theorem 2, we introduce some notation. For every non-empty subset $\mathcal{W}\subseteq\Theta$, we use the following counterpart of Definition 3. Let

$$ $\mathcal{W}^{\delta}:=\left\{\theta\in\mathcal{W}:L(\theta)\leq\min_{\theta^{ \prime}\in\mathcal{W}}L(\theta^{\prime})+\delta\right\}$ (45) $$

be the set of $\delta$-optimal solutions in the restricted decision space $\mathcal{W}$, and

$$ $\widehat{\mathcal{W}}_{k}^{\delta}:=\left\{\theta\in\mathcal{W}:\frac{1}{k} \sum_{i=1}^{k}l(\theta,z_{i})=\min_{\theta^{\prime}\in\mathcal{W}}\frac{1}{k} \sum_{i=1}^{k}l(\theta^{\prime},z_{i})+\delta\right\}$ (46) $$

be the set of $\delta$-optimal solutions of the SAA with an i.i.d. data set of size $k$.

Proof of Theorem 2 for $\mathsf{ROVEs}$. Given the retrieved solution set $\mathcal{S}$ and the chosen $\epsilon$, the rest of Phase II of Algorithm 2 exactly performs Algorithm 3 on the restricted problem $\min_{\theta\in\mathcal{S}}\mathbb{E}\left[l(\theta,z)\right]$ to obtain $\hat{\theta}_{n}$ with the data $\bm{z}_{\lfloor n/2\rfloor+1:n}$, the set-valued learning algorithm (24), the chosen $\epsilon$ value and $k=k_{2},B=B_{2}$.

To show the upper bound for the unconditional convergence probability $\mathbb{P}\left(\hat{\theta}_{n}\notin\Theta^{2\delta}\right)$, note that

$$ $\left\{\mathcal{S}\cap\Theta^{\delta}\neq\emptyset\right\}\cap\left\{L(\hat{ \theta}_{n})\leq\min_{\theta\in\mathcal{S}}L(\theta)+\delta\right\}\subseteq \left\{\hat{\theta}_{n}\in\Theta^{2\delta}\right\},$ $$

and hence by union bound we can write

$$ $\mathbb{P}\left(\hat{\theta}_{n}\notin\Theta^{2\delta}\right)\leq\mathbb{P} \left(\mathcal{S}\cap\Theta^{\delta}=\emptyset\right)+\mathbb{P}\left(L(\hat{ \theta}_{n})>\min_{\theta\in\mathcal{S}}L(\theta)+\delta\right).$ (47) $$

$\mathbb{P}\left(\mathcal{S}\cap\Theta^{\delta}=\emptyset\right)$ has a bound from Lemma 6. We focus on the second probability.

For a fixed retrieved subset $\mathcal{S}\subseteq\Theta$, define the tail of the maximum deviation on $\mathcal{S}$

$$ $T_{k}^{\mathcal{S}}(\cdot):=\mathbb{P}\left(\sup_{\theta\in\mathcal{S}}\left| \frac{1}{k}\sum_{i=1}^{k}l(\theta,z_{i})-L(\theta)\right|>\cdot\right).$ $$

It is straightforward that $T_{k}^{\mathcal{S}}(\cdot)\leq T_{k}(\cdot)$ where $T_{k}$ is the tail of the maximum deviation over the whole space $\Theta$. Since $\mathbb{P}\left(\epsilon\in[\underline{\epsilon},\overline{\epsilon}]\right)=1$, we have

$$ $1-T_{k_{2}}^{\mathcal{S}}\left(\dfrac{\epsilon}{2}\right)-T_{k_{2}}^{\mathcal{ S}}\left(\dfrac{\delta-\epsilon}{2}\right)\geq 1-T_{k_{2}}^{\mathcal{S}}\left( \dfrac{\underline{\epsilon}}{2}\right)-T_{k_{2}}^{\mathcal{S}}\left(\dfrac{ \delta-\overline{\epsilon}}{2}\right).$ $$

If $T_{k_{2}}\left((\delta-\overline{\epsilon})/2\right)+T_{k_{2}}\left(\underline
{\epsilon}/2\right)<1/5$, we have $T_{k_{2}}^{\mathcal{S}}\left((\delta-\overline{\epsilon})/2\right)+T_{k_{2}}^{
\mathcal{S}}\left(\underline{\epsilon}/2\right)<1/5$ and subsequently $1-T_{k_{2}}^{\mathcal{S}}\left((\delta-\epsilon)/2\right)-T_{k_{2}}^{\mathcal{
S}}\left(\epsilon/2\right)>4/5$, and hence $\bar{\eta}_{k_{2},\eta}\geq 1-T_{k_{2}}^{\mathcal{S}}\left((\delta-\epsilon)/2
\right)-T_{k_{2}}^{\mathcal{S}}\left(\epsilon/2\right)>4/5$ by Lemma 7 for Phase II of $\mathsf{ROVEs}$ conditioned on $\mathcal{S}$ and $\epsilon$, therefore the bound (31) from Theorem 3 applies. Using the inequalities (39) and (40) to upper bound the $\min(1-p_{k}^{\max},p_{k}^{\max}-\bar{\eta}_{k,\delta})$ term in (31) gives

$$ $\displaystyle\mathbb{P}\left(L(\hat{\theta}_{n})>\min_{\theta\in\mathcal{S}}L( \theta)+\delta\big{|}\mathcal{S},\epsilon\right)$ $\displaystyle\leq$ $\displaystyle\left|\mathcal{S}\right|\left(3\min\left(e^{-2/5},C_{1}\max\left( T_{k_{2}}^{\mathcal{S}}\left(\dfrac{\underline{\epsilon}}{2}\right),T_{k_{2}}^ {\mathcal{S}}\left(\dfrac{\delta-\overline{\epsilon}}{2}\right)\right)\right)^ {\frac{n}{2C_{2}k_{2}}}+\exp\left(-\dfrac{B_{2}}{C_{3}}\right)\right)$ $\displaystyle=$ $\displaystyle\left|\mathcal{S}\right|\left(3\min\left(e^{-2/5},C_{1}T_{k_{2}}^ {\mathcal{S}}\left(\dfrac{\min(\underline{\epsilon},\delta-\overline{\epsilon} )}{2}\right)\right)^{\frac{n}{2C_{2}k_{2}}}+\exp\left(-\dfrac{B_{2}}{C_{3}} \right)\right)$ $\displaystyle\leq$ $\displaystyle\left|\mathcal{S}\right|\left(3\min\left(e^{-2/5},C_{1}T_{k_{2}} \left(\dfrac{\min(\underline{\epsilon},\delta-\overline{\epsilon})}{2}\right) \right)^{\frac{n}{2C_{2}k_{2}}}+\exp\left(-\dfrac{B_{2}}{C_{3}}\right)\right).$ $$

Further relaxing $\left|\mathcal{S}\right|$ to $B_{1}$ and taking full expectation on both sides give

$$ $\mathbb{P}\left(L(\hat{\theta}_{n})>\min_{\theta\in\mathcal{S}}L(\theta)+ \delta\right)\leq B_{1}\left(3\min\left(e^{-2/5},C_{1}T_{k_{2}}\left(\dfrac{ \min(\underline{\epsilon},\delta-\overline{\epsilon})}{2}\right)\right)^{\frac {n}{2C_{2}k_{2}}}+\exp\left(-\dfrac{B_{2}}{C_{3}}\right)\right).$ $$

This leads to the desired bound (10) after the above bound is plugged into (47) and the bound (37) from Lemma 6 is applied with $k=k_{1}$.
$\Box$

Proof of Theorem 2 for $\mathsf{ROVE}$. For every non-empty subset $\mathcal{W}\subseteq\Theta$ and $k_{2}$, we consider the indicator

$$ $\mathbbm{1}_{k_{2}}^{\theta,\mathcal{W},\epsilon}(z_{1},\ldots,z_{k_{2}}):= \mathbbm{1}\left(\frac{1}{k_{2}}\sum_{i=1}^{k_{2}}l(\theta,z_{i})\leq\min_{ \theta^{\prime}\in\mathcal{W}}\frac{1}{k_{2}}\sum_{i=1}^{k_{2}}l(\theta^{ \prime},z_{i})+\epsilon\right)\quad\text{for }\theta\in\mathcal{W},\epsilon\in [0,\delta/2],$ $$

which indicates whether a solution $\theta\in\mathcal{W}$ is $\epsilon$-optimal for the SAA formed by $\left\{z_{1},\ldots,z_{k_{2}}\right\}$. Here we add $\epsilon$ and $\mathcal{W}$ to the superscript to emphasize its dependence on them. The counterparts of the solution probabilities $p_{k},\hat{p}_{k},\bar{p}_{k}$ for $\mathbbm{1}_{k_{2}}^{\theta,\mathcal{W},\epsilon}$ are

$$ $\displaystyle p_{k_{2}}^{\mathcal{W},\epsilon}(\theta)$ $\displaystyle:=\mathbb{E}\left[\mathbbm{1}_{k_{2}}^{\theta,\mathcal{W}, \epsilon}(z_{1},\ldots,z_{k_{2}})\right],$ $\displaystyle\hat{p}_{k_{2}}^{\mathcal{W},\epsilon}(\theta)$ $\displaystyle:=\mathbb{E}_{*}\left[\mathbbm{1}_{k_{2}}^{\theta,\mathcal{W}, \epsilon}(z_{1}^{*},\ldots,z_{k_{2}}^{*})\right],$ $\displaystyle\bar{p}_{k_{2}}^{\mathcal{W},\epsilon}(\theta)$ $\displaystyle:=\dfrac{1}{B_{2}}\sum_{b=1}^{B_{2}}\mathbbm{1}_{k_{2}}^{\theta, \mathcal{W},\epsilon}(z_{1}^{b},\ldots,z_{k_{2}}^{b}).$ $$

We need to show the uniform convergence of these probabilities for $\epsilon\in[0,\delta/2]$. To do so, we define a slighted modified version of $\mathbbm{1}_{k_{2}}^{\theta,\mathcal{W},\epsilon}$

$$ $\mathbbm{1}_{k_{2}}^{\theta,\mathcal{W},\epsilon-}(z_{1},\ldots,z_{k_{2}}):= \mathbbm{1}\left(\frac{1}{k_{2}}\sum_{i=1}^{k_{2}}l(\theta,z_{i})<\min_{\theta ^{\prime}\in\mathcal{W}}\frac{1}{k_{2}}\sum_{i=1}^{k_{2}}l(\theta^{\prime},z_{ i})+\epsilon\right)\quad\text{for }\theta\in\mathcal{W},\epsilon\in[0,\delta/2],$ $$

which indicates a strict $\epsilon$-optimal solution, and let $p_{k_{2}}^{\mathcal{W},\epsilon-},\hat{p}_{k_{2}}^{\mathcal{W},\epsilon-},\bar
{p}_{k_{2}}^{\mathcal{W},\epsilon-}$ be the corresponding counterparts of solution probabilities. For any integer $m>1$ we construct brackets of size at most $1/m$ to cover the family of indicator functions $\{\mathbbm{1}_{k_{2}}^{\theta,\mathcal{W},\epsilon}:\epsilon\in[0,\delta/2]\}$, i.e., let $m^{\prime}=\lfloor p_{k_{2}}^{\mathcal{W},\delta/2}(\theta)m\rfloor$ and

$$ $\displaystyle\epsilon_{0}$ $\displaystyle:=0,$ $\displaystyle\epsilon_{i}$ $\displaystyle:=\inf\left\{\epsilon\in[0,\delta/2]:p_{k_{2}}^{\mathcal{W}, \epsilon}(\theta)\geq i/m\right\}\quad\text{for }1\leq i\leq m^{\prime},$ $\displaystyle\epsilon_{m^{\prime}+1}$ $\displaystyle:=\dfrac{\delta}{2},$ $$

where we assume that $\epsilon_{i},i=0,\ldots,m^{\prime}+1$ are strictly increasing without loss of generality (otherwise we can delete duplicated values). Then for any $\epsilon\in[\epsilon_{i},\epsilon_{i+1})$, we have that

$$ $\displaystyle\bar{p}_{k_{2}}^{\mathcal{W},\epsilon}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon}(\theta)$ $\displaystyle\leq$ $\displaystyle\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i+1}-}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon_{i}}(\theta)$ $\displaystyle\leq$ $\displaystyle\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i+1}-}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon_{i+1}-}(\theta)+p_{k_{2}}^{\mathcal{W},\epsilon_{i+1}-}( \theta)-p_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)$ $\displaystyle\leq$ $\displaystyle\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i+1}-}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon_{i+1}-}(\theta)+\dfrac{1}{m}$ $$

and that

$$ $\displaystyle\bar{p}_{k_{2}}^{\mathcal{W},\epsilon}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon}(\theta)$ $\displaystyle\geq$ $\displaystyle\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon_{i+1}-}(\theta)$ $\displaystyle\geq$ $\displaystyle\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon_{i}}(\theta)+p_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta) -p_{k_{2}}^{\mathcal{W},\epsilon_{i+1}-}(\theta)$ $\displaystyle\geq$ $\displaystyle\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon_{i}}(\theta)-\dfrac{1}{m}.$ $$

Therefore

$$ $\displaystyle\sup_{\epsilon\in[0,\delta/2]}\left|\bar{p}_{k_{2}}^{\mathcal{W}, \epsilon}(\theta)-p_{k_{2}}^{\mathcal{W},\epsilon}(\theta)\right|$ (48) $\displaystyle\leq$ $\displaystyle\max_{0\leq i\leq m^{\prime}+1}\max\left(\left|\bar{p}_{k_{2}}^{ \mathcal{W},\epsilon_{i}}(\theta)-p_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta) \right|,\left|\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}-}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon_{i}-}(\theta)\right|\right)+\dfrac{1}{m}.$ $$

To show that the random variable in (48) converges to 0 0 in probability, we note that the U-statistic has the minimum variance among all unbiased estimators, in particular the following simple sample average estimators based on the first $\lfloor n/k_{2}\rfloor\cdot k_{2}$ data

$$ $\displaystyle\tilde{p}_{k_{2}}^{\mathcal{W},\epsilon}(\theta)$ $\displaystyle:=\dfrac{1}{\lfloor n/k_{2}\rfloor}\sum_{i=1}^{\lfloor n/k_{2} \rfloor}\mathbbm{1}_{k_{2}}^{\theta,\mathcal{W},\epsilon}(z_{k_{2}(i-1)+1}, \ldots,z_{k_{2}i}),$ $\displaystyle\tilde{p}_{k_{2}}^{\mathcal{W},\epsilon-}(\theta)$ $\displaystyle:=\dfrac{1}{\lfloor n/k_{2}\rfloor}\sum_{i=1}^{\lfloor n/k_{2} \rfloor}\mathbbm{1}_{k_{2}}^{\theta,\mathcal{W},\epsilon-}(z_{k_{2}(i-1)+1}, \ldots,z_{k_{2}i}).$ $$

Therefore we can write

$$ $\displaystyle\mathbb{E}\left[\left(\max_{0\leq i\leq m^{\prime}+1}\max\left( \left|\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)-p_{k_{2}}^{\mathcal{W },\epsilon_{i}}(\theta)\right|,\left|\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i} -}(\theta)-p_{k_{2}}^{\mathcal{W},\epsilon_{i}-}(\theta)\right|\right)\right)^ {2}\right]$ $\displaystyle\leq$ $\displaystyle\sum_{0\leq i\leq m^{\prime}+1}\left(\mathbb{E}\left[\left(\bar{p }_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)-p_{k_{2}}^{\mathcal{W},\epsilon_{ i}}(\theta)\right)^{2}\right]+\mathbb{E}\left[\left(\bar{p}_{k_{2}}^{\mathcal{ W},\epsilon_{i}-}(\theta)-p_{k_{2}}^{\mathcal{W},\epsilon_{i}-}(\theta)\right) ^{2}\right]\right)$ $\displaystyle\leq$ $\displaystyle\sum_{0\leq i\leq m^{\prime}+1}\left(\mathbb{E}\left[\left(\bar{p }_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)-\hat{p}_{k_{2}}^{\mathcal{W}, \epsilon_{i}}(\theta)\right)^{2}\right]+\mathbb{E}\left[\left(\hat{p}_{k_{2}}^ {\mathcal{W},\epsilon_{i}}(\theta)-p_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta )\right)^{2}\right]\right)+$ $\displaystyle\sum_{0\leq i\leq m^{\prime}+1}\left(\mathbb{E}\left[\left(\bar{p }_{k_{2}}^{\mathcal{W},\epsilon_{i}-}(\theta)-\hat{p}_{k_{2}}^{\mathcal{W}, \epsilon_{i}-}(\theta)\right)^{2}\right]+\mathbb{E}\left[\left(\hat{p}_{k_{2}} ^{\mathcal{W},\epsilon_{i}-}(\theta)-p_{k_{2}}^{\mathcal{W},\epsilon_{i}-}( \theta)\right)^{2}\right]\right)$ since $\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)$ and $\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}-}(\theta)$ are conditionally unbiased for $\hat{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)$ and $\hat{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}-}(\theta)$ $\displaystyle\leq$ $\displaystyle\sum_{0\leq i\leq m^{\prime}+1}\left(\mathbb{E}\left[\mathbb{E}_{ *}\left[\left(\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)-\hat{p}_{k_{2 }}^{\mathcal{W},\epsilon_{i}}(\theta)\right)^{2}\right]\right]+\mathbb{E}\left [\left(\tilde{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}}(\theta)-p_{k_{2}}^{ \mathcal{W},\epsilon_{i}}(\theta)\right)^{2}\right]\right)+$ $\displaystyle\sum_{0\leq i\leq m^{\prime}+1}\left(\mathbb{E}\left[\mathbb{E}_{ *}\left[\left(\bar{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}-}(\theta)-\hat{p}_{k_{ 2}}^{\mathcal{W},\epsilon_{i}-}(\theta)\right)^{2}\right]\right]+\mathbb{E} \left[\left(\tilde{p}_{k_{2}}^{\mathcal{W},\epsilon_{i}-}(\theta)-p_{k}^{ \mathcal{W},\epsilon_{i}-}(\theta)\right)^{2}\right]\right)$ $\displaystyle\leq$ $\displaystyle(m^{\prime}+2)\left(\dfrac{2}{B_{2}}+\dfrac{2}{\lfloor n/k_{2} \rfloor}\right)\leq(m+2)\left(\dfrac{2}{B_{2}}+\dfrac{4}{n/k_{2}}\right).$ $$

By Minkowski inequality, the supremum satisfies

$$ $\mathbb{E}\left[\sup_{\epsilon\in[0,\delta/2]}\left|\bar{p}_{k_{2}}^{\mathcal{ W},\epsilon}(\theta)-p_{k_{2}}^{\mathcal{W},\epsilon}(\theta)\right|\right] \leq\sqrt{(m+2)\left(\dfrac{2}{B_{2}}+\dfrac{4}{n/k_{2}}\right)}+\frac{1}{m}.$ $$

Choosing $m$ such that $m\to\infty$, $m/B_{2}\to 0$ and $mk_{2}/n\to 0$ leads to the convergence $\sup_{\epsilon\in[0,\delta/2]}\left|\bar{p}_{k_{2}}^{\mathcal{W},\epsilon}(
\theta)-p_{k_{2}}^{\mathcal{W},\epsilon}(\theta)\right|\to 0$ in probability. Since $\Theta$ has finite cardinality and has a finite number of subsets, it also holds that

$$ $\sup_{\mathcal{W}\subseteq\Theta,\theta\in\mathcal{W},\epsilon\in[0,\delta/2]} \left|\bar{p}_{k_{2}}^{\mathcal{W},\epsilon}(\theta)-p_{k_{2}}^{\mathcal{W}, \epsilon}(\theta)\right|\to 0\text{\ in probability}.$ (49) $$

Recall the bound (43) from the proof of Lemma 7. Here we have the similar bound $\max_{\theta\in\mathcal{W}\backslash\mathcal{W}^{\delta}}p_{k_{2}}^{\mathcal{W
},\epsilon}(\theta)\leq\mathbb{P}\left(\widehat{\mathcal{W}}_{k_{2}}^{\epsilon
}\not\subseteq\mathcal{W}^{\delta}\right)$, and hence

$$ $\sup_{\epsilon\in[0,\delta/2]}\max_{\theta\in\mathcal{W}\backslash\mathcal{W}^ {\delta}}p_{k}^{\mathcal{W},\epsilon}(\theta)\leq\sup_{\epsilon\in[0,\delta/2] }\mathbb{P}\left(\widehat{\mathcal{W}}_{k_{2}}^{\epsilon}\not\subseteq\mathcal {W}^{\delta}\right)=\mathbb{P}\left(\widehat{\mathcal{W}}_{k_{2}}^{\delta/2} \not\subseteq\mathcal{W}^{\delta}\right).$ $$

We bound the probability $\mathbb{P}\left(\widehat{\mathcal{W}}_{k_{2}}^{\delta/2}\not\subseteq\mathcal{
W}^{\delta}\right)$ more carefully. We let

$$ $\displaystyle\Delta_{o}$ $\displaystyle:=\min\left\{L(\theta^{\prime})-L(\theta):\theta,\theta^{\prime} \in\Theta,L(\theta^{\prime})>L(\theta)\right\}>0,$ $\displaystyle\hat{L}_{k_{2}}(\theta)$ $\displaystyle:=\frac{1}{k_{2}}\sum_{i=1}^{k_{2}}l(\theta,z_{i}),$ $$

and have

$$ $\displaystyle\left\{\widehat{\mathcal{W}}_{k_{2}}^{\delta/2}\not\subseteq \mathcal{W}^{\delta}\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta,\theta^{\prime}\in\mathcal{W}\text{ s.t. }L( \theta^{\prime})-L(\theta)>\delta}\left\{\hat{L}_{k_{2}}(\theta^{\prime})\leq \hat{L}_{k_{2}}(\theta)+\dfrac{\delta}{2}\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta,\theta^{\prime}\in\Theta\text{ s.t. }L(\theta^{ \prime})-L(\theta)>\delta}\left\{\hat{L}_{k_{2}}(\theta^{\prime})-L(\theta^{ \prime})+L(\theta^{\prime})-L(\theta)\leq\hat{L}_{k_{2}}(\theta)-L(\theta)+ \dfrac{\delta}{2}\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta,\theta^{\prime}\in\Theta\text{ s.t. }L(\theta^{ \prime})-L(\theta)>\delta}\left\{\hat{L}_{k_{2}}(\theta^{\prime})-L(\theta^{ \prime})+\max(\Delta,\delta)\leq\hat{L}_{k_{2}}(\theta)-L(\theta)+\dfrac{ \delta}{2}\right\}$ by the definition of $\Delta_{o}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta,\theta^{\prime}\in\Theta}\left\{\hat{L}_{k_{2}}( \theta^{\prime})-L(\theta^{\prime})+\max\left(\Delta_{o}-\dfrac{\delta}{2}, \dfrac{\delta}{2}\right)\leq\hat{L}_{k_{2}}(\theta)-L(\theta)\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta,\theta^{\prime}\in\Theta}\left\{\hat{L}_{k_{2}}( \theta^{\prime})-L(\theta^{\prime})\leq-\max\left(\dfrac{\Delta_{o}}{2}-\dfrac {\delta}{4},\dfrac{\delta}{4}\right)\text{ or }\hat{L}_{k_{2}}(\theta)-L( \theta)\geq\max\left(\dfrac{\Delta_{o}}{2}-\dfrac{\delta}{4},\dfrac{\delta}{4} \right)\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta\in\Theta}\left\{\left|\hat{L}_{k_{2}}(\theta)-L( \theta)\right|\geq\max\left(\dfrac{\Delta_{o}}{2}-\dfrac{\delta}{4},\dfrac{ \delta}{4}\right)\right\}$ $\displaystyle\subseteq$ $\displaystyle\bigcup_{\theta\in\Theta}\left\{\left|\hat{L}_{k_{2}}(\theta)-L( \theta)\right|\geq\dfrac{\Delta_{o}}{4}\right\}$ $\displaystyle\subseteq$ $\displaystyle\left\{\sup_{\theta\in\Theta}\left|\hat{L}_{k_{2}}(\theta)-L( \theta)\right|\geq\dfrac{\Delta_{o}}{4}\right\},$ $$

where the last line holds because $\max\left(\Delta_{o}/2-\delta/4,\delta/4\right)\geq\Delta_{o}/4$. This gives

$$ $\sup_{\epsilon\in[0,\delta/2]}\max_{\theta\in\mathcal{W}\backslash\mathcal{W}^ {\delta}}p_{k_{2}}^{\mathcal{W},\epsilon}(\theta)\leq T_{k_{2}}\left(\dfrac{ \Delta_{o}}{4}\right)\to 0\text{\ \ as\ }k_{2}\to\infty.$ $$

We also have the trivial bound $\inf_{\epsilon\in[0,\delta/2]}\max_{\theta\in\mathcal{W}}p_{k_{2}}^{\mathcal{W
},\epsilon}(\theta)=\max_{\theta\in\mathcal{W}}p_{k_{2}}^{\mathcal{W},0}(
\theta)\geq 1/\left|\mathcal{W}\right|$, where the inequality comes from the fact that $\sum_{\theta\in\mathcal{W}}p_{k_{2}}^{\mathcal{W},0}(\theta)\geq 1$. Now choose a $\underline{k}<\infty$ such that

$$ $T_{k_{2}}\left(\dfrac{\Delta_{o}}{4}\right)\leq\frac{1}{2\left|\Theta\right|} \text{\ \ for all $k_{2}\geq\underline{k}$}$ $$

and we have for all $k_{2}\geq\underline{k}$ and all non-empty $\mathcal{W}\subseteq\Theta$ that

$$ $\displaystyle\inf_{\epsilon\in[0,\delta/2]}\left(\max_{\theta\in\mathcal{W}}p_ {k_{2}}^{\mathcal{W},\epsilon}(\theta)-\max_{\theta\in\mathcal{W}\backslash \mathcal{W}^{\delta}}p_{k_{2}}^{\mathcal{W},\epsilon}(\theta)\right)$ $\displaystyle\geq$ $\displaystyle\inf_{\epsilon\in[0,\delta/2]}\max_{\theta\in\mathcal{W}}p_{k_{2} }^{\mathcal{W},\epsilon}(\theta)-\sup_{\epsilon\in[0,\delta/2]}\max_{\theta\in \mathcal{W}\backslash\mathcal{W}^{\delta}}p_{k_{2}}^{\mathcal{W},\epsilon}(\theta)$ $\displaystyle\geq$ $\displaystyle\dfrac{1}{\left|\mathcal{W}\right|}-\frac{1}{2\left|\Theta\right| }\geq\dfrac{1}{2\left|\Theta\right|}.$ $$

Due to the uniform convergence (49), we have

$$ $\min_{\mathcal{W}\subseteq\Theta}\inf_{\epsilon\in[0,\delta/2]}\left(\max_{ \theta\in\mathcal{W}}\bar{p}_{k_{2}}^{\mathcal{W},\epsilon}(\theta)-\max_{ \theta\in\mathcal{W}\backslash\mathcal{W}^{\delta}}\bar{p}_{k_{2}}^{\mathcal{W },\epsilon}(\theta)\right)\to\min_{\mathcal{W}\subseteq\Theta}\inf_{\epsilon \in[0,\delta/2]}\left(\max_{\theta\in\mathcal{W}}p_{k_{2}}^{\mathcal{W}, \epsilon}(\theta)-\max_{\theta\in\mathcal{W}\backslash\mathcal{W}^{\delta}}p_{ k_{2}}^{\mathcal{W},\epsilon}(\theta)\right)$ $$

in probability, and hence

$$ $\mathbb{P}\left(\min_{\mathcal{W}\subseteq\Theta}\inf_{\epsilon\in[0,\delta/2] }\left(\max_{\theta\in\mathcal{W}}\bar{p}_{k_{2}}^{\mathcal{W},\epsilon}( \theta)-\max_{\theta\in\mathcal{W}\backslash\mathcal{W}^{\delta}}\bar{p}_{k_{2 }}^{\mathcal{W},\epsilon}(\theta)\right)\leq 0\right)\to 0.$ (50) $$

Finally, we combine all the pieces to get

$$ $\displaystyle\left\{\hat{\theta}_{n}\not\in\Theta^{2\delta}\right\}$ $\displaystyle\subseteq$ $\displaystyle\left\{\mathcal{S}\cap\Theta^{\delta}=\emptyset\right\}\cup\left \{\hat{\theta}_{n}\not\in\mathcal{S}^{\delta}\right\}$ $\displaystyle\subseteq$ $\displaystyle\left\{\mathcal{S}\cap\Theta^{\delta}=\emptyset\right\}\cup\left \{\max_{\theta\in\mathcal{S}}\bar{p}_{k_{2}}^{\mathcal{S},\epsilon}(\theta)- \max_{\theta\in\mathcal{S}\backslash\mathcal{S}^{\delta}}\bar{p}_{k_{2}}^{ \mathcal{S},\epsilon}(\theta)\leq 0\right\}$ $\displaystyle\subseteq$ $\displaystyle\left\{\mathcal{S}\cap\Theta^{\delta}=\emptyset\right\}\cup\left \{\epsilon>\dfrac{\delta}{2}\right\}\cup\left\{\inf_{\epsilon\in[0,\delta/2]} \left(\max_{\theta\in\mathcal{S}}\bar{p}_{k_{2}}^{\mathcal{S},\epsilon}(\theta )-\max_{\theta\in\mathcal{S}\backslash\mathcal{S}^{\delta}}\bar{p}_{k_{2}}^{ \mathcal{S},\epsilon}(\theta)\right)\leq 0\right\}$ $\displaystyle\subseteq$ $\displaystyle\left\{\mathcal{S}\cap\Theta^{\delta}=\emptyset\right\}\cup\left \{\epsilon>\dfrac{\delta}{2}\right\}\cup\left\{\min_{\mathcal{W}\subseteq \Theta}\inf_{\epsilon\in[0,\delta/2]}\left(\max_{\theta\in\mathcal{W}}\bar{p}_ {k_{2}}^{\mathcal{W},\epsilon}(\theta)-\max_{\theta\in\mathcal{W}\backslash \mathcal{W}^{\delta}}\bar{p}_{k_{2}}^{\mathcal{W},\epsilon}(\theta)\right)\leq 0 \right\}.$ $$

By Lemma 6 we have $\mathbb{P}\left(\mathcal{S}\cap\Theta^{\delta}=\emptyset\right)\to 0$ under the conditions that $\limsup_{k\to\infty}\mathcal{E}_{k,\delta}<1$ and $k_{1},n/k_{1},B_{1}\to\infty$. Together with the condition $\mathbb{P}\left(\epsilon\geq\delta/2\right)\to 0$ and (50), we conclude $\mathbb{P}\left(\hat{\theta}_{n}\not\in\Theta^{2\delta}\right)\to 0$ by the union bound.
$\Box$

###### Lemma 6 (Quality of retrieved solutions in Algorithm 2 )

###### Lemma 7 (Bounds of η ¯ k , δ subscript ¯ 𝜂 𝑘 𝛿 \bar{\eta}_{k,\delta} over¯ start_ARG italic_η end_ARG start_POSTSUBSCRIPT italic_k , italic_δ end_POSTSUBSCRIPT for the set-valued learning algorithm ( 24 ))

## Appendix D Supplementary material for numerical experiments

This section supplements Section 3. We first provide details for the architecture of the neural networks in Section D.1, and the considered stochastic programs in Section D.2. Section D.3 presents a comprehensive profiling of hyperparameters of our methods, and Section D.4 provides additional experimental results.

### D.1 MLP architecture

The input layer of our MLPs has the same number of neurons as the input dimension, and the output layer is a single neuron that gives the final prediction. All activations are ReLU. For experiments on synthetic data, the architecture of hidden layers is as follows under different numbers of hidden layers $H$:

- •
$H=2$: Each hidden layer has 50 neurons.
- •
$H=4$: There are 50, 300, 300, 50 neurons from the first to the fourth hidden layer.
- •
$H=6$: There are 50, 300, 500, 500 300, 50 neurons from the first to the sixth hidden layer.
- •
$H=8$: There are 50, 300, 500, 800, 800 500 300, 50 neurons from the first to the eighth hidden layer.

For experiments on real data, the MLP with 4 hidden layers has 100, 300, 300, 100 neurons from the first to the fourth hidden layer.

### D.2 Stochastic programming problems

Resource allocation (Kleywegt et al., )   The decision maker wants to choose a subset of $m$ projects.
A quantity $q$ of low-cost resource is available to be allocated, and any additional resource can be obtained at an incremental unit cost $c$.
Each project $i$ has an expected reward $r_{i}$.
The amount of resource required by each project $i$ is a random variable, denoted by $W_{i}$.
We can formulate the problem as

$$ $\max_{\theta\in\{0,1\}^{m}}\left\{\sum_{i=1}^{m}r_{i}\theta_{i}-c\mathbb{E} \left[\sum_{i=1}^{m}W_{i}\theta_{i}-q\right]^{+}\right\}.$ (51) $$

In the experiment, we consider the three-product scenario, i.e., $m=3$, and assume that the random variable $W_{i}$ follows the Pareto distribution.

Supply chain network design (Shapiro et al., , Chapter 1.5)
Consider a network of suppliers, processing facilities, and customers, where the goal is to optimize the overall supply chain efficiency.
The supply chain design problem can be formulated as a two-stage stochastic optimization problem

$$ $\min_{\theta\in\{0,1\}^{\left|P\right|}}\sum_{p\in P}c_{p}\theta_{p}+\mathbb{E }[Q(\theta,z)],$ (52) $$

where $P$ is the set of processing facilities, $c_{p}$ is the cost of opening facility $p$, and $z$ is the vector of (random) parameters, i.e., $(h,q,d,s,R,M)$ in (53).
Function $Q(\theta,z)$ represents the total processing and transportation cost, and it is equal to the optimal objective value of the following second-stage problem:

$$ $\begin{array}[]{cl}\min_{y\geq 0,z\geq 0}&q^{\top}y+h^{\top}z\[5.0pt] \text{ s.t. }&Ny=0,\[5.0pt] &Cy+z\geq d,\[5.0pt] &Sy\leq s,\[5.0pt] &Ry\leq M\theta,\end{array}$ (53) $$

where $N,C,S$ are appropriate matrices that describe the network flow constraints.
More details about this example can be found in (Shapiro et al., , Chapter 1.5).
In our experiment, we consider the scenario of 3 suppliers, 2 facilities, 3 consumers, and 5 products.
We choose supply $s$ and demand $d$ as random variables that follow the Pareto distribution.

Maximum weight matching and stochastic linear program
We explore both the maximum weight matching problem and the linear program that arises from it.
Let $G=(V,E)$ be a general graph, where each edge $e\in E$ is associated with a (possibly) random weight $w_{e}$. For each node $v\in V$, denote $E(v)$ as the set of edges incident to $v$.
Based on this setup, we consider the following linear program

$$ $\begin{array}[]{lll}\max_{\theta\in[0,1]^{\left|E\right|}}&\mathbb{E}\left[ \sum_{e\in E}w_{e}\theta_{e}\right]&\[5.0pt] \text{subject to }&\sum_{e\in E(v)}a_{e}\theta_{e}\leq 1,&\forall v\in V,\end{array}$ (54) $$

where $a_{e}$ is some positive coefficient.
When $a_{e}=1$ for all $e\in E$ and $\theta$ is restricted to the discrete set $\{0,1\}^{\left|E\right|}$, (54) is equivalent to the maximum weight matching problem.
For the maximum weight matching, we consider a complete bipartite graph with 5 nodes on each side (the dimension is 25).
The weights of nine edges are Pareto distributed and the remaining are prespecified constants.
For the linear programming problem, we consider a 28-dimensional instance (the underlying graph is an 8-node complete graph), where all $w_{e}$ follows the Pareto distribution.

Mean-variance portfolio optimization
Consider constructing a portfolio based on $m$ assets. Each asset $i$ has a rate of return $r_{i}$ that is random with mean $\mu_{i}$. The goal is to minimize the variance of the portfolio while ensuring that the expected rate of return surpasses a target level $b$.
The problem can be formulated as

$$ $\begin{array}[]{ll}\min_{\theta}&\mathbb{E}\left[(\sum_{i=1}^{m}(r_{i}-\mu_{i} )\theta_{i})^{2}\right]\[5.0pt] \text{subject to }&\sum_{i=1}^{m}\mu_{i}\theta_{i}\geq b,\[5.0pt] &\sum_{i=1}^{m}\theta_{i}=1,\[5.0pt] &\theta_{i}\geq 0\ \forall i=1,\ldots,m\end{array}$ (55) $$

where $\theta$ is the decision variable and each $\mu_{i}$ is assumed known.
In the experiment, we consider a scenario with 10 assets, i.e., $m=10$, where each rate of return $r_{i}$ is a linear combination of the rates of return of 100 underlying assets in the form $r_{i}=\tilde{r}_{10(i-1)+1}/2+\sum_{j=1}^{100}\tilde{r}_{j}/200$. Each of these underlying assets has a Pareto rate of return $\tilde{r}_{j},j=1,\ldots,100$.

### D.3 Hyperparameter profiling

We test the effect of different hyperparameters in our ensemble methods, including subsample sizes $k,k_{1},k_{2}$, ensemble sizes $B,B_{1},B_{2}$, and threshold $\epsilon$.
Throughout this profiling stage, we use the sample average approximation (SAA) as the base algorithm.
To profile the effect of subsample sizes and ensemble sizes, we consider the resource allocation problem.

Figure: (a) Profiling for $k$ ($\mathsf{MoVE}$).
Refer to caption: extracted/5899831/img/profiling_alg1_k_2.png

Figure: (a) Profiling for $k$ (instance 1).
Refer to caption: extracted/5899831/img/profiling_alg1_k_1.png

Figure: (a) Figure 5(a): suboptimal probability.
Refer to caption: extracted/5899831/img/SSKP_prob_suboptimal_2.png

Figure: (a) $\mathsf{MoVE}$.
Refer to caption: extracted/5899831/img/profiling_alg1_B.png

Figure: (a) Multiple optima.
Refer to caption: extracted/5899831/img/profiling_alg4_epsilon_multiple.png

Subsample size
We explored scenarios where $k$ (equivalently $k_{1}$ and $k_{2}$) is both dependent on and independent of the total sample size $n$ (see Figures 5(a), 6(a), and 6(b)).
The results suggest that a constant $k$ generally suffices, although the optimal $k$ varies by problem instance.
For example, Figures 6(a) and 6(b) show that $k=2$ yields the best performance; increasing $k$ degrades results.
Conversely, in Figure 5(a), $k=2$ proves inadequate, with larger $k$ delivering good results.
The underlying reason is that the effective performance of $\mathsf{MoVE}$ requires $\theta^{*}\in\operatorname*{arg\,max}_{\theta\in\Theta}p_{k}(\theta)$.
In the former, this is achieved with only two samples, enabling $\mathsf{MoVE}$ to identify $\theta^{*}$ with a subsample size of $2$.
For the latter, a higher number of samples is required to meet this condition, explaining the suboptimal performance at $k=2$.
In Figure 7, we simulate $p_{k}(\theta)$ for the two cases, which further explains the influence of the subsample size.

Ensemble size
In Figure 8, we illustrate the performance of $\mathsf{MoVE}$ and $\mathsf{ROVE}$ under different $B,B_{1},B_{2}$, where we set $k=k_{1}=k_{2}=10$ and $\epsilon=0.005$.
From the figure, we find that the performance of our ensemble methods is improving in the ensemble sizes.

Threshold $\epsilon$
The optimal choice of $\epsilon$ in $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ is problem-dependent and related to the number of (near) optimal solutions.
This dependence is illustrated by the performance of $\mathsf{ROVE}$ shown in Figures 5(b) and 5(c).
Hence, we propose an adaptive strategy defined as follows: Let $g(\epsilon):=1/B_{2}\cdot\sum_{b=1}^{B_{2}}\mathbbm{1}(\hat{\theta}_{n}(
\epsilon)\in\widehat{\Theta}^{\epsilon,b}_{k_{2}})$, where we use $\hat{\theta}_{n}(\epsilon)$ to emphasize the dependency of $\hat{\theta}_{n}$ on $\epsilon$.
Then, we select $\epsilon^{*}:=\min\left\{\epsilon:g(\epsilon)\geq 1/2\right\}$.
By definition, $g(\epsilon)$ is the proportion of times that $\hat{\theta}_{n}(\epsilon)$ is included in the “near optimum set” $\widehat{\Theta}^{\epsilon,b}_{k_{2}}$.
The choice of $\epsilon^{*}$ makes it more likely for the true optimal solution to be included in the “near optimum set”, instead of being ruled out by suboptimal solutions.
Practically, $\epsilon^{*}$ can be efficiently determined using a binary search as an intermediate step between Phases I and II.
To prevent data leakage, we compute $\epsilon^{*}$ using $\mathbf{z}_{1:\lfloor\frac{n}{2}\rfloor}$ (Phase I data) for $\mathsf{ROVEs}$.
From Figure 5, we observe that this adaptive strategy exhibits decent performance for all scenarios.
Similar behaviors can also be observed for $\mathsf{ROVEs}$ in Figure 9.

### D.4 Additional experimental results

Here, we present additional figures that supplement the experiments and discussions in Section 3.
Recall that $\mathsf{MoVE}$ refers to Algorithm 1, $\mathsf{ROVE}$ refers to Algorithm 2 without data splitting, and $\mathsf{ROVEs}$ refers to Algorithm 2 with data splitting.
We briefly introduce each figure below and refer the reader to the figure caption for detailed discussions. Figures 10-16 all follow the recommended configuration listed in Section 3.

- •
Figure 10 supplements the results in Figure 1 with MLPs with $H=2,4$ hidden layers.
- •
Figure 11 shows results for MLP regression on a slightly different synthetic example than in Section 3.1.
- •
Figures 12 and 13 show results for regression on synthetic data with least squares regression and Ridge regression as the base learning algorithms respectively.
- •
Figure 14 shows results on the stochastic linear program example with light-tailed uncertainties.
- •
Figure 15 contains additional results on the supply chain network design example for different choices of hyperparameters and a different problem instance with strong correlation between solutions.
- •
In Figure 16, we apply our ensemble methods to resource allocation and maximum weight matching using DRO with Wasserstein metric as the base algorithm. This result, together with Figure 3 where the base algorithm is SAA, demonstrates that the benefit of our ensemble methods is agnostic to the underlying base algorithm.
- •
In Figure 17, we simulate the generalization sensitivity $\bar{\eta}_{k,\delta}$, defined in (29), which explains the superior performance of $\mathsf{ROVE}$ and $\mathsf{ROVEs}$ in the presence of multiple optimal solutions.

Figure: (a) Pareto noise, $H=2$.
Refer to caption: extracted/5899831/img/NN_d50_l2_ES_plots/NN_d50_l2_ES_avg_ylog.png

Figure: (a) Average performance with $95\%$ CIs.
Refer to caption: extracted/5899831/img/NN_d30_l4_ES_WB_plots/NN_d30_l4_ES_WB_avg_ylog.png

Figure: (a) $d=10$.
Refer to caption: extracted/5899831/img/LR_d10_plots/LR_d10_avg.png

Figure: (a) $d=10$.
Refer to caption: extracted/5899831/img/Ridge_d10_plots/Ridge_d10_avg.png

Figure: (a) Instance 1.
Refer to caption: extracted/5899831/img/SAA_comparison_LP_light_1.png

Figure: (a) $k=k_{1}=k_{2}=\max(10,n/10)$.
Refer to caption: extracted/5899831/img/OptCIPlots/NetworkHeavy3Data_avg.png

Figure: (a) Resource allocation.
Refer to caption: extracted/5899831/img/OptCIPlots/SSKP_DRO_error.png

Figure: Figure 17: Comparison of $\bar{\eta}_{k,\delta}$ for $\mathsf{MoVE}$ and $\mathsf{ROVE}$ in a linear program with multiple optima (corresponds to the instance in Figure 3(e)). Threshold $\epsilon$ is chosen as $\epsilon=4$ when $k=k_{1}=k_{2}=10$ and $\epsilon=2.5$ when $k=k_{1}=k_{2}=50$, according to the adaptive strategy. Note that $\bar{\eta}_{k,\delta}=\max_{\theta\in\Theta}p_{k}(\theta)-\max_{\theta\in \Theta\backslash\Theta^{\delta}}p_{k}(\theta)$ by (29), which measures the generalization sensitivity. For $\mathsf{MoVE}$, we have $p_{k}(\theta)=\mathbb{P}(\hat{\theta}_{k}^{SAA}=\theta)$; and for $\mathsf{ROVE}$, we have $p_{k}(\theta)=\mathbb{P}(\theta\in\widehat{\Theta}_{k}^{\epsilon})$, where $\widehat{\Theta}_{k}^{\epsilon}$ is the $\epsilon$-optimal set of SAA defined in (28). From the figure, we can observe that the issue brought by the presence of multiple optimal solutions can be alleviated using the two-phase strategy in $\mathsf{ROVE}$.
Refer to caption: extracted/5899831/img/eta_comparison_LP.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 1
questions: - What happens in the case that the data is not heavy-tailed? Do the gains from MoVE and ROVE become constant factor as before? Can you provide the bounds and the correct selection of k for this case?
- It appears that for the method to work, k (sub-sampled data) has to be selected at a lower rate than n which implies that the ensemble size (B = O(n/k)) will be very large. How do existing generalization bounds for ensembles scale with the ensemble size and what are the consequences for them when B is not constant but rather a function of n?
- It seems like in the real-world dataset settings, ROVE is underperforming the base method. Can you detail out why this is the case and what assumptions seem to be breaking?
- Would it be possible to test out this method with slightly more real world datasets (like MNIST, CIFAR and converting them to a regression problem)? I would like to see the convergence rate of ROVE when compared with traditional ensembles as a function of the number of learners in the ensemble.
- Theorems 1 and 2 are quite dense and hard to parse. Would it be possible to instantiate both the theorems for a few example function classes and loss functions and have them as corollaries?
- In the discrete setting, what is the dependence of the generalization error on the size of the parameter class for the traditional ensembling methods?
soundness: 2
strengths: - Provide a new approach to ensemble learning with exponentially decaying tails for excess risk. 
- The paper presents experimental evaluations to support the theoretical claims.
summary: This paper presents a novel approach to ensemble learning called MoVE (Majority Vote Ensembling) and ROVE (Retrieval and ϵ-Optimality Vote Ensembling). The approach leverages subsampling to create ensembles that can produce exponentially decaying tails for excess risk, overcoming the polynomial decay limitations typically faced in heavy-tailed distributions. By iteratively training on subsets of the data and combining models through majority voting or likelihood-based selection, the method aims to provide stronger generalization than traditional ensembling techniques, especially in contexts with heavy-tailed data distributions. Empirical evaluations demonstrate improved out-of-sample performance across synthetic and real datasets and various stochastic programming problems.
weaknesses: - The theoretical results are quite dense and hard to parse. Need to provide example instantiations to allow the reader to get a better sense of them.
- The proposed algorithms look quite computationally expensive and yet underperform the baseline on real-world tasks.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: **1)** What does "majority vote" mean for the models? What does it imply in practice? How can models trained on different subsets of data end up producing the same outputs that allow for voting, as seen in the final lines of Alg-1 and Alg-2? In standard ensemble methods, majority voting is based on the predictions of the models. I might not be familiar with this specific approach to ensembling, so I would appreciate a clearer explanation.

**2)** In Algorithm 2 line 231-233, how to get $\theta_{\prime}$ in practice?

**3)** For more questions, see Weaknesses.
soundness: 3
strengths: The paper presents a fresh view on ensemble learning, focusing on exponential improvement in generalization tail bounds rather than the traditional variance reduction approach. 

The paper offers valuable theoretical contributions, e.g. developing a sharper concentration result for U-statistics with binary kernels, learners, performing a sensitivity analysis on the regret, and developing a uniform law of large numbers (LLN) for the class of events of being $\epsilon$-optimal. 

Theoretical insights are complemented by experiments on both real and synthetic data using neural networks and stochastic programs, demonstrating relevance of the theory.
summary: This paper presents an ensemble learning technique that enhances the generalization performance by selecting the best models trained on subsamples via majority voting, resulting in exponentially decaying tails ($C_{2}\gamma^{n/k}$) for excess risk, even when base laerners have slow (polynomial i.e. $C_{1}n^{-\alpha}$) decay rates. The method is agnostic to the underlying base learner  and offers a stronger improvement than traditional variance reduction, providing rate improvement rather than just a constant factor improvement. Experimetns show that the method also effective in scenarios with heavy-tailed data and slow convergence rates.
weaknesses: **Only one baseline used**: The paper could benefit from including additional baselines beyond just the 'base' model. Since the authors mention variance reduction, it would be useful to incorporate it as a baseline, as well as to provide more detailed comparisons with state-of-the-art ensemble methods across a wider range of datasets and problem types.

**Section 3.1**: The paper only addresses regression problems, using synthetic and real data from UCI. Why not include classification tasks like CIFAR-10, CIFAR-100, or ImageNet? Does the proposed method work for these? If it doesn’t, could you explain why that wasn’t mentioned?

The models used appear to be quite small, with only a few-layer MLPs employed. Exploring more advanced architectures, such as ResNet or WideResNet, could enhance the practical applicability of the findings.

**Real data is too old**:  UCI Machine Learning Repository was launched in 1987 and is too old a benchmark. How about considering latest more challenging regression benchmarks e.g. OpenML Benchmark Datasets (2013), PMLB (Penn Machine Learning Benchmark) (2018) etc.?

**Section 3.2**: The choice using SSA isn’t clearly justified. What is the reason behind this choice? How about doing an abaltion on other methods e.g.  Distributionally Robust Optimization (DRO), RO etc.?

**Computational cost not discussed**: The paper doesn’t clearly address the computational cost of the proposed method, which could be a major concern for large-scale applications.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: It would be helpful if the authors could address my questions regarding the practicality, implementation, and evaluation of the proposed method, as detailed above.
soundness: 2
strengths: The authors propose to study an important problem in ensemble learning.

The idea of using majority voting to select the best model in ensemble learning seems interesting.
summary: This paper presents a method for ensemble learning. After training B models on subsampled data from the training set, the authors propose two strategies (MoVE and ROVE / ROVEs) to select the best model from the B ensemble members. The authors proved that their method can achieve an exponentially decaying rate for excess risk. Experiments results on synthetic data and six datasets from UCI repository validate the proposed method’s performance in out-of-sample testing.
weaknesses: =============

After rebuttal: 

Thank you to the authors for their further response. I appreciate the linear regression example provided in the discussion with Reviewer aq9J. The derivation on the polynomial tail in the excess risk of the base learner and the exponential tail of the proposed method helps me better understand the technical contribution of the work and addresses my previous concerns about the base learner. Additionally, the discussion on mode estimation in models sampling is insightful. I would recommend including this discussion in the final paper, as it strengthens the rationale behind the proposed ensemble method. Taking all this into account, I believe the theoretical contributions in this work outweigh limitations in experiments. I am happy to raise my score.

=============

I have the following concerns/questions regarding the work:

1. Practicality of Theorem Assumptions: in Theorem 1, Eq. (9) is under the assumption that $n_{k, \delta} > 4/5$, which indicate that $p_k^{\text{max}} > 4/5$. This assumption appears quite strong, as it requires the most frequently subsampled model to have a subsample rate above 80%. Could the authors clarify the rationale behind this requirement?

2. In contrast to MoVE which seems to be based on a strong assumption, ROVE/ROVEs presents a more practical strategy for selecting the best model. However, the practical insights of Theorem 2 regarding the excess risk of ROVE/ROVEs remain somewhat unclear. Could the authors provide more detailed practical insights into Theorem 2?

Regarding the implementation details and evaluation, I have the following additional questions:

3. In line 193, the authors mention that "Therefore, $n_{k, \delta}$ taking large values signals the situation where the base learner already generalizes well." However, this statement depends not only on the base learner but also on the choice of $\delta$, i.e., a larger $\delta$ leads to a smaller $\mathcal{E}\_{k,\delta}$, thus a larger $n_{k, \delta}$. Could the authors provide further clarification on this claim?

4. In line 271, the authors discussed the choice of $\epsilon$ as "In our experiments, we find it a good strategy to choose an ε that leads to a maximum likelihood around 1/2." Could the authors elaborate further on this?

5. Details of the Base Learning Algorithm $\mathcal{A}$: The paper assumes that the ensemble members are trained using a generic learning algorithm $\mathcal{A}$ based on the subsampled data. Could the authors clarify details about $\mathcal{A}$? For instance, would the settings in $\mathcal{A}$ (such as learning rate, batch size in gradient-based methods) impact the performance of the proposed framework?

6. Experimental Validation and Comparisons: The current experimental results may be insufficient to support the claims fully. Firstly, one of the main contributions is the proof of an exponential decay rate of excess risk; could the authors present numerical results to validate this aspect? Secondly, the current experiments focus only on comparisons between the proposed method (MoVE and ROVE / ROVEs) and the base model. Including additional ensemble methods, as discussed in Section 4, would strengthen the experimental comparisons. Finally, the appendix mentions that the ensemble members use MLPs with 2-8 hidden layers; it would be helpful to explore whether the framework is applicable to other architectures, such as CNNs or transformers.

---

confidence: 4
contribution: 3
initial_rating: 3
presentation: 3
questions: 1) Can you give a simple self-contained example of a "reasonable" base learner that provably does not have exponential tails but has polynomial ones?  It seems to me that this means that you either have to consider infinite VC cases, or losses with unbounded moments, or algorithms that are far inferior to ERM in the limit.  There might be a way to make a convincing case for the last option, that is, inferior to ERM in the limit, but nevertheless practically reasonable in practical settings.

2) A brief review of the proofs suggests that they might repurposed to argue that ensembling can considerably improve the tails long before the asymptotic regime.  Is this correct?
soundness: 2
strengths: The paper is clearly written.  I learned something from this paper (but maybe not what the authors claim.)
summary: The paper presents an asymptotic theory of ensembling that focuses on the tail distribution of the excess error, that is, the probability that the generalization error exceeds the best possible one by more than a specified constant.  The paper shows that the proposed ensembles can have exponential tails even if the base algorithm has polynomial tails.  Experiments show that ensembles do not only improve the variance, but can have a dramatic effect on the tails.



/* Raised score after discussion */
weaknesses: I believe that the asymptotic viewpoint undermines the theoretical result. Even without going to the refined empirical process literature, we know from Vapnik (Vapnik 1998, Wiley) that if the function L(theta,z) has a finite capacity (VC dim for classification, extensions of VC dim for continuous losses with either bounded range or bounded p-th moments), we have relative uniform convergence results with asymptotically exponential tails.  Hence in that setup, any reasonable base learner will also have exponential tails.  On the other hand, if the VC dim is infinite, we know that it will be very difficult to get any guarantee at all (see also JMLR 11 (2010) 2635-2670). 

I can see therefore two ways to interpret the results of the paper:

* We can consider that the author prove that ensembling alone is as good as empirical risk optimization and variants in the sense that it is going to give exponential tails using possibly bad base algorithms.  However this amounts to saying that averaging is optimization, and we know that averaging can be superior to optimization. 

* We can also recall that the VC exponential bounds are often meaningless for practical data sizes. Okay they're exponential in the very large sample limit, but for anything practical, they don't look exponential.  Ensembling can be a way to improve dramatically these practical tails (the experiments support this point of view). But that also means that the asymptotic perspective of the theory proposed in the paper is not the right way to look at this.

Meta-review: The paper introduces MoVE (Majority Vote Ensembling) and ROVE (Retrieval and ϵ-Optimality Vote Ensembling), two ensemble learning techniques designed to achieve exponentially decaying tails for excess risk, improving generalization even when base learners exhibit slower polynomial decay rates. By selecting the best model trained on subsampled datasets and combining them using majority voting or likelihood-based selection, the methods enhance performance on out-of-sample tasks, particularly for heavy-tailed data distributions. The paper proves exponential improvement in tail bounds for excess risk, surpassing traditional variance reduction techniques. The authors validate the methods on synthetic and real-world datasets, demonstrating improved generalization over base learners. The paper proposes an innovative approach to ensemble learning with exponential tail improvements, with a rigorous theoretical foundation complemented by experiments.  However, the paper also has some weaknesses, as indicated by the reviewers. It uses outdated datasets, lacks comparisons with state-of-the-art methods, and explores only small model architectures.  Regarding the theoretical aspect, the results are abstract and hard to interpret without concrete examples.   The reviewers give mixed scores (ranging from reject to accept) highlighting the novelty of the theoretical contribution and the shortcomings in experimental validation and practical insights.  Recommendations include testing on more diverse datasets, improving clarity, and explaining computational costs in the paper. Summing up, while the paper introduces promising theoretical advancements in ensemble learning, I believe that it requires significant revisions in experiments, presentation, and practical applications to fully demonstrate its impact.

Decision: Reject