# Bagging Improves Generalization Exponentially

Huajie Qian
DAMO Academy
Alibaba Group
Bellevue, WA 98004
<EMAIL_ADDRESS>
&Donghao Ying∗
IEOR Department
UC Berkeley
Berkeley, CA 94720
<EMAIL_ADDRESS>
&Henry Lam
IEOR Department
Columbia University
New York, NY 10027
<EMAIL_ADDRESS>
&Wotao Yin
DAMO Academy
Alibaba Group
Bellevue, WA 98004
<EMAIL_ADDRESS>
Equal contribution.

###### Abstract

Bagging is a popular ensemble technique to improve the accuracy of machine
learning models. It hinges on the well-established rationale that, by
repeatedly retraining on resampled data, the aggregated model exhibits lower
variance and hence higher stability, especially for discontinuous base
learners. In this paper, we provide a new perspective on bagging: By suitably
aggregating the base learners at the parametrization instead of the output
level, bagging improves generalization performances exponentially, a strength
that is significantly more powerful than variance reduction. More precisely,
we show that for general stochastic optimization problems that suffer from
slowly (i.e., polynomially) decaying generalization errors, bagging can
effectively reduce these errors to an exponential decay. Moreover, this power
of bagging is agnostic to the solution schemes, including common empirical
risk minimization, distributionally robust optimization, and various
regularizations. We demonstrate how bagging can substantially improve
generalization performances in a range of examples involving heavy-tailed data
that suffer from intrinsically slow rates.

##  1 Introduction

Bootstrap aggregating, or bagging, is a popular ensemble technique to improve
the accuracy of machine learning models [1]. It comprises repeated resampling
of the data to retrain models (the “base learners”), which are then aggregated
through averaging or majority vote. In the literature, the main justification
for bagging pertains to variance reduction or higher stability thanks to its
smoothing effect. This justification has been shown to be particularly
relevant for certain U-statistics [2], and models with hard-thresholding rules
[3] such as linear regression with variable selection and decision trees that
give rise to random forests [4].

Contrary to the established understanding, in this paper we present a new view
of bagging in offering an arguably stronger power than variance reduction: By
suitably aggregating the base learners at the parametrization instead of the
output level, bagging can provide an _exponential_ improvement in
generalization. More precisely, we show that, for general stochastic
optimization problems that suffer from a slow, namely polynomial, decay in
generalization errors, bagging can reduce these errors to an exponential
decay. Thus, instead of the typical constant factor of improvement exhibited
by variance reduction, bagging 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 bagging perspective, proposed procedures, and the
technicality involved in a full analysis. Extensive experiments show that our
improvement is particularly significant for discrete-decision problems, which
is in a sense consistent with past literature on tree-based models where
conventional bagging has routinely been applied.

#### Main results at a high level.

We begin by introducing a generic stochastic optimization problem

$\min_{x\in\mathcal{X}}Z(x):=\mathbb{E}\left[h(x,\xi)\right],$ (1)

where $x$ is the decision variable on space $\mathcal{X}$, $\xi\in\Xi$ denotes
the randomness governed by a probability distribution $F$, and $h$ is the cost
function. Typically, $F$ is only observed from noisy data, say i.i.d. samples
$\\{\xi_{1},\ldots,\xi_{n}\\}$. In machine learning, $x$ corresponds to model
parameters (e.g., regression coefficients), $\\{\xi_{1},\ldots,\xi_{n}\\}$ the
training data, $h$ the loss function, and $Z$ the population-level expected
loss. More generally, (1) encapsulates data-driven decision-making problems,
namely the integration of data on $\xi$ into a downstream optimization task
with overall cost function $h$ and prescriptive decision $x$. These problems
are increasingly prevalent in various industrial applications [5, 6, 7], such
as in supply chain network design where $x$ may represent the decision to open
processing facilities, $\xi$ the variability in supply and demand, and $h$ the
total cost of processing and transportation.

Given the data, we can train the model or decision by approaches such as
empirical risk minimization or sample average approximation (SAA) [8, 9],
distributionally robust optimization (DRO) [10, 11], and various
regularizations. Our proposal and theory described below are agnostic to which
empirical approach to use, as long as it satisfies the presented
prerequisites.

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

$\mathbb{P}\Big{(}Z(\hat{x})>\min_{x\in\mathcal{X}}Z(x)+\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 [12, 13, 14] due to slow concentration,
which frequently arise in machine learning applications such as large language
models (e.g, [15, 16, 17] among others), finance [18, 19] and physics [20,
21], and are proved to be tight [22] for empirical risk minimization. As our
key insight, our proposed bagging methodology can improve the above to an
exponential decay, i.e.,

$\mathbb{P}\Big{(}Z(\hat{x})>\min_{x\in\mathcal{X}}Z(x)+\delta\Big{)}\leq
C_{2}\gamma^{n/k},$ (3)

where $k$ is the resample data size used in the bagging 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 is truly
exponential. This exponential acceleration is qualitatively different from the
well-known variance reduction benefit of bagging in several aspects. First,
the variance reduction refers to the smaller variability in predictions from
models trained on independent data sets and thus has a more direct impact on
the expected regret than the tail decay rate. Second, the improvement by
variance reduction is typically of a constant factor (e.g., [3] 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
$\mathcal{X}$. Our bagging methodology uses a majority-vote mechanism at the
model level: After we repeatedly resample data to retrain models, we output
the model that appears most frequently among $\mathcal{X}$. This implicitly
solves a “meta-optimization” problem that shares the same decision space
$\mathcal{X}$ with (1) but maximizes the probability of being selected as the
trained model from the data. This conversion of the original training problem
from the 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 resample and original data,
this in turn results in exponentially decaying generalization errors for the
overall regret.

For more general problems with continuous solution space, we can replace the
simple majority vote presented above with a vote on the highest likelihood of
falling into an $\epsilon$-optimal set. This avoids the degeneracy issue in
using simple majority vote for continuous problems, while retains 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 empirical performance is especially significant in this latter case.

The rest of the paper is organized as follows. Section 2 presents our bagging
methodologies and their finite-sample bounds. Section 3 presents promising
experimental results. Technical proofs, review of related work, and additional
technical and experimental results can be found in Appendix.

##  2 Methodology and theoretical guarantees

To solve (1) using data, we consider the empirical problem

$\min_{x\in\mathcal{X}}\ \hat{Z}_{n}(x;\xi_{1},\ldots,\xi_{n}),$ (4)

where $\hat{Z}_{n}(x;\xi_{1},\ldots,\xi_{n})$ is an empirical estimate of
$Z(x)$ based on i.i.d. observations $(\xi_{1},\ldots,\xi_{n})$. This
formulation includes SAA, DRO and various regularizations, and gives rise to
our base learner. For convenience, we shorthand
$\hat{Z}_{n}(x;\xi_{1},\ldots,\xi_{n})$ as $\hat{Z}_{n}(x)$ when no confusion
arises.

###  2.1 A basic procedure

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

Algorithm 1 Bagging Models via Majority Vote (BAG)

1: Input: $n$ i.i.d. observations $\bm{\xi}_{1:n}=(\xi_{1},\ldots,\xi_{n})$,
subsample size $k<n$ and Monte Carlo size $B$.

2: for $b=1$ to $B$ do

3: Randomly sample $\bm{\xi}_{k}^{b}=(\xi_{1}^{b},\ldots,\xi_{k}^{b})$
uniformly from $\bm{\xi}_{1:n}$ without replacement, and obtain
$\hat{x}_{k}^{b}\in\operatorname{argmin}_{x\in\mathcal{X}}\hat{Z}_{k}(x;\xi_{1},\ldots,\xi_{k})$.

4: end for

5: Output:
$\hat{x}^{BAG}_{n,k}\in\operatorname{argmax}_{x\in\mathcal{X}}\sum_{b=1}^{B}\mathbbm{1}(\hat{x}^{b}_{k}=x)$.

To understand Algorithm 1, denote by $\hat{x}_{k}(\xi_{1},\ldots,\xi_{k})$ the
solution to (4) with input data $(\xi_{1},\ldots,\xi_{k})$. Associated with
(1), we define a meta-optimization problem

$\max_{x\in\mathcal{X}}\
p_{k}(x):=\mathbb{P}\left(x=\hat{x}_{k}(\xi_{1},\ldots,\xi_{k})\right),$ (5)

which maximizes the probability of a model being output by the weak learner
trained using $k$ observations. If $B=\infty$, Algorithm 1 essentially
maximizes an empirical approximation of (5), i.e.

$\max_{x\in\mathcal{X}}\
\mathbb{P}_{*}\left(x=\hat{x}_{k}(\xi^{*}_{1},\ldots,\xi^{*}_{k})\right),$ (6)

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

$\max_{x\in\mathcal{X}}\
\frac{1}{B}\sum_{b=1}^{B}\mathbbm{1}(x=\hat{x}_{k}(\xi_{1}^{b},\ldots,\xi_{k}^{b})),$
(7)

which gives exactly the output of Algorithm 1. In other words, Algorithm 1 is
a _bootstrap approximation_ to the solution of (5). The following result
materializes the intuition explained in the introduction on the conversion of
the original potentially heavy-tailed problem (4) into a probability
maximization (7) that possess exponential bounds:

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

Consider discrete decision space $\mathcal{X}$. Recall $p_{k}(x)$ defined in
(5). Let $p_{k}^{\max}:=\max_{x\in\mathcal{X}}p_{k}(x)$,
$\mathcal{X}^{\delta}:=\left\\{x\in\mathcal{X}:Z(x)\leq\min_{x^{\prime}\in\mathcal{X}}Z(x^{\prime})+\delta\right\\}$
be the set of $\delta$-optimal solutions of problem (1),
$D_{\operatorname{KL}}(p\|q):=p\ln\frac{p}{q}+(1-p)\ln\frac{1-p}{1-q}$ be the
Kullback–Leibler divergence between two Bernoulli distributions with means $p$
and $q$, and

$\eta_{k,\delta}:=\max_{x\in\mathcal{X}^{\delta}}p_{k}(x)-\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x),$
(8)

where $\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)$ evaluates
to $0$ if $\mathcal{X}\backslash\mathcal{X}^{\delta}$ is empty. For every
$k\leq n$ and $\delta\geq 0$ such that $\eta_{k,\delta}>0$, the solution
output by Algorithm 1 satisfies that

$\displaystyle\
\mathbb{P}\left(Z(\hat{x}_{n,k}^{BAG})>\min_{x\in\mathcal{X}}Z(x)+\delta\right)$
(9) $\displaystyle\leq$
$\displaystyle\left|\mathcal{X}\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$, (9) is further bounded by

$\left|\mathcal{X}\right|\left(3\min\left(e^{-2/5},C_{1}\max(1-p_{k}^{\max},\;p_{k}^{\max}-\eta_{k,\delta})\right)^{\frac{n}{C_{2}k}}+\exp\left(-\frac{B}{C_{3}}\right)\right),$
(10)

where $C_{1},C_{2},C_{3}>0$ are universal constants.

Theorem 1 states that the generalization bound using Algorithm 7 decays
exponentially in the ratio $n/k$ and the number of base learners $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 (5) with (6). The second part quantifies the Monte Carlo error in
approximating (6) with finite $B$. The third part comes from the interaction
between the two sources of errors and is typically of higher order. The
multiplier $\left|\mathcal{X}\right|$ in the bound is avoidable, e.g., via a
space reduction as in our next algorithm.

The condition $\eta_{k,\delta}>0$ on the bound plays two roles. First, it
measures how much the optimality of the original problem (1) can be propagated
to that of the meta problem (5), in the sense that every solution that is
$\eta_{k,\delta}$-optimal for (5) is $\delta$-optimal for (1). Second, the
term $\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)$ itself
resembles the generalization bound of the base learner and
$\max_{x\in\mathcal{X}^{\delta}}p_{k}(x)$ captures the concentration on
$\delta$-optimal solutions, so that $\eta_{k,\delta}$ implicitly encodes the
generalization performance of the base learner. Indeed, when $p_{k}^{\max}$
and $\eta_{k,\delta}$ both take large values, this signals the situation where
the base learner already generalizes well. In this case, (9) can be simplified
to (10). The bound (10) suggests that our bagging does not hurt the
performance of an already high-performing base learner as the performance of
the base learner is captured by the
$\max(1-p_{k}^{\max},\;p_{k}^{\max}-\eta_{k,\delta})$ term in the bound. See
Appendix B for a more detailed comparison.

The quantity $\eta_{k,\delta}$ also hints on how to choose subsample size $k$.
As long as $\eta_{k,\delta}$ is bounded away from $0$, our bound decays
exponentially fast. Therefore, $k$ can be chosen in such a way that the base
learner starts to show some generalization power 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 (9) 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
a $B=\mathcal{O}(n/k)$ is sufficient to control the Monte Carlo error to a
similar magnitude as the data error.

With more careful treatment, the generalization sensitivity $\eta_{k,\delta}$
and the maximum probability $p_{k}^{\max}$ can be more explicitly bounded in
terms of the tail probability for the maximum deviation of the empirical
estimate $\hat{Z}_{n}$, defined as

$T_{n}(\epsilon):=\mathbb{P}\left(\sup_{x\in\mathcal{X}}\left|\hat{Z}_{n}(x)-Z(x)\right|\geq\epsilon\right),\text{\
\ for }\epsilon\geq 0,$ (11)

which connects the performance of Algorithm 1 to solution multiplicity of (1).
This is deferred to Appendix B.

###  2.2 A more general procedure

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

Algorithm 2 Bagging Models via $\epsilon$-Optimality Vote (ReBAG / ReBAG-S)

Input: $n$ i.i.d. observations $\bm{\xi}_{1:n}=(\xi_{1},\ldots,\xi_{n})$,
subsample size $k<n/2$, Monte Carlo sizes $B_{1}$ and $B_{2}$.

Phase I: Model Candidate Retrieval

for $b=1$ to $B_{1}$ do

Randomly sample $\bm{\xi}_{k}^{b}=(\xi_{1}^{b},\ldots,\xi_{k}^{b})$ uniformly
from $\bm{\xi}_{1:n}$ (if no split) or
$\bm{\xi}_{1:\lfloor\frac{n}{2}\rfloor}$ (if split) without replacement, and
obtain
$\hat{x}^{b}_{k}\in\operatorname{argmin}_{x\in\mathcal{X}}\hat{Z}_{k}(x;\xi_{1}^{b},\ldots,\xi_{k}^{b})$.

end for

Let $\mathcal{S}:=\\{\hat{x}^{b}_{k}:b=1,\ldots,B_{1}\\}$ be the set of all
retrieved models.

Phase II: $\epsilon$-Optimality Vote

Choose $\epsilon\geq 0$ using the data $\bm{\xi}_{1:n}$ (if no split) or
$\bm{\xi}_{1:\lfloor\frac{n}{2}\rfloor}$ (if split).

for $b=1$ to $B_{2}$ do

Randomly sample $\bm{\xi}_{k}^{b}=(\xi_{1}^{b},\ldots,\xi_{k}^{b})$ uniformly
from $\bm{\xi}_{1:n}$ (if no split) or
$\bm{\xi}_{\lfloor\frac{n}{2}\rfloor+1:n}$ (if split) without replacement, and
calculate
$\widehat{\mathcal{X}}^{\epsilon,b}_{k}:=\left\\{x\in\mathcal{S}:\hat{Z}_{k}(x;\xi_{1}^{b},\ldots,\xi_{k}^{b})\leq\min_{x^{\prime}\in\mathcal{S}}\hat{Z}_{k}(x^{\prime};\xi_{1}^{b},\ldots,\xi_{k}^{b})+\epsilon\right\\}.$

end for

Output:
$\hat{x}^{BAG}_{n,k}\in\operatorname{argmax}_{x\in\mathcal{S}}\sum_{b=1}^{B_{2}}\mathbbm{1}(x\in\widehat{\mathcal{X}}^{\epsilon,b}_{k})$.

This procedure, displayed in Algorithm 2, proceeds initially the same as
Algorithm 1 in repeatedly resampling data and retraining the model using (4).
However, in the aggregation step, instead of using simple majority vote,
Algorithm 2 outputs, among all the retrained models, the one that has the
highest likelihood of being $\epsilon$-optimal. This $\epsilon$-optimality
avoids the degeneracy of majority vote and, moreover, since we have restricted
our output to the collection of retrained models, the corresponding likelihood
maximization is readily doable by simple enumeration. We have the following
theoretical guarantees for Algorithm 2:

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

Consider Algorithm 2 with data splitting, i.e., ReBAG-S. Recall the tail
probability $T_{n}$ in (11). 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}\left((\delta-\overline{\epsilon})/2\right)+T_{k}\left(\underline{\epsilon}/2\right)<1/5$,
then

$\displaystyle\mathbb{P}\Big{(}Z(\hat{x}_{n,k}^{BAG})>\min_{x\in\mathcal{X}}Z(x)+2\delta\Big{)}\leq$
$\displaystyle
B_{1}\left[3\min\left(e^{-2/5},C_{1}T_{k}\left(\dfrac{\min(\underline{\epsilon},\delta-\overline{\epsilon})}{2}\right)\right)^{\frac{n}{2C_{2}k}}+e^{-{B_{2}}/{C_{3}}}\right]$
(12) $\displaystyle+$
$\displaystyle\min\left(e^{\frac{-(1-T_{k}(\delta/2))}{C_{4}}},C_{5}T_{k}\left({\delta}/{2}\right)\right)^{\frac{n}{2C_{6}k}}+e^{-\frac{B_{1}(1-T_{k}(\delta/2))}{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., ReBAG, and discrete
decision space $\mathcal{X}$. Assume $\lim_{n\to\infty}T_{n}(\epsilon)=0$ for
all $\epsilon>0$. Then, for every fixed $\delta>0$, we have
$\mathbb{P}\left(Z(\hat{x}_{n,k}^{BAG})>\min_{x\in\mathcal{X}}Z(x)+2\delta\right)\to
0$ as $n\to\infty$, if $\mathbb{P}\left(\epsilon>\delta/2\right)\to 0$,
$k\to\infty$, $n/k\to\infty$ and $B_{1},B_{2}\to\infty$ as $n\to\infty$.

Theorem 2 provides an exponential generalization bound, regardless of discrete
or continuous space. The first line in the bound (12) is inherited from the
bound (10) for Algorithm 1 from majority to $\epsilon$-optimality vote. In
particular, the multiplier $\left|\mathcal{X}\right|$ in (10) is now replaced
by $B_{1}$, the number of retrieved models. The second line in (12) bounds the
performance sacrifice due to the restriction of Phase I model candidates.
Algorithm 2 is found to be more robust than Algorithm 1, as its bound is not
affected the possible presence of multiple (nearly) optimal solutions for
which Algorithm 1 may not concentrate as well due to the spread of the votes
to each individual optimum. More details can be found in Appendix B.

Algorithm 2 may be carried out with (ReBAG-S) or without (ReBAG) splitting the
data between the two phases. Data splitting makes the procedure theoretically
more tractable by avoiding inter-dependency between the phases, but sacrifices
some statistical power by effectively halving the data size. Empirically we
find ReBAG 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 established
through several novel techniques. The first is a sharper concentration result
for U-statistics with binary kernels than well-known Bernstein-type
inequalities (e.g., [23, 24]). This is needed to recover the correct order of
the bound, in particular the that captures the convergence of not only the
bootstrap approximation in bagging, but also the base learner, especially
those already possessing exponential decay rates. Thus it provides sharp
insights for the robustness of bagging for fast convergent base learners. The
second is a new sensitivity analysis on the regret for the original problem
with respect to that of the accompanying probability maximization, which is
needed to translate the superior generalization performance in the
accompanying problem into the accelerated convergence for the original
problem. Lastly, in establishing asymptotic consistency for Algorithm 2
without data splitting, in order to allow a randomly chosen $\epsilon$
tolerance, we develop a uniform law of large numbers (LLNs) for the class of
events of being included in the $\epsilon$-optimal set via a direct analysis
of the second moment of the maximum deviation. Uniform LLNs are particularly
challenging for this class because, unlike the standard settings where the
class is fixed, this particular class itself depends on the subsample size $k$
as $n\to\infty$ and therefore is of a dynamic nature.

##  3 Numerical experiments

In this section, we numerically validate Algorithm 1 (BAG), Algorithm 2
(ReBAG), and Algorithm 2 with data splitting (ReBAG-S). Specifically, we
consider the following problems: resource allocation, supply chain network
design, portfolio optimization, model selection, maximum weight matching, and
linear programming. To ensure uniformity in presentation, we convert all
problems into minimization forms. All experiments are repeated more than 50
times with the average performance being reported, and they are carried out on
a personal computer with Apple M1 Pro chip and 16 GB memory. Due to space
limitations, we defer the detailed introduction about the problems and some
experimental results to Appendix D.

###  3.1 Hyperparameter profiling

We first test the effect of different hyperparameters in the bagging
framework, including subsample size $k$, number of subsamples $B,B1,B2$, and
threshold $\epsilon$. Throughout this profiling stage, we use the sample
average approximation (SAA) as the base model. To profile the effect of
subsample size and numbers, we consider the resource allocation problem.

(a) Profiling for $k$ (BAG).

(b) Profiling for $\epsilon$ (multiple optima).

(c) Profiling for $\epsilon$ (unique optima).

Figure 1: Profiling for subsample size $k$ and threshold $\epsilon$. (a):
Resource allocation problem, where $B=200$; (b) and (c): Linear program, where
$k=\max(10,0.005n)$, $B_{1}=20$, and $B_{2}=200$.

Subsample size We explored scenarios where $k$ is both dependent on and
independent of the total sample size $n$ (see Figures 1(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 1(a), $k=2$ proves inadequate, with larger $k$
delivering good results. The underlying reason is that the effective
performance of BAG requires
$x^{*}\in\operatorname{argmax}_{x\in\mathcal{X}}p_{k}(x)$. In the former, this
is achieved with only two samples, enabling BAG to identify $x^{*}$ 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 9, we simulate $p_{k}(x)$ for the two cases, which further explains the
influence of the subsample size.

Number of subsamples In Figure 7, we illustrate the performance of BAG and
ReBAG under different $B,B_{1},B2$, where we set $k=10$ for both methods and
$\epsilon=0.005$ for ReBAG. From the figure, we find that the performance of
bagging is increasing in the number of subsample numbers.

Threshold $\epsilon$ The optimal choice of $\epsilon$ in ReBAG and ReBAG-S is
problem-dependent and related to the number of (near) optimal solutions. This
dependence is illustrated by the performance of ReBAG shown in Figures 1(b)
and 1(c). Hence, we propose an adaptive strategy defined as follows: Let
$g(\epsilon):=1/B_{2}\cdot\sum_{b=1}^{B_{1}}\mathbbm{1}(\hat{x}^{BAG}_{n,k}(\epsilon)\in\widehat{\mathcal{X}}^{\epsilon,b}_{k})$,
where we use $\hat{x}^{BAG}_{n,k}(\epsilon)$ to emphasize the dependency of
$\hat{x}^{BAG}_{n,k}$ 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{x}^{BAG}_{n,k}(\epsilon)$ is included in the “near optimum set”
$\widehat{\mathcal{X}}^{\epsilon,b}_{k}$. The choice of $\epsilon^{*}$ makes
it more likely for the true optimal solution to be included in the evaluation
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 $\bm{\xi}_{1:\lfloor\frac{n}{2}\rfloor}$ (Phase I data)
for ReBAG-S. From Figure 1, we observe that this adaptive strategy exhibits
decent performance for all scenarios. Similar behaviors can also be observed
for ReBAG-S in Figure 8.

###  3.2 Algorithm comparison

In this section, we demonstrate how the bagging framework enhances the
performance of base models. Based on the profiling results, we choose
$k=\max(10,0.005n)$, $B=200$, $B_{1}=20$, and $B_{2}=200$ as default
hyperparameters and adopt the proposed adaptive strategy to determine
$\epsilon$. Note that the choice of $k$ ensures efficient use of samples,
particularly when the total sample size $n$ is large.

(a) Resource allocation.

(b) Network design.

(c) Portfolio optimization.

(d) Model selection.

(e) Maximum weight matching.

(f) Linear program (multiple optima).

Figure 2: Performance of bagging approaches using SAA as the base model.
Unmarked problems have a single optimum. Hyperparameters: $k=\max(10,0.005n)$,
$B=200$, $B_{1}=20$, and $B_{2}=200$.

We first use SAA as the base model. Figure 2 displays the comparison results
in the six problems considered in the paper, where we observe that bagging
approaches consistently outperform the base model. Particularly, when there is
a single optimal solution, we find that BAG performs slightly better than
ReBAG and ReBAG-S in general (see Figures 2(a)-2(e)). Conversely, with
multiple optima, ReBAG and ReBAG-S demonstrate dominant performances. The
underlying reason is that the presence of multiple optima would decrease the
probability for the base model to output any specific optimal solution,
thereby making optimal and suboptimal solutions less distinguishable, which
results in the slightly worse performance of BAG. On the other side, ReBAG and
ReBAG-S alleviate this issue by including more candidate solutions into a
solution ranking phase. In Figure 10, we plot the generalization sensitivity
$\eta_{k,\delta}$ for different bagging approaches, further explaining the
superiority of ReBAG and ReBAG-S when multiple optimal occur. In Figure 11, we
also compare the bagging approaches (using SAA as the base model) with
Wasserstein metric-based distributionally robust optimization (DRO). We
observe that though DRO performs slightly better than SAA under some choices
of the ambiguity set radius, the improvement is marginal compared to the
bagging approaches.

Finally, in Figure 3, we apply the bagging framework using DRO with
Wasserstein metric as the base model. The experimental result demonstrates
that our bagging framework is not bound to some specific base model but
applies to all suitable ones.

(a) Resource allocation (Figure 2(a)).

(b) Maximum weight matching (Figure 2(e)).

Figure 3: Performance of bagging approaches using DRO with Wasserstein metric
as the base model. Hyperparameters: $k=\max(10,0.005n)$, $B=200$, $B_{1}=20$,
and $B_{2}=200$.

###  3.3 Further discussion

Based on the experimental results, we provide further discussions and offer
practical recommendations about the proposed bagging approaches. Firstly,
based on the comparison results in Figures 2 and 3, we recommend using BAG
when it is known a prior that the problem admits a unique solution. Otherwise,
using ReBAG and ReBAG-S is a better choice as these algorithms deliver stable
results across all scenarios and perform exceptionally well when there exist
multiple optima.

(a) Instance 1.

(b) Instance 2.

(c) Instance 3.

Figure 4: Performance of bagging approaches using SAA as the base model in
linear programs with light-tailed objectives. Hyperparameters:
$k=\max(10,0.005n)$, $B=200$, $B_{1}=20$, and $B_{2}=200$.

(a) Influence of tail heaviness.

(b) Running time comparison.

Figure 5: (a): Influence of tail heaviness when using SAA as the base model in
a linear program with multiple optima. Hyperparameters: $n=1000000$, $k=50$,
$B=2000$, $B_{1}=200$, $B_{2}=5000$. The tail heaviness parameter corresponds
to the mean of the Pareto random variable. (b): Running time of bagging
approaches using SAA as the base model in a network design problem.
Hyperparameters: $k=10$, $B=200$, $B_{1}=20$, $B_{2}=200$. “Sequential” refers
to using a simple “for” loop for the bagging process; “Parallel” refers to
using 8 processors to solve the bagging process in parallel.

Secondly, we note that our bagging framework is not limited to heavy-tailed
problems only. In light-tailed scenarios where the base model may already
enjoy fast convergence, our bagging approaches can also demonstrate comparable
or superior performance (see Figure 4). Moreover, the advantage of our methods
becomes more prominent as the underlying random variable goes more heavy-
tailed. In Figure 5(a), we observe that while the base model SAA is sensitive
to the tail heaviness, bagging approaches, especially ReBAG and ReBAG-S, are
generally much more stable (note that the relatively weak performance of BAG
is due to the precense of multiple optima). Hence, when the noises/randomness
are heavy-tailed, it is recommended to use bagging approaches to improve and
stabilize the performance of base models.

Finally, we note that our bagging framework, which requires solving multiple
sample-based optimization problems, does not impose a significantly higher
computational burden compared to plain data-driven methods that solve only one
(see Figure 5(b)), and it can even be computationally advantageous for certain
problems, such as mixed integer programs (especially when the sample size is
large). For instance, methods like DRO [25, 10] and two-stage stochastic
optimization often expand the decision space/constraints, leading to an actual
problem size that grows at least linearly with the sample size. Conversely,
the subproblems in the bagging process, due to the subsampling, are usually
smaller and quicker to solve. Moreover, our theories reveal that the number of
subproblems need not exceed $\mathcal{O}(n/k)$, as solving more does not
further improve generalization performance. Since bagging processes can be
easily parallelized, the total computational effort remains comparable to
distributing the full data set across $n/k$ machines, each handling a
subsample of size $k$.

##  4 Conclusion and limitation

In this paper, we present a novel perspective on bagging, demonstrating its
ability to exponentially enhance generalization performance by aggregating
base learners at the parameter level. Our framework, which operates
independently of the solution method, is particularly effective for stochastic
optimization problems with polynomially decaying generalization errors. We
validate this approach through extensive numerical examples involving heavy-
tailed data, highlighting bagging’s substantial benefits in scenarios with
intrinsically slow convergence rates. This work underscores the powerful
potential of bagging across a broad range of machine learning applications.

For the limitation, we remark that although our bagging framework accelerates
tail convergence of the regret, it does not necessarily reduce the expected
regret except for discrete problems as discussed before. In addition, like
existing bagging procedures, our framework may increase the model bias due to
the use of subsampling [2], therefore it’s best suited for problems with
relatively low biases.

## References

  * [1] Leo Breiman.  Bagging predictors.  Machine learning, 24:123–140, 1996.
  * [2] Andreas Buja and Werner Stuetzle.  Observations on bagging.  Statistica Sinica, pages 323–351, 2006.
  * [3] Peter Bühlmann and Bin Yu.  Analyzing bagging.  The annals of Statistics, 30(4):927–961, 2002.
  * [4] Leo Breiman.  Random forests.  Machine learning, 45:5–32, 2001.
  * [5] Sachin S Kamble, Angappa Gunasekaran, and Shradha A Gawankar.  Achieving sustainable performance in a data-driven agriculture supply chain: A review for research and applications.  International Journal of Production Economics, 219:179–194, 2020.
  * [6] Dimitris Bertsimas, Shimrit Shtern, and Bradley Sturt.  A data-driven approach to multistage stochastic linear optimization.  Management Science, 69(1):51–74, 2023.
  * [7] Shubhechyya Ghosal, Chin Pang Ho, and Wolfram Wiesemann.  A unifying framework for the capacitated vehicle routing problem under risk and ambiguity.  Operations Research, 72(2):425–443, 2024.
  * [8] Vladimir Vapnik.  Principles of risk minimization for learning theory.  Advances in neural information processing systems, 4, 1991.
  * [9] Alexander Shapiro, Darinka Dentcheva, and Andrzej Ruszczynski.  Lectures on stochastic programming: modeling and theory.  SIAM, 2021.
  * [10] Peyman Mohajerin Esfahani and Daniel Kuhn.  Data-driven distributionally robust optimization using the wasserstein metric: Performance guarantees and tractable reformulations.  Mathematical Programming, 171(1):115–166, 2018.
  * [11] Rui Gao and Anton Kleywegt.  Distributionally robust stochastic optimization with wasserstein distance.  Mathematics of Operations Research, 48(2):603–655, 2023.
  * [12] Vlasta Kaňková and Michal Houda.  Thin and heavy tails in stochastic programming.  Kybernetika, 51(3):433–456, 2015.
  * [13] Jie Jiang, Zhiping Chen, and Xinmin Yang.  Rates of convergence of sample average approximation under heavy tailed distributions.  To preprint on Optimization Online, 2020.
  * [14] Jie Jiang and Shengjie Li.  On complexity of multistage stochastic programs under heavy tailed distributions.  Operations Research Letters, 49(2):265–269, 2021.
  * [15] Hamid Jalalzai, Pierre Colombo, Chloé Clavel, Eric Gaussier, Giovanna Varni, Emmanuel Vignon, and Anne Sabourin.  Heavy-tailed representations, text polarity classification & data augmentation.  Advances in Neural Information Processing Systems, 33:4295–4307, 2020.
  * [16] Jingzhao Zhang, Sai Praneeth Karimireddy, Andreas Veit, Seungyeon Kim, Sashank Reddi, Sanjiv Kumar, and Suvrit Sra.  Why are adaptive methods good for attention models?  Advances in Neural Information Processing Systems, 33:15383–15393, 2020.
  * [17] Ashok Cutkosky and Harsh Mehta.  High-probability bounds for non-convex stochastic optimization with heavy tails.  Advances in Neural Information Processing Systems, 34:4883–4895, 2021.
  * [18] Georg Mainik, Georgi Mitov, and Ludger Rüschendorf.  Portfolio optimization for heavy-tailed assets: Extreme risk index vs. markowitz.  Journal of Empirical Finance, 32:115–134, 2015.
  * [19] Manfred Gilli and Evis Këllezi.  An application of extreme value theory for measuring financial risk.  Computational Economics, 27:207–228, 2006.
  * [20] Jean-Yves Fortin and Maxime Clusel.  Applications of extreme value statistics in physics.  Journal of Physics A: Mathematical and Theoretical, 48(18):183001, 2015.
  * [21] Anna PM Michel and Alan D Chave.  Analysis of laser-induced breakdown spectroscopy spectra: the case for extreme value statistics.  Spectrochimica Acta Part B: Atomic Spectroscopy, 62(12):1370–1378, 2007.
  * [22] Olivier Catoni.  Challenging the empirical mean and empirical variance: A deviation study.  Annales de l’IHP Probabilités et statistiques, 48(4):1148–1185, 2012.
  * [23] Miguel A Arcones.  A bernstein-type inequality for u-statistics and u-processes.  Statistics & probability letters, 22(3):239–247, 1995.
  * [24] Thomas Peel, Sandrine Anthoine, and Liva Ralaivola.  Empirical bernstein inequalities for u-statistics.  Advances in Neural Information Processing Systems, 23, 2010.
  * [25] Aharon Ben-Tal, Dick Den Hertog, Anja De Waegenaere, Bertrand Melenberg, and Gijs Rennen.  Robust solutions of optimization problems affected by uncertain probabilities.  Management Science, 59(2):341–357, 2013.
  * [26] Max Biggs, Rim Hariss, and Georgia Perakis.  Constrained optimization of objective functions determined from random forests.  Production and Operations Management, 32(2):397–415, 2023.
  * [27] Georgia Perakis and Leann Thayaparan.  Umotem: Upper bounding method for optimizing over tree ensemble models.  Available at SSRN 3972341, 2021.
  * [28] Keliang Wang, Leonardo Lozano, Carlos Cardonha, and David Bergman.  Optimizing over an ensemble of neural networks.  arXiv preprint arXiv:2112.07007, 2021.
  * [29] Max Biggs and Georgia Perakis.  Tightness of prescriptive tree-based mixed-integer optimization formulations.  arXiv preprint arXiv:2302.14744, 2023.
  * [30] Edward Anderson and Harrison Nguyen.  When can we improve on sample average approximation for stochastic optimization?  Operations Research Letters, 48(5):566–572, 2020.
  * [31] John R Birge.  Uses of sub-sample estimates to reduce errors in stochastic optimization models.  arXiv preprint arXiv:2310.07052, 2023.
  * [32] Henry Lam and Huajie Qian.  Assessing solution quality in stochastic optimization via bootstrap aggregating.  In Proceedings of the 2018 Winter Simulation Conference, pages 2061–2071. IEEE, 2018.
  * [33] Henry Lam and Huajie Qian.  Bounding optimality gap in stochastic optimization via bagging: Statistical efficiency and stability.  arXiv preprint arXiv:1810.02905, 2018.
  * [34] Xiaotie Chen and David L Woodruff.  Distributions and bootstrap for data-based stochastic programming.  Computational Management Science, 21(1):33, 2024.
  * [35] Xiaotie Chen and David L Woodruff.  Software for data-based stochastic programming using bootstrap estimation.  INFORMS Journal on Computing, 35(6):1218–1224, 2023.
  * [36] Andreas Eichhorn and Werner Römisch.  Stochastic integer programming: Limit theorems and confidence intervals.  Mathematics of Operations Research, 32(1):118–135, 2007.
  * [37] Miles Lopes, Shusen Wang, and Michael Mahoney.  Error estimation for randomized least-squares algorithms via the bootstrap.  In International Conference on Machine Learning, pages 3217–3226. PMLR, 2018.
  * [38] Jessie XT Chen and Miles Lopes.  Estimating the error of randomized newton methods: A bootstrap approach.  In International Conference on Machine Learning, pages 1649–1659. PMLR, 2020.
  * [39] Yixin Fang, Jinfeng Xu, and Lei Yang.  Online bootstrap confidence intervals for the stochastic gradient descent estimator.  Journal of Machine Learning Research, 19(78):1–21, 2018.
  * [40] Yanjie Zhong, Todd Kuffner, and Soumendra Lahiri.  Online bootstrap inference with nonconvex stochastic gradient descent estimator.  arXiv preprint arXiv:2306.02205, 2023.
  * [41] Sungjin Im, Benjamin Moseley, and Kirk Pruhs.  Stochastic scheduling of heavy-tailed jobs.  In 32nd International Symposium on Theoretical Aspects of Computer Science (STACS 2015). Schloss-Dagstuhl-Leibniz Zentrum für Informatik, 2015.
  * [42] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.  Language models are few-shot learners.  Advances in neural information processing systems, 33:1877–1901, 2020.
  * [43] Josh Achiam, Steven Adler, Sandhini Agarwal, Lama Ahmad, Ilge Akkaya, Florencia Leoni Aleman, Diogo Almeida, Janko Altenschmidt, Sam Altman, Shyamal Anadkat, et al.  Gpt-4 technical report.  arXiv preprint arXiv:2303.08774, 2023.
  * [44] Roberto I Oliveira and Philip Thompson.  Sample average approximation with heavier tails i: non-asymptotic bounds with weak assumptions and stochastic constraints.  Mathematical Programming, 199(1):1–48, 2023.
  * [45] Shahar Mendelson.  Learning without concentration for general loss functions.  Probability Theory and Related Fields, 171(1):459–502, 2018.
  * [46] Shahar Mendelson.  Learning without concentration.  Journal of the ACM (JACM), 62(3):1–25, 2015.
  * [47] Abhishek Roy, Krishnakumar Balasubramanian, and Murat A Erdogdu.  On empirical risk minimization with dependent and heavy-tailed data.  Advances in Neural Information Processing Systems, 34:8913–8926, 2021.
  * [48] Vu C Dinh, Lam S Ho, Binh Nguyen, and Duy Nguyen.  Fast learning rates with heavy-tailed losses.  Advances in neural information processing systems, 29, 2016.
  * [49] Stanislav Minsker.  Geometric median and robust estimation in banach spaces.  Bernoulli, 21(4):2308–2335, 2015.
  * [50] Gábor Lugosi and Shahar Mendelson.  Mean estimation and regression under heavy-tailed distributions: A survey.  Foundations of Computational Mathematics, 19(5):1145–1190, 2019.
  * [51] Gábor Lugosi and Shahar Mendelson.  Sub-Gaussian estimators of the mean of a random vector.  The Annals of Statistics, 47(2):783–794, 2019.
  * [52] Guillaume Lecué and Matthieu Lerasle.  Learning from mom’s principles: Le cam’s approach.  Stochastic Processes and their applications, 129(11):4385–4410, 2019.
  * [53] Gabor Lugosi and Shahar Mendelson.  Risk minimization by median-of-means tournaments.  Journal of the European Mathematical Society, 22(3):925–965, 2019.
  * [54] Guillaume Lecué and Matthieu Lerasle.  Robust machine learning by median-of-means: Theory and practice.  The Annals of Statistics, 48(2):906–931, 2020.
  * [55] Daniel Hsu and Sivan Sabato.  Loss minimization and parameter estimation with heavy tails.  Journal of Machine Learning Research, 17(18):1–40, 2016.
  * [56] Daniel Hsu and Sivan Sabato.  Heavy-tailed regression with a generalized median-of-means.  In International Conference on Machine Learning, pages 37–45. PMLR, 2014.
  * [57] Jean-Yves Audibert and Olivier Catoni.  Robust linear least squares regression.  The Annals of Statistics, 39(5):2766–2794, 2011.
  * [58] Lijun Zhang and Zhi-Hua Zhou.  $\ell_{1}$-regression with heavy-tailed distributions.  Advances in Neural Information Processing Systems, 31, 2018.
  * [59] Christian Brownlees, Edouard Joly, and Gábor Lugosi.  Empirical risk minimization for heavy-tailed losses.  The Annals of Statistics, 43(6):2507–2536, 2015.
  * [60] Eduard Gorbunov, Marina Danilova, and Alexander Gasnikov.  Stochastic optimization with heavy-tailed noise via accelerated gradient clipping.  Advances in Neural Information Processing Systems, 33:15042–15053, 2020.
  * [61] Nikita Puchkin, Eduard Gorbunov, Nickolay Kutuzov, and Alexander Gasnikov.  Breaking the heavy-tailed noise barrier in stochastic optimization problems.  In International Conference on Artificial Intelligence and Statistics, pages 856–864. PMLR, 2024.
  * [62] Arnaud Deza and Elias B Khalil.  Machine learning for cutting planes in integer programming: a survey.  In Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence, pages 6592–6600, 2023.
  * [63] Yunhao Tang, Shipra Agrawal, and Yuri Faenza.  Reinforcement learning for integer programming: Learning to cut.  In International conference on machine learning, pages 9367–9376. PMLR, 2020.
  * [64] Elias Khalil, Pierre Le Bodic, Le Song, George Nemhauser, and Bistra Dilkina.  Learning to branch in mixed integer programming.  Proceedings of the AAAI Conference on Artificial Intelligence, 30(1), 2016.
  * [65] He He, Hal Daume III, and Jason M Eisner.  Learning to search in branch and bound algorithms.  Advances in neural information processing systems, 27, 2014.
  * [66] Lara Scavuzzo, Feng Chen, Didier Chételat, Maxime Gasse, Andrea Lodi, Neil Yorke-Smith, and Karen Aardal.  Learning to branch with tree mdps.  Advances in Neural Information Processing Systems, 35:18514–18526, 2022.
  * [67] Timo Berthold and Gregor Hendel.  Learning to scale mixed-integer programs.  Proceedings of the AAAI Conference on Artificial Intelligence, 35(5):3661–3668, 2021.
  * [68] Yunzhuang Shen, Yuan Sun, Andrew Eberhard, and Xiaodong Li.  Learning primal heuristics for mixed integer programs.  In 2021 international joint conference on neural networks (ijcnn), pages 1–8. IEEE, 2021.
  * [69] Irwan Bello, Hieu Pham, Quoc V Le, Mohammad Norouzi, and Samy Bengio.  Neural combinatorial optimization with reinforcement learning.  arXiv preprint arXiv:1611.09940, 2016.
  * [70] Tianlong Chen, Xiaohan Chen, Wuyang Chen, Howard Heaton, Jialin Liu, Zhangyang Wang, and Wotao Yin.  Learning to optimize: A primer and a benchmark.  Journal of Machine Learning Research, 23(189):1–59, 2022.
  * [71] Xiaohan Chen, Jialin Liu, and Wotao Yin.  Learning to optimize: A tutorial for continuous and mixed-integer optimization.  Science China Mathematics, pages 1–72, 2024.
  * [72] Yoshua Bengio, Andrea Lodi, and Antoine Prouvost.  Machine learning for combinatorial optimization: a methodological tour d’horizon.  European Journal of Operational Research, 290(2):405–421, 2021.
  * [73] Jiayi Zhang, Chang Liu, Xijun Li, Hui-Ling Zhen, Mingxuan Yuan, Yawen Li, and Junchi Yan.  A survey for solving mixed integer programming via machine learning.  Neurocomputing, 519:205–217, 2023.
  * [74] Wassily Hoeffding.  Probability inequalities for sums of bounded random variables.  Journal of the American Statistical Association, 58(301):13–30, 1963.
  * [75] Anton J Kleywegt, Alexander Shapiro, and Tito Homem-de Mello.  The sample average approximation method for stochastic discrete optimization.  SIAM Journal on optimization, 12(2):479–502, 2002.
  * [76] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267–288, 1996.

Supplemental materials

The appendices are organized as follows. In Appendix A, we provide a
literature review for the related work. Appendix B presents additional
discussion and technical results for Algorithm 1. Next, in Appendix C, we
document the proofs of theoretical results in our paper. Specifically, we
introduce some preliminary definitions and lemmas in Appendix C.1. Then, the
proofs of Theorem 1 and Proposition 1 can be found in Appendix C.2, and the
proof of Theorem 2 can be found in Appendix C.3. Finally, we provide
additional numerical experiments in Appendix D.

##  Appendix A Related work

This work is closely connected to various topics in optimization and machine
learning, and we only review the three most relevant research directions.

Bagging for stochastic optimization: Bagging has been adopted in stochastic
optimization for various purposes. The most relevant line of works [26, 27,
28, 29] 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. [30] empirically evaluates
several statistical techniques including bagging against the plain SAA and
finds bagging advantageous for portfolio optimization problems. [31]
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 [32, 33, 34, 35, 36]
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 [37], randomized Newton methods [38], and
stochastic gradient descent [39, 40], which is orthogonal to our focus on
generalization performance.

Optimization and learning with heavy tails: Optimization with heavy-tailed
noises has recently attracted significant attention due to its prevalence in
conventional applications such as portfolio [18] and scheduling [41], as well
as emerging domains like large language models [42, 43]. Tail bounds of most
existing algorithms in optimization and learning are guaranteed to decay
exponentially under sub-Gaussian or uniformly bounded costs, but deteriorate
to a slow polynomial decay under heavy-tailedness [12, 13, 14, 44]. For SAA or
empirical risk minimization, faster rates are possible under the small-ball
[45, 46, 47] or Bernstein’s condition [48] on the involved function class,
while our framework is free from such conditions. Considerable effort has been
made to mitigate the adverse effects of heavy-tailedness with more
sophisticated and robust procedures than SAA, among which the geometric median
[49], or more generally, median-of-means (MOM) [50, 51] approach is most
similar to our framework. The basic idea there is to estimate a true mean by
dividing the data into disjoint subsamples, computing an estimate based on
each subsample, and then taking the median of the estimates. [52, 53, 54] use
MOM in estimating the expected cost and establish exponential tail bounds for
the mean squared loss and convex function classes. [55, 56] 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 [22] and allows heavy
tails for linear regression [57, 58] or problems with uniformly bounded
function classes [59], but is computationally intractable due to the
truncation and thus more of theoretical interest. In contrast, our framework
is a meta algorithm that can act on any training algorithm, thus providing
exponential tail bounds for general costs and decision spaces and inheriting
the computational efficiency of the underlying training algorithm. Relatedly,
various techniques such as gradient clipping [17, 60] and MOM [61] 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., [62, 63]), search strategies (e.g., [64, 65, 66]), scaling
[67], and primal heuristics [68] in mixed integer optimization, or even
directly generate high-quality solutions (e.g., neural combinatorial
optimization pioneered by [69]). See [70, 71, 72, 73] for comprehensive
surveys on L2O. This line of research is in parallel with our goal, and L2O
techniques can work as part of or directly serve as the training algorithm
within our framework.

##  Appendix B Additional technical results

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

###### Proposition 1

Assume the condition of Theorem 1. Let $Z^{*}:=\min_{x\in\mathcal{X}}Z(x)$,
$\Delta:=\min\left\\{Z(x)-Z^{*}:Z(x)\neq Z^{*},x\in\mathcal{X}\right\\}$, and
$\mathcal{X}^{\delta}:=\left\\{x\in\mathcal{X}:Z(x)\leq
Z^{*}+\delta\right\\}$. For Algorithm 1 we have
$\eta_{k,\delta}\geq\sup_{\delta^{\prime}\geq
0}(1-T_{k}(\delta^{\prime}/2))/\lvert\mathcal{X}^{\delta^{\prime}}\rvert-
T_{k}(\delta/2)$ whenever $\sup_{\delta^{\prime}\geq
0}(1-T_{k}(\delta^{\prime}/2))/\lvert\mathcal{X}^{\delta^{\prime}}\rvert-
T_{k}(\delta/2)>0$. In particular, if (1) has a unique optimal solution and
$T_{k}(\Delta/2)+T_{k}(\delta/2)<1/5$, the bound (10) reduces to

$\mathbb{P}\left(Z(\hat{x}_{n,k}^{BAG})>\min_{x\in\mathcal{X}}Z(x)+\delta\right)\leq\left|\mathcal{X}\right|\left(3\min\left(e^{-2/5},C_{1}T_{k}\left(\frac{\Delta}{2}\right)\right)^{\frac{n}{C_{2}k}}+\exp\left(-\dfrac{B}{C_{3}}\right)\right).$

Proposition 1 shows that, in the case of multiple optimal solutions from (1),
the generalization sensitivity $\eta_{k,\delta}$ is bounded away from $1$, and
hence the bound (10) that fully exploits the generalization performance of the
training algorithm is no longer applicable. The presence of the problem gap
$\Delta$ in the bound also indicates a similar effect when there exist nearly
optimal solutions. Therefore Algorithm 1 may exhibit suboptimal overall
performance in the case of multiple or nearly optimal solutions.

##  Appendix C Technical proofs

###  C.1 Preliminaries

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

###### Definition 1

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

$U(\xi_{1},\ldots,\xi_{n})=\frac{k!}{n(n-1)\cdots(n-k+1)}\sum_{1\leq
i_{1}<i_{2}<\cdots<i_{k}\leq n}\kappa(\xi_{i_{1}},\ldots,\xi_{i_{k}}).$ (13)

######  Lemma 1 (MGF dominance of U-statistics from [74])

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

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

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(\xi_{1},\ldots,\xi_{n})=\frac{1}{n!}\sum_{\mathrm{bijection}\
\pi:[n]\rightarrow[n]}\bar{\kappa}(\xi_{\pi(1)},\ldots,\xi_{\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}(\xi_{\pi(1)},\ldots,\xi_{\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}(\xi_{\pi(1)},\ldots,\xi_{\pi(n)})\right)\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\exp\left(t\cdot\bar{\kappa}(\xi_{1},\ldots,\xi_{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(\xi_{1},\dots,\xi_{k})$ be a $\\{0,1\\}$-valued symmetric kernel
of order $k\leq n$, and $U(\xi_{1},\dots,\xi_{n})$ be the associated
U-statistic. 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}$ be the
following random variable

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

where we use the shorthand notation $\hat{n}:=\lfloor\frac{n}{k}\rfloor$.
Then, 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)$
(15)
$\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],$

where we apply the Markov inequality in $(i)$, and step $(ii)$ is due to Lemma
1. 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]$
(16)
$\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. Substitute (16) into (15), 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).$
(17)

Now, we consider minimizing $f(t)$ over $\mathbb{R}_{+}$. 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)$
(18) $\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 (18) into (17) 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)}$
(19) $\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 (19).

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\dfrac{(p-q)^{2}}{2\max_{r\in[\min(p,q),\max(p,q)]}r(1-r)},$
(20)

and that

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

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).$ (22)

Proof of Lemma 3. To prove inequality (20), we consider the case $0\leq
q<p\leq 1$ without loss of generality. We perform the second-order Taylor
expansion of $D_{\operatorname{KL}}(p\|q)$ at $q$:

$D_{\operatorname{KL}}(p\|q)=D_{\operatorname{KL}}(q\|q)+\dfrac{\partial
D_{\operatorname{KL}}(q\|q)}{\partial
p}\cdot(p-q)+\dfrac{\partial^{2}D_{\operatorname{KL}}(r\|q)}{\partial
p^{2}}\cdot\dfrac{(p-q)^{2}}{2},$ (23)

where $r\in[q,p]$. It is easy to compute that

$\dfrac{\partial D_{\operatorname{KL}}(p\|q)}{\partial
p}=\ln\dfrac{p}{q}-\ln\dfrac{1-p}{1-q},\qquad\dfrac{\partial^{2}D_{\operatorname{KL}}(p\|q)}{\partial
p^{2}}=\dfrac{1}{p(1-p)}.$

Since $D_{\operatorname{KL}}(q\|q)=\frac{\partial
D_{\operatorname{KL}}(q\|q)}{\partial p}=0$, it follows from (23) that

$D_{\operatorname{KL}}(p\|q)\geq\min_{r\in[q,p]}\left\\{\dfrac{\partial^{2}D_{\operatorname{KL}}(r\|q)}{\partial
p^{2}}\right\\}\cdot\dfrac{(p-q)^{2}}{2}=\dfrac{(p-q)^{2}}{2\max_{r\in[q,p]}r(1-r)}.$
(24)

The other case $0\leq p<q\leq 1$ of (20) can be similarly derived.

To show (21), 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 (21)
immediately.

Finally, the lower bound (22) 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 symmetric set-valued mapping
$X_{k}(\xi_{1},\ldots,\xi_{k}):\Xi^{k}\to 2^{\mathcal{X}}$ that is invariant
with respect to the permutation of its inputs. Each of our proposed bagging
algorithms attempts to solve the following probability-maximization problem
for a certain choice of $X_{k}$

$\max_{x\in\mathcal{X}}\ \hat{p}_{k}(x):=\mathbb{P}_{*}\left(x\in
X_{k}(\xi^{*}_{1},\ldots,\xi^{*}_{k})\right),$ (25)

where $\left\\{\xi^{*}_{1},\ldots,\xi^{*}_{k}\right\\}$ is resampled from the
i.i.d. data $\left\\{\xi_{1},\ldots,\xi_{n}\right\\}$ uniformly without
replacement, and $\mathbb{P}_{*}\left(\cdot\right)$ denotes the probability
with respect to the resampling randomness conditioned on the data. Note that
this problem is an empirical approximation of the problem

$\max_{x\in\mathcal{X}}\ p_{k}(x):=\mathbb{P}\left(x\in
X_{k}(\xi_{1},\ldots,\xi_{k})\right).$ (26)

The problem actually solved with a finite number of resamples is

$\max_{x\in\mathcal{X}}\
\bar{p}_{k}(x):=\frac{1}{B}\sum_{b=1}^{B}\mathbbm{1}(x\in
X_{k}(\xi_{1}^{b},\ldots,\xi_{k}^{b})).$ (27)

Specifically, Algorithm 1 solves (25) with

$X_{k}(\xi^{*}_{1},\ldots,\xi^{*}_{k})=\left\\{\hat{x}_{k}(\xi^{*}_{1},\ldots,\xi^{*}_{k})\right\\}$
(28)

where $\hat{x}_{k}(\cdot)$ denotes the solution output by the training
algorithm based on a data set of size $k$, whereas Algorithm 2 uses

$X_{k}(\xi^{*}_{1},\ldots,\xi^{*}_{k})=\left\\{x\in\mathcal{S}:\hat{Z}_{k}(x;\xi^{*}_{1},\ldots,\xi^{*}_{k})\leq\min_{x^{\prime}\in\mathcal{S}}\hat{Z}_{k}(x^{\prime};\xi^{*}_{1},\ldots,\xi^{*}_{k})+\epsilon\right\\},$
(29)

conditioned on the retrieved solution set $\mathcal{S}$. We define:

###### Definition 2

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

$\mathcal{P}^{\delta}_{k}:=\\{x\in\mathcal{X}:p_{k}(x)\geq\max_{x^{\prime}\in\mathcal{X}}p_{k}(x^{\prime})-\delta\\}$
(30)

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

$x_{k}^{\max}\in\operatorname{argmax}_{x\in\mathcal{X}}p_{k}(x)$

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}:=\\{x\in\mathcal{X}:\hat{p}_{k}(x)\geq\hat{p}_{k}(x_{k}^{\max})-\delta\\}$
(31)

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

and

###### Definition 3

Let

$\mathcal{X}^{\delta}:=\left\\{x\in\mathcal{X}:Z(x)\leq\min_{x^{\prime}\in\mathcal{X}}Z(x^{\prime})+\delta\right\\}$
(32)

be the set of $\delta$-optimal solutions of problem (1), and

$\widehat{\mathcal{X}}_{k}^{\delta}:=\left\\{x\in\mathcal{X}:\hat{Z}_{k}(x;\xi_{1},\ldots,\xi_{k})\leq\min_{x^{\prime}\in\mathcal{X}}\hat{Z}_{k}(x^{\prime};\xi_{1},\ldots,\xi_{k})+\delta\right\\}$
(33)

be the set of $\delta$-optimal solutions of (4) with i.i.d. data
$(\xi_{1},\ldots,\xi_{k})$. For every $\epsilon\geq 0$ and $\delta\geq 0$
define

$q_{k}^{\epsilon,\delta}:=\mathbb{P}\left(\widehat{\mathcal{X}}_{k}^{\epsilon}\subseteq\mathcal{X}^{\delta}\right),$
(34)

and

$r_{k}^{\epsilon}:=\mathbb{P}\left(\mathcal{X}^{0}\subseteq\widehat{\mathcal{X}}_{k}^{\epsilon}\right).$
(35)

We have the following lower bound of $q_{k}^{\epsilon,\delta}$ and
$r_{k}^{\epsilon}$:

###### Lemma 4

Assuming $0\leq\epsilon\leq\delta$. For general space $\mathcal{X}$, it holds
that

$q_{k}^{\epsilon,\delta}\geq
1-T_{k}\left(\dfrac{\delta-\epsilon}{2}\right),r_{k}^{\epsilon}\geq
1-T_{k}\left(\dfrac{\epsilon}{2}\right).$

where $T_{k}$ is the tail probability defined in (11).

Proof of Lemma 4. We have

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

therefore

$q_{k}^{\epsilon,\delta}\geq\mathbb{P}\left(\max_{x\in\mathcal{X}}\left|\hat{Z}_{k}(x)-Z(x)\right|\leq\dfrac{\delta-\epsilon}{2}\right)\geq
1-T_{k}\left(\dfrac{\delta-\epsilon}{2}\right).$

For $r_{k}^{\epsilon}$, we have

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

therefore

$r_{k}^{\epsilon}\geq\mathbb{P}\left(\max_{x\in\mathcal{X}}\left|\hat{Z}_{k}(x)-Z(x)\right|\leq\dfrac{\epsilon}{2}\right)\geq
1-T_{k}\left(\dfrac{\epsilon}{2}\right).$

This completes the proof. $\Box$

###  C.2 Proofs for Theorem 1 and Proposition 1

We consider a more general version of Algorithm 1 that operates on set
estimators, which reduces to exactly Algorithm 1 by using the set estimator
(28):

Algorithm 3 Generic Bagging with Set Estimator

1: Input: $n$ i.i.d. observations $\bm{\xi}_{1:n}=(\xi_{1},\ldots,\xi_{n})$,
positive integers $k<n$ and $B$, and a symmetric set estimator
$X_{k}:\Xi^{k}\to 2^{\mathcal{X}}$.

2: for $b=1$ to $B$ do

3: Randomly sample $\bm{\xi}_{k}^{b}=(\xi_{1}^{b},\ldots,\xi_{k}^{b})$
uniformly from $\bm{\xi}_{1:n}$ without replacement, and obtain
$X_{k}^{b}=X_{k}(\xi_{1}^{b},\ldots,\xi_{k}^{b})$

4: end for

5: Output
$\hat{x}^{BAG}_{n,k}\in\operatorname{argmax}_{x\in\mathcal{X}}\sum_{b=1}^{B}\mathbbm{1}(\hat{x}^{b}_{k}\in
X_{k}^{b})$.

We have the following finite-sample result for Algorithm 3

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

Consider discrete decision space $\mathcal{X}$. Recall $p_{k}(x)$ defined in
(26). Let $p_{k}^{\max}:=\max_{x\in\mathcal{X}}p_{k}(x)$,
$D_{\operatorname{KL}}(p\|q):=p\ln\frac{p}{q}+(1-p)\ln\frac{1-p}{1-q}$ be the
Kullback–Leibler divergence between two Bernoulli distributions with means $p$
and $q$, and

$\eta_{k,\delta}:=\max_{x\in\mathcal{X}^{\delta}}p_{k}(x)-\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x),$
(36)

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

$\displaystyle\mathbb{P}\left(Z(\hat{x}_{n,k}^{BAG})>\min_{x\in\mathcal{X}}Z(x)+\delta\right)$
(37) $\displaystyle\leq$
$\displaystyle\left|\mathcal{X}\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,\eta_{k,\delta}]$. In particular, if
$\eta_{k,\delta}>4/5$, (37) is further bounded by

$\left|\mathcal{X}\right|\left(3\min\left(e^{-2/5},C_{1}\max(1-p_{k}^{\max},\;p_{k}^{\max}-\eta_{k,\delta})\right)^{\frac{n}{C_{2}k}}+\exp\left(-\frac{B}{C_{3}}\right)\right),$
(38)

where $C_{1},C_{2},C_{3}>0$ are universal constants.

Proof of Theorem 3. We first prove generalization tail bounds for the problem
(26), splitted into two lemmas, Lemmas 5 and 6 below.

###### Lemma 5

Consider discrete decision space $\mathcal{X}$. Recall from Definition 2 that
$p_{k}^{\max}=p_{k}(x_{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|\mathcal{X}\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 5. By Definition 2, we observe the following equivalence

$\left\\{\widehat{\mathcal{P}}^{\epsilon}_{k}\not\subseteq\mathcal{P}^{\delta}_{k}\right\\}=\bigcup_{x\in\mathcal{X}\backslash\mathcal{P}^{\delta}_{k}}\left\\{x\in\widehat{\mathcal{P}}^{\epsilon}_{k}\right\\}=\bigcup_{x\in\mathcal{X}\backslash\mathcal{P}^{\delta}_{k}}\left\\{{\hat{p}}_{k}(x)\geq{\hat{p}}_{k}\left(x_{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_{x\in\mathcal{X}\backslash\mathcal{P}^{\delta}_{k}}\mathbb{P}\left(\left\\{{\hat{p}}_{k}(x)\geq{\hat{p}}_{k}\left(x_{k}^{\max}\right)-\epsilon\right\\}\right).$

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

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

On one hand, the first probability in (39) is solely determined by and
increasing in $p_{k}(x)=\mathbb{E}\left[{\hat{p}}_{k}(x)\right]$. On the other
hand, we have $p_{k}(x)<p_{k}\left(x_{k}^{\max}\right)-\delta$ for every
$x\in\mathcal{X}\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(\left\\{{\hat{p}}_{k}(x)\geq{\hat{p}}_{k}\left(x_{k}^{\max}\right)-\epsilon\right\\}\right)$
$\displaystyle\leq$ $\displaystyle\mathbb{P}\left(\left\\{{\hat{p}}_{k}(x)\geq
p_{k}(x_{k}^{\max})-\dfrac{\delta+\epsilon}{2}\right\\}\Big{|}p_{k}(x)=p_{k}(x_{k}^{\max})-\delta\right)$
$\displaystyle\hskip
8.61108pt+\mathbb{P}\left(\left\\{{\hat{p}}_{k}\left(x_{k}^{\max}\right)\leq
p_{k}(x_{k}^{\max})-\dfrac{\delta-\epsilon}{2}\right\\}\right)$
$\displaystyle\leq$
$\displaystyle\mathbb{P}\left(\left\\{{\hat{p}}_{k}(x)-p_{k}(x)\geq\dfrac{\delta-\epsilon}{2}\right\\}\Big{|}p_{k}(x)=p_{k}(x_{k}^{\max})-\delta\right)$
$\displaystyle\hskip
8.61108pt+\mathbb{P}\left(\left\\{{\hat{p}}_{k}\left(x_{k}^{\max}\right)-p_{k}(x_{k}^{\max})\leq-\dfrac{\delta-\epsilon}{2}\right\\}\right).$

Note that the probability ${\hat{p}}_{k}(x)$ can be viewed as a U-statistic
with a binary kernel

$\displaystyle\hat{p}_{k}(x)$
$\displaystyle=\frac{k!}{n(n-1)\cdots(n-k+1)}\sum_{1\leq
i_{1}<\ldots<i_{k}\leq n}\mathbbm{1}(x\in
X_{k}(\xi_{i_{1}},\ldots,\xi_{i_{k}})).$

A similar representation holds for ${\hat{p}}_{k}\left(x_{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_{x\in\mathcal{X}\backslash\mathcal{P}^{\delta}_{k}}\mathbb{P}\left(\left\\{{\hat{p}}_{k}(x)\geq{\hat{p}}_{k}\left(x_{k}^{\max}\right)-\epsilon\right\\}\right)$
$\displaystyle\leq\left|\mathcal{X}\backslash\mathcal{P}^{\delta}_{k}\right|\left[\mathbb{P}\left({\hat{p}}_{k}(x)-p_{k}(x)\geq\dfrac{\delta-\epsilon}{2}\Big{|}p_{k}(x)=p_{k}(x_{k}^{\max})-\delta\right)\right.$
$\displaystyle\hskip
51.6665pt+\left.\mathbb{P}\left(p_{k}\left(x_{k}^{\max}\right)-{\hat{p}}_{k}\left(x_{k}^{\max}\right)\leq-\dfrac{\delta-\epsilon}{2}\right)\right]$
$\displaystyle\leq\left|\mathcal{X}\right|\left[\exp\left(-\dfrac{n}{2k}\cdot
D_{\operatorname{KL}}\left(p_{k}\left(x_{k}^{\max}\right)-\delta+\dfrac{\delta-\epsilon}{2}\Big{\|}p_{k}\left(x_{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(x_{k}^{\max}\right)-\dfrac{\delta-\epsilon}{2}\Big{\|}p_{k}\left(x_{k}^{\max}\right)\right)\right)\right],$

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

###### Lemma 6

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

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

where $\left|\cdot\right|$ denotes the cardinality of a set.

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

$\bar{p}_{k}(x)-\bar{p}_{k}(x_{k}^{\max})=\dfrac{1}{B}\sum_{b=1}^{B}\left[\mathbbm{1}(x\in
X_{k}(\xi_{1}^{b},\ldots,\xi_{k}^{b}))-\mathbbm{1}(x_{k}^{\max}\in
X_{k}(\xi_{1}^{b},\ldots,\xi_{k}^{b}))\right],$

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

$\mathbbm{1}_{x}^{*}:=\mathbbm{1}(x\in
X_{k}(\xi_{1}^{*},\ldots,\xi_{k}^{*}))\text{\ for }x\in\mathcal{X}$

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

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

Since

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

and

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

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

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

Therefore, we have that

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

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

$\displaystyle\mathbb{P}_{*}\left(\hat{x}_{n,k}^{BAG}\notin\widehat{\mathcal{P}}^{\epsilon}_{k}\right)$
$\displaystyle\leq$
$\displaystyle\left|\mathcal{X}\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|\mathcal{X}\backslash\widehat{\mathcal{P}}^{\epsilon}_{k}\right|\leq\left|\mathcal{X}\right|$
completes the proof of Lemma 6. $\Box$

We are now ready for the proof of Theorem 3. If $\eta_{k,\delta}>0$, it
follows from Definition 2 that

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

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

$\displaystyle\mathbb{P}\left(\hat{x}_{n,k}^{BAG}\notin\mathcal{X}^{\delta}\right)\leq\mathbb{P}\left(\hat{x}_{n,k}^{BAG}\notin\mathcal{P}^{\eta}_{k}\right)$
$\displaystyle\leq$
$\displaystyle\mathbb{P}\left(\left\\{\hat{x}_{n,k}^{BAG}\notin\widehat{\mathcal{P}}^{\eta/2}_{k}\right\\}\cup\left\\{\widehat{\mathcal{P}}^{\eta/2}_{k}\not\subseteq\mathcal{P}^{\eta}_{k}\right\\}\right)$
(41) $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(\hat{x}_{n,k}^{BAG}\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 (41). Lemma
5 gives that

$\displaystyle\mathbb{P}\left(\widehat{\mathcal{P}}^{\eta/2}_{k}\not\subseteq\mathcal{P}^{\eta}_{k}\right)\leq\left|\mathcal{X}\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)$
(42) $\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 6 with $\epsilon=\eta/2$, we can bound the first
probability on the right-hand side of (41) as

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

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

$\displaystyle\mathbb{E}\left[\exp\left(-\dfrac{B}{24}\cdot\dfrac{\eta^{2}}{\min\left(\hat{p}_{k}(x_{k}^{\max}),1-\hat{p}_{k}(x_{k}^{\max})\right)+\eta/2}\right)\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{P}\left(\hat{p}_{k}(x_{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}(x_{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}(x_{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}(x_{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}(x_{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}(x_{k}^{\max})$, the U-statistic estimate for $p_{k}^{\max}$.
Together, the above equations imply that

$\displaystyle\mathbb{P}\left(\hat{x}_{n,k}^{BAG}\notin\mathcal{X}^{\delta}\right)$
$\displaystyle\leq$
$\displaystyle\left|\mathcal{X}\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,\;\eta_{k,\delta})$, the validity of the case
$\eta=\eta_{k,\delta}$ follows from pushing $\eta$ to the limit
$\eta_{k,\delta}$. This gives (37).

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

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

and that $p_{k}^{\max}-\eta_{k,\delta}\leq 1-\eta_{k,\delta}\leq 1/5$,
therefore by the bound (22) 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\eta_{k,\delta}}{4}\Big{\|}p_{k}^{\max}-\eta_{k,\delta}\right)\right)$
$\displaystyle\leq$
$\displaystyle\exp\left(\dfrac{n}{2k}\ln\left(2((p_{k}^{\max}-\eta_{k,\delta})(1-p_{k}^{\max}+\eta_{k,\delta}))^{1/5}\right)\right)$
$\displaystyle=$
$\displaystyle\left(2((p_{k}^{\max}-\eta_{k,\delta})(1-p_{k}^{\max}+\eta_{k,\delta}))^{1/5}\right)^{n/(2k)}$
$\displaystyle\leq$
$\displaystyle\left(2(p_{k}^{\max}-\eta_{k,\delta})^{1/5}\right)^{n/(2k)}$
$\displaystyle=$
$\displaystyle\left(2^{5}(p_{k}^{\max}-\eta_{k,\delta})\right)^{n/(10k)},$

and similarly

$\displaystyle\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\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\eta_{k,\delta}/4\|p_{k}^{\max}-\eta_{k,\delta}\right)$
and
$D_{\operatorname{KL}}\left(p_{k}^{\max}-\eta_{k,\delta}/4\|p_{k}^{\max}\right)$
are bounded below by $\eta_{k,\delta}^{2}/8$, therefore

$\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)\leq\exp\left(-\dfrac{n}{2k}\cdot\dfrac{\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}-\eta_{k,\delta}/4\|p_{k}^{\max}\right)\right)$.
For the third term in the bound (37) we have

$\dfrac{\eta_{k,\delta}^{2}}{\min(p_{k}^{\max},1-p_{k}^{\max})+3\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{\eta_{k,\delta}^{2}}{\min\left(p_{k}^{\max},1-p_{k}^{\max}\right)+3\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 set
estimator (28), therefore the results of Theorem 3 automatically apply. $\Box$

######  Lemma 7 (Bound of $\eta_{k,\delta}$ for set estimator (28))

Consider discrete decision space $\mathcal{X}$. If the set estimator (28) is
used, it holds that

$p_{k}^{\max}=\max_{x\in\mathcal{X}}p_{k}(x)\geq\sup_{\delta\geq
0}\frac{q_{k}^{0,\delta}}{\lvert\mathcal{X}^{\delta}\rvert},$ (44)

and that

$\displaystyle\eta_{k,\delta}$ $\displaystyle\geq\sup_{\delta^{\prime}\geq
0}\frac{q_{k}^{0,\delta^{\prime}}}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}+q_{k}^{0,\delta}-1\quad\text{whenever
}\sup_{\delta^{\prime}\geq
0}\frac{q_{k}^{0,\delta^{\prime}}}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}+q_{k}^{0,\delta}-1>0,$
(45) $\displaystyle\eta_{k,\delta}$ $\displaystyle\geq
2p_{k}^{\max}-1\quad\text{whenever }\eta_{k,\delta}>0.$

Proof of Lemma 7. For any $\delta\geq 0$, we can write that

$\displaystyle\max_{x\in\mathcal{X}}p_{k}(x)$ $\displaystyle\geq$
$\displaystyle\max_{x\in\mathcal{X}^{\delta}}p_{k}(x)$ $\displaystyle\geq$
$\displaystyle\frac{\sum_{x\in\mathcal{X}^{\delta}}p_{k}(x)}{\lvert\mathcal{X}^{\delta}\rvert}$
$\displaystyle=$
$\displaystyle\frac{\mathbb{P}(\hat{x}_{k}(\xi_{1},\ldots,\xi_{k})\in\mathcal{X}^{\delta})}{\lvert\mathcal{X}^{\delta}\rvert}$
$\displaystyle\geq$
$\displaystyle\frac{\mathbb{P}(\widehat{\mathcal{X}}_{k}^{0}\subseteq\mathcal{X}^{\delta})}{\lvert\mathcal{X}^{\delta}\rvert}=\frac{q_{k}^{0,\delta}}{\lvert\mathcal{X}^{\delta}\rvert},$

which directly implies the first inequality in the lemma, since $\delta$ can
be any non-negative number. For a fixed $\delta\geq 0$, it holds that for any
$x\in\mathcal{X}\backslash\mathcal{X}^{\delta}$

$p_{k}(x)=\mathbb{P}(\hat{x}_{k}(\xi_{1},\ldots,\xi_{k})=x)\leq\mathbb{P}(\widehat{\mathcal{X}}_{k}^{0}\not\subseteq\mathcal{X}^{\delta})=1-q_{k}^{0,\delta},$

therefore $\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)\leq
1-q_{k}^{0,\delta}$. Combining this upper bound for
$\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)$ and the bound
(44) gives

$\max_{x\in\mathcal{X}}p_{k}(x)-\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)\geq\sup_{\delta^{\prime}\geq
0}\frac{q_{k}^{0,\delta^{\prime}}}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}+q_{k}^{0,\delta}-1.$

Whenever $\sup_{\delta^{\prime}\geq
0}\frac{q_{k}^{0,\delta^{\prime}}}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}+q_{k}^{0,\delta}-1>0$
we have
$\max_{x\in\mathcal{X}}p_{k}(x)>\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)$
and hence
$\max_{x\in\mathcal{X}}p_{k}(x)=\max_{x\in\mathcal{X}^{\delta}}p_{k}(x)$,
therefore
$\eta_{k,\delta}=\max_{x\in\mathcal{X}}p_{k}(x)-\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)\geq\sup_{\delta^{\prime}\geq
0}\frac{q_{k}^{0,\delta^{\prime}}}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}+q_{k}^{0,\delta}-1$.
To show the second bound in (45), note that $\eta_{k,\delta}>0$ implies that
$p_{k}^{\max}$ is achieved by some $x^{\prime}\in\mathcal{X}^{\delta}$. On the
other hand, the particular choice of the indicator (28) ensures that
$\sum_{x\in\mathcal{X}}p_{k}(x)=1$. Thus we can write

$p_{k}^{\max}-\eta_{k,\delta}\leq\sum_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)\leq\sum_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)+p_{k}(x^{\prime})-p_{k}(x^{\prime})\leq
1-p_{k}(x^{\prime})=1-p_{k}^{\max},$

which immediately implies the desired bound. $\Box$

Proof of Proposition 1. By Lemma 4 we have

$\sup_{\delta^{\prime}\geq
0}\frac{q_{k}^{0,\delta^{\prime}}}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}+q_{k}^{0,\delta}-1\geq\sup_{\delta^{\prime}\geq
0}\frac{1-T_{k}(\delta^{\prime}/2)}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}-T_{k}(\delta/2),$

therefore whenever $\sup_{\delta^{\prime}\geq
0}\frac{1-T_{k}(\delta^{\prime}/2)}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}-T_{k}(\delta/2)>0$,
we have $\sup_{\delta^{\prime}\geq
0}\frac{q_{k}^{0,\delta^{\prime}}}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}+q_{k}^{0,\delta}-1>0$,
and hence by Lemma 7

$\eta_{k,\delta}\geq\sup_{\delta^{\prime}\geq
0}\frac{q_{k}^{0,\delta^{\prime}}}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}+q_{k}^{0,\delta}-1\geq\sup_{\delta^{\prime}\geq
0}\frac{1-T_{k}(\delta^{\prime}/2)}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}-T_{k}(\delta/2).$

This proves the first part. If (1) has a unique solution, then
$\lvert\mathcal{X}^{\delta^{\prime}}\rvert=1$ for all
$\delta^{\prime}<\Delta$, therefore

$\sup_{\delta^{\prime}\geq
0}\frac{1-T_{k}(\delta^{\prime}/2)}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}\geq\lim_{\delta^{\prime}\to\Delta-}\frac{1-T_{k}(\delta^{\prime}/2)}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}=1-T_{k}(\Delta/2).$

When $T_{k}(\Delta/2)+T_{k}(\delta/2)<1/5$, i.e.,
$1-T_{k}(\Delta/2)-T_{k}(\delta/2)>4/5$, we also have $\eta_{k,\delta}>4/5$,
therefore the bound (10) applies. The desired bound then follows by noticing
that $p_{k}^{\max}-\eta_{k,\delta}\leq
p_{k}^{\max}-(2p_{k}^{\max}-1)=1-p_{k}^{\max}$ by Lemma 7 and that
$1-p_{k}^{\max}\leq 1-\sup_{\delta^{\prime}\geq
0}\frac{q_{k}^{0,\delta^{\prime}}}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}\leq
1-\sup_{\delta^{\prime}\geq
0}\frac{1-T_{k}(\delta^{\prime}/2)}{\lvert\mathcal{X}^{\delta^{\prime}}\rvert}\leq
T_{k}(\Delta/2)$ where the first inequality is by (44). $\Box$

###  C.3 Proof for Theorem 2

######  Lemma 8 (Gap translation for set estimator (29))

Consider discrete decision space $\mathcal{X}$. If the estimator (29) is used
with $\epsilon\geq 0$, it holds that

$p_{k}^{\max}=\max_{x\in\mathcal{X}}p_{k}(x)\geq r_{k}^{\epsilon}.$ (46)

For any $\delta\geq 0$, whenever
$r_{k}^{\epsilon}+q_{k}^{\epsilon,\delta}-1>0$, it holds that

$\eta_{k,\delta}\geq r_{k}^{\epsilon}+q_{k}^{\epsilon,\delta}-1,$ (47)

and that

$p_{k}^{\max}-\eta_{k,\delta}\leq 1-q_{k}^{\epsilon,\delta}.$ (48)

Proof of Lemma 8. Let $x^{*}$ be an optimal solution of (1). By the definition
of $r_{k}^{\epsilon}$ we have

$\max_{x\in\mathcal{X}}p_{k}(x)\geq
p_{k}(x^{*})=\mathbb{P}\left(x^{*}\in\widehat{\mathcal{X}}_{k}^{\epsilon}\right)\geq\mathbb{P}\left(\mathcal{X}^{0}\in\widehat{\mathcal{X}}_{k}^{\epsilon}\right)=r_{k}^{\epsilon}.$

This proves (46). For any $x\in\mathcal{X}\backslash\mathcal{X}^{\delta}$ it
holds that

$p_{k}(x)=\mathbb{P}\left(x\in\widehat{\mathcal{X}}_{k}^{\epsilon}\right)\leq\mathbb{P}\left(\widehat{\mathcal{X}}_{k}^{\epsilon}\not\subseteq\mathcal{X}^{\delta}\right)=1-q_{k}^{\epsilon,\delta},$

hence $\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)\leq
1-q_{k}^{\epsilon,\delta}$. Therefore, whenever
$r_{k}^{\epsilon}-(1-q_{k}^{\epsilon,\delta})>0$, we have
$\max_{x\in\mathcal{X}}p_{k}(x)>\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)$
and thus
$\max_{x\in\mathcal{X}}p_{k}(x)=\max_{x\in\mathcal{X}^{\delta}}p_{k}(x)$,
i.e.,
$\eta_{k,\delta}=\max_{x\in\mathcal{X}}p_{k}(x)-\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)\geq
r_{k}^{\epsilon}-(1-q_{k}^{\epsilon,\delta})$. This proves (47). To show (48),
note that

$p_{k}^{\max}-\eta_{k,\delta}=\max_{x\in\mathcal{X}^{\delta}}p_{k}(x)-\eta_{k,\delta}=\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x),$

and that for each $x\in\mathcal{X}\backslash\mathcal{X}^{\delta}$ it holds
that
$p_{k}(x)=\mathbb{P}\left(x\in\widehat{\mathcal{X}}_{k}^{\epsilon}\right)\leq\mathbb{P}\left(\widehat{\mathcal{X}}_{k}^{\epsilon}\not\subseteq\mathcal{X}^{\delta}\right)=1-q_{k}^{\epsilon,\delta}$.
Therefore $p_{k}^{\max}-\eta_{k,\delta}\leq 1-q_{k}^{\epsilon,\delta}$. This
completes the proof. $\Box$

######  Lemma 9 (Near optimality of retrieved solutions in Algorithms 2)

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

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

where $q_{k}^{0,\delta}$ is from (34), and $C_{4},C_{5},C_{6},C_{7}>0$ are
universal constants. The same holds true for Algorithm 2 without data
splitting if $n$ is replaced by $n/2$.

Proof of Lemma 9. Denote by
$\hat{x}_{k}^{*}:=\hat{x}_{k}(\xi_{1}^{*},\ldots,\xi_{k}^{*})$ an optimal
solution output by the training algorithm on a resampled data set of size $k$,
and by $\hat{x}_{k}:=\hat{x}_{k}(\xi_{1},\ldots,\xi_{k})$ an optimal solution
output by the training algorithm based on an i.i.d. data set of size $k$. Let
$\mathbb{P}_{*}$ denote the probability with respect to the resampling
randomness conditioned on the data. Consider the two probabilities

$\mathbb{P}\left(\hat{x}_{k}\in\mathcal{X}^{\delta}\right),\
\mathbb{P}_{*}\left(\hat{x}_{k}^{*}\in\mathcal{X}^{\delta}\right).$

We have

$p_{k,\delta}:=\mathbb{P}\left(\hat{x}_{k}\in\mathcal{X}^{\delta}\right)\geq\mathbb{P}\left(\widehat{\mathcal{X}}_{k}^{0}\subseteq\mathcal{X}^{\delta}\right)=q_{k}^{0,\delta}$

by definition of $q_{k}^{0,\delta}$. We also have the conditional probability

$\mathbb{P}\left(\mathcal{S}\cap\mathcal{X}^{\delta}=\emptyset\Big{\rvert}\mathbb{P}_{*}\left(\hat{x}_{k}^{*}\in\mathcal{X}^{\delta}\right)\right)=\left(1-\mathbb{P}_{*}\left(\hat{x}_{k}^{*}\in\mathcal{X}^{\delta}\right)\right)^{B_{1}}.$

Therefore we can write

$\displaystyle\mathbb{P}\left(\mathcal{S}\cap\mathcal{X}^{\delta}=\emptyset\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(1-\mathbb{P}_{*}\left(\hat{x}_{k}^{*}\in\mathcal{X}^{\delta}\right)\right)^{B_{1}}\right]$
(50) $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(\mathbb{P}_{*}\left(\hat{x}_{k}^{*}\in\mathcal{X}^{\delta}\right)<\frac{p_{k,\delta}}{e}\right)+\left(1-\frac{p_{k,\delta}}{e}\right)^{B_{1}}$

where $e$ is the base of the natural logarithm. Applying Lemma 2 with
$\kappa(\xi_{1},\ldots,\xi_{k})=\mathbbm{1}\left(\hat{x}_{k}(\xi_{1},\ldots,\xi_{k})\in\mathcal{X}^{\delta}\right)$
gives

$\mathbb{P}\left(\mathbb{P}_{*}\left(\hat{x}_{k}^{*}\in\mathcal{X}^{\delta}\right)<\frac{p_{k,\delta}}{e}\right)\leq\exp\left(-\frac{n}{2k}\cdot
D_{\operatorname{KL}}\left(\frac{p_{k,\delta}}{e}\Big{\|}p_{k,\delta}\right)\right).$

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

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

and

$\displaystyle
D_{\operatorname{KL}}\left(\frac{p_{k,\delta}}{e}\Big{\|}p_{k,\delta}\right)$
$\displaystyle=$ $\displaystyle
D_{\operatorname{KL}}\left(1-\frac{p_{k,\delta}}{e}\Big{\|}1-p_{k,\delta}\right)$
$\displaystyle\geq$
$\displaystyle\left(1-\frac{p_{k,\delta}}{e}\right)\ln\frac{1-p_{k,\delta}/e}{1-p_{k,\delta}}-p_{k,\delta}+\frac{p_{k,\delta}}{e}\text{\
\ by bound \eqref{KL lower bound: ratio}}$ $\displaystyle\geq$
$\displaystyle\left(1-\frac{p_{k,\delta}}{e}\right)\ln\left(1-\frac{p_{k,\delta}}{e}\right)-\left(1-\frac{1}{e}\right)\ln\left(1-p_{k,\delta}\right)-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\left(1-p_{k,\delta}\right)-1+\frac{1}{e}$
$\displaystyle=$
$\displaystyle\left(1-\frac{1}{e}\right)\ln\frac{e-1}{e^{2}(1-p_{k,\delta})}.$

Combining the two bounds for the KL divergence we have

$\mathbb{P}\left(\mathbb{P}_{*}\left(\hat{x}_{k}^{*}\in\mathcal{X}^{\delta}\right)<\frac{p_{k,\delta}}{e}\right)\leq\min\left(\exp\left(-\frac{n}{2k}\cdot\left(1-\frac{2}{e}\right)p_{k,\delta}\right),\left(\frac{e^{2}(1-p_{k,\delta})}{e-1}\right)^{(1-1/e)\frac{n}{2k}}\right).$

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

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

The conclusion then follows by setting $C_{4},C_{5},C_{6},C_{7}$ to be the
appropriate constants, and noticing that $p_{k,\delta}\geq q_{k}^{0,\delta}$
and that the bound above decreases in $p_{k,\delta}$. $\Box$

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

$\mathcal{W}^{\delta}:=\left\\{x\in\mathcal{W}:Z(x)\leq\min_{x^{\prime}\in\mathcal{W}}Z(x^{\prime})+\delta\right\\}$
(51)

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

$\widehat{\mathcal{W}}_{k}^{\delta}:=\left\\{x\in\mathcal{W}:\hat{Z}_{k}(x;\xi_{1},\ldots,\xi_{k})=\min_{x^{\prime}\in\mathcal{W}}\hat{Z}_{k}(x^{\prime};\xi_{1},\ldots,\xi_{k})+\delta\right\\}$
(52)

be the set of $\delta$-optimal solutions of the restricted training
optimization problem with an i.i.d. data set of size $k$. For every
$\epsilon\geq 0$ and $\delta\geq 0$ let

$q_{k}^{\epsilon,\delta,\mathcal{W}}:=\mathbb{P}\left(\widehat{\mathcal{W}}_{k}^{\epsilon}\subseteq\mathcal{W}^{\delta}\right),$
(53)

and

$r_{k}^{\epsilon.\mathcal{W}}:=\mathbb{P}\left(\mathcal{W}^{0}\subseteq\widehat{\mathcal{W}}_{k}^{\epsilon}\right),$
(54)

be the counterparts of $q_{k}^{\epsilon,\delta}$ and $r_{k}^{\epsilon}$
respectively.

Proof of Theorem 2 for ReBAG-S. 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_{x\in\mathcal{S}}\mathbb{E}\left[h(x,\xi)\right]$ to obtain
$\hat{x}_{n,k}^{BAG}$ with the data $\bm{\xi}_{\lfloor n/2\rfloor+1:n}$, the
set estimator (29), the chosen $\epsilon$ value and $B=B_{2}$.

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

$\left\\{\mathcal{S}\cap\mathcal{X}^{\delta}\neq\emptyset\right\\}\cap\left\\{Z(\hat{x}_{n,k}^{BAG})\leq\min_{x\in\mathcal{S}}Z(x)+\delta\right\\}\subseteq\left\\{\hat{x}_{n,k}^{BAG}\in\mathcal{X}^{2\delta}\right\\},$

and hence by union bound we can write

$\mathbb{P}\left(\hat{x}_{n,k}^{BAG}\notin\mathcal{X}^{2\delta}\right)\leq\mathbb{P}\left(\mathcal{S}\cap\mathcal{X}^{\delta}=\emptyset\right)+\mathbb{P}\left(Z(\hat{x}_{n,k}^{BAG})>\min_{x\in\mathcal{S}}Z(x)+\delta\right).$
(55)

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

Given a retrieved solution set $\mathcal{S}$, we derive lower bounds for
probabilities $q_{k}^{\epsilon,\delta,\mathcal{S}}$ and
$r_{k}^{\epsilon,\mathcal{S}}$. For $\epsilon\leq\delta$, we have by applying
Lemma (4) to the restricted problem that

$q_{k}^{\epsilon,\delta,\mathcal{S}}\geq\mathbb{P}\left(\max_{x\in\mathcal{S}}\left|\hat{Z}_{k}(x)-Z(x)\right|<\dfrac{\delta-\epsilon}{2}\right)\geq\mathbb{P}\left(\max_{x\in\mathcal{X}}\left|\hat{Z}_{k}(x)-Z(x)\right|<\dfrac{\delta-\epsilon}{2}\right)=1-T_{k}\left(\dfrac{\delta-\epsilon}{2}\right).$

where in the probability
$\mathbb{P}\left(\max_{x\in\mathcal{S}}\left|\hat{Z}_{k}(x)-Z(x)\right|<\dfrac{\delta-\epsilon}{2}\right)$
the set $\mathcal{S}$ is viewed as given and fixed. For
$r_{k}^{\epsilon,\mathcal{S}}$, we similarly have

$r_{k}^{\epsilon,\mathcal{S}}\geq\mathbb{P}\left(\max_{x\in\mathcal{S}}\left|\hat{Z}_{k}(x)-Z(x)\right|<\dfrac{\epsilon}{2}\right)\geq\mathbb{P}\left(\max_{x\in\mathcal{X}}\left|\hat{Z}_{k}(x)-Z(x)\right|<\dfrac{\epsilon}{2}\right)=1-T_{k}\left(\dfrac{\epsilon}{2}\right),$

where in the probability
$\mathbb{P}\left(\max_{x\in\mathcal{S}}\left|\hat{Z}_{k}(x)-Z(x)\right|<\dfrac{\epsilon}{2}\right)$
the set $\mathcal{S}$ is viewed as given and fixed. Since
$\mathbb{P}\left(\epsilon\in[\underline{\epsilon},\overline{\epsilon}]\right)=1$,
we have

$q_{k}^{\epsilon,\delta,\mathcal{S}}\geq
q_{k}^{\overline{\epsilon},\delta,\mathcal{S}}\geq
1-T_{k}\left(\dfrac{\delta-\overline{\epsilon}}{2}\right),$ (56)

and

$r_{k}^{\epsilon,\mathcal{S}}\geq r_{k}^{\underline{\epsilon},\mathcal{S}}\geq
1-T_{k}\left(\dfrac{\underline{\epsilon}}{2}\right).$ (57)

If
$T_{k}\left((\delta-\overline{\epsilon})/2\right)+T_{k}\left(\underline{\epsilon}/2\right)<1/5$,
we have
$r_{k}^{\epsilon,\mathcal{S}}+q_{k}^{\epsilon,\delta,\mathcal{S}}-1>4/5$ and
hence $\eta_{k,\eta}>4/5$ by Lemma 8, therefore the bound (38) from Theorem
(3) applies. Using the inequalities (46) and (48) to handle the
$\min(1-p_{k}^{\max},p_{k}^{\max}-\eta_{k,\delta})$ term in (38) gives

$\displaystyle\mathbb{P}\left(Z(\hat{x}_{n,k}^{BAG})>\min_{x\in\mathcal{S}}Z(x)+\delta\big{|}\mathcal{S},\epsilon\right)$
$\displaystyle\leq$
$\displaystyle\left|\mathcal{S}\right|\left(3\min\left(e^{-2/5},C_{1}\left(1-\min(r_{k}^{\epsilon,\mathcal{S}},q_{k}^{\epsilon,\delta,\mathcal{S}})\right)\right)^{\frac{n}{2C_{2}k}}+\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}\max\left(T_{k}\left(\dfrac{\underline{\epsilon}}{2}\right),T_{k}\left(\dfrac{\delta-\overline{\epsilon}}{2}\right)\right)\right)^{\frac{n}{2C_{2}k}}+\exp\left(-\dfrac{B_{2}}{C_{3}}\right)\right)\text{\
\ by \eqref{eq: uniform lower bound for q} and \eqref{eq: uniform lower bound
for r}}$ $\displaystyle=$
$\displaystyle\left|\mathcal{S}\right|\left(3\min\left(e^{-2/5},C_{1}T_{k}\left(\dfrac{\min(\underline{\epsilon},\delta-\overline{\epsilon})}{2}\right)\right)^{\frac{n}{2C_{2}k}}+\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(Z(\hat{x}_{n,k}^{BAG})>\min_{x\in\mathcal{S}}Z(x)+\delta\right)\leq
B_{1}\left(3\min\left(e^{-2/5},C_{1}T_{k}\left(\dfrac{\min(\underline{\epsilon},\delta-\overline{\epsilon})}{2}\right)\right)^{\frac{n}{2C_{2}k}}+\exp\left(-\dfrac{B_{2}}{C_{3}}\right)\right).$

This leads to the desired bound (12) after the above bound is plugged into
(55) and the bound (49) from Lemma 9 is applied with $q_{k}^{0,\delta}$
replaced by its lower bound $1-T_{k}(\delta/2)$. $\Box$

Proof of Theorem 2 for ReBAG. For every non-empty subset
$\mathcal{W}\subseteq\mathcal{X}$ and $k$, we consider the indicator

$\mathbbm{1}_{k}^{x,\mathcal{W},\epsilon}(\xi_{1},\ldots,\xi_{k}):=\mathbbm{1}\left(\hat{Z}_{k}(x;\xi_{1},\ldots,\xi_{k})\leq\min_{x^{\prime}\in\mathcal{W}}\hat{Z}_{k}(x^{\prime};\xi_{1},\ldots,\xi_{k})+\epsilon\right)\quad\text{for
}x\in\mathcal{W},\epsilon\in[0,\delta/2],$

which indicates whether a solution $x\in\mathcal{W}$ is $\epsilon$-optimal for
the training problem formed by $\left\\{\xi_{1},\ldots,\xi_{k}\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}^{x,\mathcal{W},\epsilon}$
are

$\displaystyle p_{k}^{\mathcal{W},\epsilon}(x)$
$\displaystyle:=\mathbb{E}\left[\mathbbm{1}_{k}^{x,\mathcal{W},\epsilon}(\xi_{1},\ldots,\xi_{k})\right],$
$\displaystyle\hat{p}_{k}^{\mathcal{W},\epsilon}(x)$
$\displaystyle:=\mathbb{E}_{*}\left[\mathbbm{1}_{k}^{x,\mathcal{W},\epsilon}(\xi_{1}^{*},\ldots,\xi_{k}^{*})\right],$
$\displaystyle\bar{p}_{k}^{\mathcal{W},\epsilon}(x)$
$\displaystyle:=\dfrac{1}{B_{2}}\sum_{b=1}^{B_{2}}\mathbbm{1}_{k}^{x,\mathcal{W},\epsilon}(\xi_{1}^{b},\ldots,\xi_{k}^{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}^{x,\mathcal{W},\epsilon}$

$\mathbbm{1}_{k}^{x,\mathcal{W},\epsilon-}(\xi_{1},\ldots,\xi_{k}):=\mathbbm{1}\left(\hat{Z}_{k}(x;\xi_{1},\ldots,\xi_{k})<\min_{x^{\prime}\in\mathcal{W}}\hat{Z}_{k}(x^{\prime};\xi_{1},\ldots,\xi_{k})+\epsilon\right)\quad\text{for
}x\in\mathcal{W},\epsilon\in[0,\delta/2],$

which indicates a strict $\epsilon$-optimal solution, and let
$p_{k}^{\mathcal{W},\epsilon-},\hat{p}_{k}^{\mathcal{W},\epsilon-},\bar{p}_{k}^{\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}^{x,\mathcal{W},\epsilon}:\epsilon\in[0,\delta/2]\\}$,
i.e., let $m^{\prime}=\lfloor p_{k}^{\mathcal{W},\delta/2}(x)m\rfloor$ and

$\displaystyle\epsilon_{0}$ $\displaystyle:=0,$ $\displaystyle\epsilon_{i}$
$\displaystyle:=\inf\left\\{\epsilon\in[0,\delta/2]:p_{k}^{\mathcal{W},\epsilon}(x)\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}^{\mathcal{W},\epsilon}(x)-p_{k}^{\mathcal{W},\epsilon}(x)$
$\displaystyle\leq$
$\displaystyle\bar{p}_{k}^{\mathcal{W},\epsilon_{i+1}-}(x)-p_{k}^{\mathcal{W},\epsilon_{i}}(x)$
$\displaystyle\leq$
$\displaystyle\bar{p}_{k}^{\mathcal{W},\epsilon_{i+1}-}(x)-p_{k}^{\mathcal{W},\epsilon_{i+1}-}(x)+p_{k}^{\mathcal{W},\epsilon_{i+1}-}(x)-p_{k}^{\mathcal{W},\epsilon_{i}}(x)$
$\displaystyle\leq$
$\displaystyle\bar{p}_{k}^{\mathcal{W},\epsilon_{i+1}-}(x)-p_{k}^{\mathcal{W},\epsilon_{i+1}-}(x)+\dfrac{1}{m}$

and that

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

Therefore

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

To show that the random variable in (58) converges to $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\rfloor\cdot k$ data

$\displaystyle\tilde{p}_{k}^{\mathcal{W},\epsilon}(x)$
$\displaystyle:=\dfrac{1}{\lfloor n/k\rfloor}\sum_{i=1}^{\lfloor
n/k\rfloor}\mathbbm{1}_{k}^{x,\mathcal{W},\epsilon}(\xi_{k(i-1)+1},\ldots,\xi_{ki}),$
$\displaystyle\tilde{p}_{k}^{\mathcal{W},\epsilon-}(x)$
$\displaystyle:=\dfrac{1}{\lfloor n/k\rfloor}\sum_{i=1}^{\lfloor
n/k\rfloor}\mathbbm{1}_{k}^{x,\mathcal{W},\epsilon-}(\xi_{k(i-1)+1},\ldots,\xi_{ki}).$

Therefore we can write

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

By Minkowski inequality, the supremum satisfies

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

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

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

The bound for $\max_{x\in\mathcal{X}\backslash\mathcal{X}^{\delta}}p_{k}(x)$
from the proof of Lemma 8 continues to hold with $\mathcal{X}$ replaced by
$\mathcal{W}$, i.e.,
$\max_{x\in\mathcal{W}\backslash\mathcal{W}^{\delta}}p_{k}^{\mathcal{W},\epsilon}(x)\leq
1-q_{k}^{\epsilon,\delta,\mathcal{W}}$. Since
$q_{k}^{\epsilon,\delta,\mathcal{W}}$ decreases in $\epsilon$ we have
$\inf_{\epsilon\in[0,\delta/2]}q_{k}^{\epsilon,\delta,\mathcal{W}}=q_{k}^{\delta/2,\delta,\mathcal{W}}$,
we can write

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

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

$\Delta_{o}:=\min\left\\{Z(x^{\prime})-Z(x):x,x^{\prime}\in\mathcal{X},Z(x^{\prime})>Z(x)\right\\}>0.$

and have

$\displaystyle\left\\{\widehat{\mathcal{W}}_{k}^{\delta/2}\not\subseteq\mathcal{W}^{\delta}\right\\}$
$\displaystyle\subseteq$
$\displaystyle\bigcup_{x,x^{\prime}\in\mathcal{W}\text{ s.t.
}Z(x^{\prime})-Z(x)>\delta}\left\\{\hat{Z}_{k}(x^{\prime})\leq\hat{Z}_{k}(x)+\dfrac{\delta}{2}\right\\}$
$\displaystyle\subseteq$
$\displaystyle\bigcup_{x,x^{\prime}\in\mathcal{X}\text{ s.t.
}Z(x^{\prime})-Z(x)>\delta}\left\\{\hat{Z}_{k}(x^{\prime})-Z(x^{\prime})+Z(x^{\prime})-Z(x)\leq\hat{Z}_{k}(x)-Z(x)+\dfrac{\delta}{2}\right\\}$
$\displaystyle\subseteq$
$\displaystyle\bigcup_{x,x^{\prime}\in\mathcal{X}\text{ s.t.
}Z(x^{\prime})-Z(x)>\delta}\left\\{\hat{Z}_{k}(x^{\prime})-Z(x^{\prime})+\max(\Delta,\delta)\leq\hat{Z}_{k}(x)-Z(x)+\dfrac{\delta}{2}\right\\}$
by the definition of $\Delta_{o}$ $\displaystyle\subseteq$
$\displaystyle\bigcup_{x,x^{\prime}\in\mathcal{X}}\left\\{\hat{Z}_{k}(x^{\prime})-Z(x^{\prime})+\max\left(\Delta_{o}-\dfrac{\delta}{2},\dfrac{\delta}{2}\right)\leq\hat{Z}_{k}(x)-Z(x)\right\\}$
$\displaystyle\subseteq$
$\displaystyle\bigcup_{x,x^{\prime}\in\mathcal{X}}\left\\{\hat{Z}_{k}(x^{\prime})-Z(x^{\prime})\leq-\max\left(\dfrac{\Delta_{o}}{2}-\dfrac{\delta}{4},\dfrac{\delta}{4}\right)\text{
or
}\hat{Z}_{k}(x)-Z(x)\geq\max\left(\dfrac{\Delta_{o}}{2}-\dfrac{\delta}{4},\dfrac{\delta}{4}\right)\right\\}$
$\displaystyle\subseteq$
$\displaystyle\bigcup_{x\in\mathcal{X}}\left\\{\left|\hat{Z}_{k}(x)-Z(x)\right|\geq\max\left(\dfrac{\Delta_{o}}{2}-\dfrac{\delta}{4},\dfrac{\delta}{4}\right)\right\\}$
$\displaystyle\subseteq$
$\displaystyle\bigcup_{x\in\mathcal{X}}\left\\{\left|\hat{Z}_{k}(x)-Z(x)\right|\geq\dfrac{\Delta_{o}}{4}\right\\}$
$\displaystyle\subseteq$
$\displaystyle\left\\{\sup_{x\in\mathcal{X}}\left|\hat{Z}_{k}(x)-Z(x)\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_{x\in\mathcal{W}\backslash\mathcal{W}^{\delta}}p_{k}^{\mathcal{W},\epsilon}(x)\leq
T_{k}\left(\dfrac{\Delta_{o}}{4}\right)\to 0\text{\ \ as\ }k\to\infty.$

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

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

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

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

Due to the uniform convergence (59), we have

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

in probability, and hence

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

Finally, we combine all the pieces to

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

By Lemma 9 we have
$\mathbb{P}\left(\mathcal{S}\cap\mathcal{X}^{\delta}=\emptyset\right)\to 0$.