Estimating the Algorithmic Variance of Randomized Ensembles via the Bootstrap

Introduction

Random forests and bagging are some of the most widely used prediction methods (Breiman 1996; Breiman 2001), and over the course of the past fifteen years, much progress has been made in analyzing their statistical performance (Bühlmann and Yu 2002; Hall and Samworth 2005; Biau, Devroye and Lugosi 2008; Biau 2012; Scornet, Biau and Vert 2015). However, from a computational perspective, relatively little is understood about the algorithmic convergence of these methods, and in practice, ad hoc criteria are generally used to assess this convergence.

To clarify the idea of algorithmic convergence, recall that when bagging and random forests are used for classification, a large collection of $t$ randomized classifiers is trained, and then new predictions are made by taking the plurality vote of the classifiers. If such a method is run several times on the same training data $\mathcal{D}$, the prediction error $\textsc{Err}_{t}$ of the ensemble will vary with each run, due to the randomized training algorithm. As the ensemble size increases $(t\to\infty)$ with $\mathcal{D}$ held fixed, the random variable $\textsc{Err}_{t}$ typically decreases and eventually stabilizes at a limiting value $\text{err}_{\infty}=\text{err}_{\infty}(\mathcal{D})$. In this way, an ensemble reaches algorithmic convergence when its prediction error nearly matches that of an infinite ensemble trained on the same data.

Meanwhile, with regard to computational cost, larger ensembles are more expensive to train, to store in memory, and to evaluate on unlabeled points. For this reason, it is desirable to have a quantitative guarantee that an ensemble of a given size will perform nearly as well as an infinite one. This type of guarantee also prevents wasted computation, and assures the user that extra classifiers are unlikely to yield much improvement in accuracy.

To measure algorithmic convergence, we propose a new bootstrap method for approximating the distribution $\mathcal{L}(\sqrt{t}(\textsc{Err}_{t}-\text{err}_{\infty})|\mathcal{D})$ as $t\to\infty$. Such an approximation allows the user to decide when the algorithmic fluctuations of $\textsc{Err}_{t}$ around $\text{err}_{\infty}$ are negligible. If particular, if we refer to the algorithmic variance

as the variance of $\textsc{Err}_{t}$ due only the training algorithm, then the parameter $\sigma_{t}$ is a concrete measure of convergence that can be estimated via the bootstrap. In addition, the computational cost of the method turns out to be quite modest, by virtue of an extrapolation technique, as described in Section 4.

Although the bootstrap is an established approach to distributional approximation and variance estimation, our work applies the bootstrap in a relatively novel way. Namely, the method is based on “bootstrapping an algorithm”, rather than “bootstrapping data” — and in essence, we are applying an inferential method in order to serve a computational purpose. The opportunities for applying this perspective to other randomized algorithms can also be seen in the papers Byrd et al. 2012; Lopes, Wang and Mahoney 2017; Lopes, Wang and Mahoney 2018, which deal with stochastic gradient methods, as well as randomized versions of matrix multiplication and least-squares.

where $\operatorname{vol}(\cdot)$ is a volume measure, the symbol $\text{div}(Z)$ denotes the divergence of the vector field $Z(\theta):=\frac{\partial}{\partial\delta}f_{\delta}(\theta)\big|_{\delta=0}$, and the symbol $d\theta$ is a volume element on $\mathcal{M}$. In our analysis, it is necessary to adapt this result to a situation where the maps $f_{\delta}$ are non-smooth, the manifold $\mathcal{M}$ is a non-smooth subset of Euclidean space, and the vector field $Z(\cdot)$ is a non-smooth Gaussian process. Furthermore, applying a version of Stokes’ theorem to the right side of equation (1.1) leads to a particular linear functional of $Z(\cdot)$, which turns out to be the Hadamard derivative relevant to understanding $\textsc{Err}_{t}$. A more detailed explanation of this connection is given below equation (B.1) in Appendix B.

Among the references just mentioned, the ones that are most closely related to the current paper are Lopes 2016 and Cannings and Samworth 2017. These works derive theoretical upper bounds on $\operatorname{var}(\textsc{Err}_{t}|\mathcal{D})$ or $\operatorname{var}(\textsc{Err}_{t,l}|\mathcal{D})$, where $\textsc{Err}_{t,l}$ is the error rate on a particular class $l$ (cf. Section 4). The paper Lopes 2016 also proposes a method to estimate the unknown parameters in such bounds. In relation to these works, the current paper differs in two significant ways. First, we offer an approximation to the full distribution $\mathcal{L}(\textsc{Err}_{t}-\text{err}_{\infty}|\mathcal{D})$, and hence provide a direct estimate of algorithmic variance, rather than a bound. Second, the method proposed here is relevant to a wider range of problems, since it can handle any number of classes, whereas the analyses in Lopes 2016 and Cannings and Samworth 2017 are specialized to the binary setting. Moreover, the theoretical analysis of the bootstrap approach is entirely different from the previous techniques used in deriving variance bounds.

Outside of the setting of randomized ensemble classifiers, the papers Sexton and Laake 2009; Arlot and Genuer 2014; Wager, Hastie and Efron 2014; Mentch and Hooker 2016; Scornet 2016a look at the algorithmic fluctuations of ensemble regression functions at a fixed test point.

2 Background and setup

We consider the general setting of a classification problem with $k\geq 2$ classes. The set of training data is denoted $\mathcal{D}:=\{(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\}$, which is contained in a sample space $\mathcal{X}\times\mathcal{Y}$. The feature space $\mathcal{X}$ is arbitrary, and the space of labels $\mathcal{Y}$ has cardinality $k$. An ensemble of $t$ classifiers is denoted by $Q_{i}:\mathcal{X}\to\mathcal{Y}$, with $i=1,\dots,t$.

The key issue in studying the algorithmic convergence of bagging and random forests is randomization. In the method of bagging, randomization is introduced by generating random sets $\mathcal{D}_{1}^{*},\dots,\mathcal{D}_{t}^{*}$, each of size $n$, via sampling with replacement from $\mathcal{D}$. For each $i=1,\dots,t$, a classifier $Q_{i}$ is trained on $\mathcal{D}_{i}^{*}$, with the same classification method being used each time. When each $Q_{i}$ is trained with a decision tree method (such as CART (Breiman et al. 1984)), the random forests procedure extends bagging by adding a randomized feature selection rule (Breiman 2001).

It is helpful to note that the classifiers in bagging and random forests can be represented in a common way. Namely, there is a deterministic function, say $g$, such that for any fixed $x\in\mathcal{X}$, each classifier $Q_{i}$ can be written as

where $\xi_{1},\xi_{2},\dots$, is an i.i.d. sequence of random objects, independent of $\mathcal{D}$, that specify the “randomizing parameters” of the classifiers (cf. Breiman 2001). For instance, in the case of bagging, the object $\xi_{i}$ specifies the randomly chosen points in $\mathcal{D}_{i}^{*}$.

Beyond bagging and random forests, our proposed method will be generally applicable to ensembles that can be represented in the form (1.2), such as those in Ho 1998; Dietterich 2000; Bühlmann and Yu 2002. This representation should be viewed abstractly, and it is not necessary for the function $g$ or the objects $\xi_{1},\xi_{2},\dots$ to be explicitly constructed in practice. Some further examples include a recent ensemble method based on random projections (Cannings and Samworth 2017), as well as the voting Gibbs classifier (Ng and Jordan 2001), which is a Bayesian ensemble method based on posterior sampling. More generally, if the functions $Q_{1},Q_{2},\dots$ are i.i.d. conditionally on $\mathcal{D}$, then the ensemble can be represented in the form (1.2), as long as the classifiers lie in a standard Borel space (Kallenberg 2006, Lemma 3.22). Lastly, it is important to note that the representation (1.2) generally does not hold for classifiers generated by boosting methods (Schapire and Freund 2012), for which the analysis of algorithmic convergence is quite different.

Let $\nu=\mathcal{L}(X,Y)$ denote the distribution of a test point $(X,Y)$ in $\mathcal{X}\times\mathcal{Y}$, drawn independently of $\mathcal{D}$ and $Q_{1},\dots,Q_{t}$. Then, for a particular realization of the classifiers $Q_{1},\dots,Q_{t}$, trained with the given set $\mathcal{D}$, the prediction error rate is defined as

where $\boldsymbol{\xi}_{t}:=(\xi_{1},\dots,\xi_{t})$. (Class-wise error rates $\textsc{Err}_{t,l}$, with $l=0,\dots,k-1$ will also be addressed in Section 4.1.) Here, it is crucial to note that $\textsc{Err}_{t}$ is a random variable, since $\bar{Q}_{t}$ is a random function. Indeed, the integral above shows that $\textsc{Err}_{t}$ is a functional of $\bar{Q}_{t}$. Moreover, there are two sources of randomness to consider: the algorithmic randomness arising from $\boldsymbol{\xi}_{t}$, and the randomness arising from the training set $\mathcal{D}$. Going forward, we will focus on the algorithmic fluctuations of $\textsc{Err}_{t}$ due to $\boldsymbol{\xi}_{t}$, and our analysis will always be conditional on $\mathcal{D}$.

Recall that the value $\text{err}_{\infty}=\text{err}_{\infty}(\mathcal{D})$ represents the ideal prediction error of an infinite ensemble trained on $\mathcal{D}$. Hence, a natural way of defining algorithmic convergence is to say that it occurs when $t$ is large enough so that the condition $|\textsc{Err}_{t}-\text{err}_{\infty}|\leq\epsilon$ holds with high probability, conditionally on $\mathcal{D}$, for some user-specified tolerance $\epsilon$. However, the immediate problem we face is that it is not obvious how to check such a condition in practice.

From the right panel of Figure 1, we see that for most $t$, the inequality $|\textsc{Err}_{t}-\text{err}_{\infty}|\leq 3\sigma_{t}$ is highly likely to hold — and this observation can be formalized using Theorem 3.1 later on. For this reason, we propose to estimate $\sigma_{t}$ as a route to measuring algorithmic convergence. It is also important to note that estimating the quantiles of $\mathcal{L}(\textsc{Err}_{t}-\text{err}_{\infty}|\mathcal{D})$ would serve the same purpose, but for the sake of simplicity, we will focus on $\sigma_{t}$. In particular, there are at least two ways that an estimate $\widehat{\sigma}_{t}$ can be used in practice:

Checking convergence for a given ensemble. If an ensemble of a given size $t_{0}$ has been trained, then convergence can be checked by asking whether or not the observable condition $3\widehat{\sigma}_{t_{0}}\leq\epsilon$ holds. Additional comments on possible choices for $\epsilon$ will be given shortly.

Selecting $t$ dynamically. In order to make the training process as computationally efficient as possible, it is desirable to select the smallest $t$ needed so that $|\textsc{Err}_{t}-\text{err}_{\infty}|\leq\epsilon$ is likely to hold. It turns out that this can be accomplished using an extrapolation technique, due to the fact that $\sigma_{t}$ tends to scale like $1/\sqrt{t}$ (cf. Theorem 3.1). More specifically, if the user trains a small initial ensemble of size $t_{0}$ and computes an estimate $\widehat{\sigma}_{t_{0}}$, then “future” values of $\sigma_{t}$ for $t\gg t_{0}$ can be estimated at no additional cost with the re-scaled estimate $\sqrt{t_{0}/t}\,\widehat{\sigma}_{t_{0}}$. In other words, it is possible to look ahead and predict how many additional classifiers are needed to achieve $3\sigma_{t}\leq\epsilon$. Additional details are given in Section 4.2.

Having described the basic formulation of the problem, it is important to identify what challenges are involved in estimating $\sigma_{t}$. First, we must keep in mind that the parameter $\sigma_{t}$ describes how $\textsc{Err}_{t}$ fluctuates over repeated ensembles generated from $\mathcal{D}$ — and so it is not obvious that it is possible to estimate $\sigma_{t}$ from the output of a single ensemble. Second, the computational cost to estimate $\sigma_{t}$ should not outweigh the cost of training the ensemble, and consequently, the proposed method should be computationally efficient. These two obstacles will be described in Sections 3.2 and 4.2 respectively.

Our proposed bootstrap method is described in Section 2, and our main consistency result is given in Section 3. Practical considerations are discussed in Section 4, numerical experiments are given in Section 5, and conclusions are stated in Section 6. The essential ideas of the proofs are explained in Appendices A and B, while the technical arguments are given Appendices C-E. Lastly, in Appendix F, we provide additional assessment of technical assumptions. All appendices are in the supplementary material.

Method

Based on the definition of $\textsc{Err}_{t}$ in equation (1.3), we may view $\textsc{Err}_{t}$ as a functional of $\bar{Q}_{t}$, denoted

From a statistical standpoint, the importance of this expression is that $\textsc{Err}_{t}$ is a functional of a sample mean, which makes it plausible that $\sigma_{t}$ is amenable to bootstrapping, provided that $\varphi$ is sufficiently smooth.

To describe the bootstrap method, let $(Q_{1}^{*},\dots,Q_{t}^{*})$ denote a random sample with replacement from the trained ensemble $(Q_{1},\dots,Q_{t})$, and put $\bar{Q}_{t}^{*}(\cdot):=\textstyle\frac{1}{t}\textstyle\sum_{i=1}^{t}Q_{i}^{*}(\cdot).$ In turn, it would be natural to regard the quantity

as a bootstrap sample of $\textsc{Err}_{t}$, but strictly speaking, this is an “idealized” bootstrap sample, because the functional $\varphi$ depends on the unknown test point distribution $\nu=\mathcal{L}(X,Y)$. Likewise, in Section 2.1 below, we explain how each value $\varphi(\bar{Q}_{t}^{*})$ can be estimated. So, in other words, if $\widehat{\varphi}$ denotes an estimate of $\varphi$, then an estimate of $\textsc{Err}_{t}$ would be written as

and the corresponding bootstrap sample is

Altogether, a basic version of the proposed bootstrap algorithm is summarized as follows.

Sample $t$ classifiers $(Q_{1}^{*},\dots,Q_{t}^{*})$ with replacement from $(Q_{1},\dots,Q_{t})$.

Compute $z_{b}:=\widehat{\varphi}(\bar{Q}_{t}^{*})$.

Return: the sample standard deviation of $z_{1},\dots,z_{B}$, denoted $\widehat{\sigma}_{t}$.

While the above algorithm is conceptually simple, it suppresses most of the implementation details, and these are explained below. Also note that in order to approximate quantiles of $\mathcal{L}(\textsc{Err}_{t}-\text{err}_{\infty}|\mathcal{D})$, rather than $\sigma_{t}$, it is only necessary to modify the last step, by returning the desired quantile of the centered values $z_{1}-\bar{z},\dots,z_{B}-\bar{z}$, with $\bar{z}=\frac{1}{B}\sum_{b=1}^{B}z_{b}$.

1 Resampling algorithm with hold-out or “out-of-bag” points

Return: the sample standard deviation of $z_{1},\dots,z_{B}$, denoted $\widehat{\sigma}_{t}$.

Since the use of a hold-out set is often undesirable in practice, we instead consider oob points — which are a special feature of bagging and random forests. To briefly review this notion, recall that each classifier $Q_{i}$ is trained on a set of $n$ points $\mathcal{D}_{i}^{*}$ obtained by sampling with replacement from $\mathcal{D}$. Consequently, each set $\mathcal{D}_{i}^{*}$ excludes approximately $(1-\textstyle\frac{1}{n})^{n}\approx$ 37% of the points in $\mathcal{D}$, and these excluded points may be used as test points for the particular classifier $Q_{i}$. If a point $X_{j}$ is excluded from $\mathcal{D}_{i}^{*}$, then we say “the point $X_{j}$ is oob for the classifer $Q_{i}$”, and we write $i\in\textsc{{oob}}(X_{j})$, where the set $\textsc{{oob}}(X_{j})\subset\{1,\dots,t\}$ indexes the classifiers for which $X_{j}$ is oob.

Main result

Our main theoretical goal is to prove that the bootstrap yields a consistent approximation of $\mathcal{L}(\sqrt{t}(\textsc{Err}_{t}-\text{err}_{\infty})|\mathcal{D})$ as $t$ becomes large. Toward this goal, we will rely on two simplifications that are customary in analyses of bootstrap and ensemble methods. First, we will exclude the Monte-Carlo error arising from the finite number of $B$ bootstrap replicates, as well as the error arising from the estimation of $\textsc{Err}_{t}$. For this reason, our results do not formally require the training or hold-out points to be i.i.d. copies of the test point $(X,Y)$ — but from a practical standpoint, it is natural to expect that this type of condition should hold in order for Algorithm 2 (or its oob version) to work well.

Second, we will analyze a simplified type of ensemble, which we will refer to as a first-order model. This type of approach has been useful in gaining theoretical insights into the behavior of complex ensemble methods in a variety of previous works Biau, Devroye and Lugosi 2008; Biau 2012; Arlot and Genuer 2014; Lin and Jeon 2006; Genuer 2012; Scornet 2016a; Scornet 2016b. In our context, the value of this simplification is that it neatly packages the complexity of the base classifiers, and clarifies the relationship between $t$ and quality of the bootstrap approximation. Also, even with such simplifications, the theoretical problem of proving bootstrap consistency still leads to considerable technical challenges. Lastly, it is important to clarify that the first-order model is introduced only for theoretical analysis, and our proposed method does not rely on this model.

Any randomized classifier $Q_{1}:\mathcal{X}\to\{\boldsymbol{e}_{0},\dots,\boldsymbol{e}_{k-1}\}$ may be viewed as a stochastic process indexed by $\mathcal{X}$. From this viewpoint, we say that another randomized classifier $T_{1}:\mathcal{X}\to\{\boldsymbol{e}_{0},\dots,\boldsymbol{e}_{k-1}\}$ is a first-order model for $Q_{1}$ if it has the same marginal distributions as $Q_{1}$, conditionally on $\mathcal{D}$, which means

Since $Q_{1}(x)$ takes values in the finite set of binary vectors $\{\boldsymbol{e}_{0},\dots,\boldsymbol{e}_{k-1}\}$, the condition (3.1) is equivalent to

where the expectation is only over the algorithmic randomness in $Q_{1}$ and $T_{1}$. A notable consequence of this matching condition is that the ensembles associated with $Q_{1}$ and $T_{1}$ have the same error rates on average. Indeed, if we let $\textsc{Err}_{t}^{\prime}$ be the error rate associated with an ensemble of $t$ independent copies of $T_{1}$, then it turns out that

for all $t\geq 1$, where $\textsc{Err}_{t}$ is the error rate for $Q_{1},\dots,Q_{t}$, as before. (A short proof is given in Appendix E.) In this sense, a first-order model $T_{1}$ is a meaningful proxy for $Q_{1}$ with regard to statistical performance — even though the internal mechanisms of $T_{1}$ may be simpler.

Having stated some basic properties that are satisfied by any first-order model, we now construct a particular version that is amenable to analysis. Interestingly, it is possible to start with an arbitrary random classifier $Q_{1}:\mathcal{X}\to\{\boldsymbol{e}_{0},\dots,\boldsymbol{e}_{k-1}\}$, and construct an associated $T_{1}$ in a relatively explicit way.

To do this, let $x\in\mathcal{X}$ be fixed, and consider the function

For any fixed $\theta\in\Delta$, there is an associated partition of the unit interval into sub-intervals

such that the width of interval $I_{l}(\theta)$ is equal to $\theta_{l}$ for $l\geq 1$. Namely, we put $I_{1}(\theta):=[0,\theta_{1}]$, and for $l=2,\dots,k-1$,

Now, if we let $x\in\mathcal{X}$ be fixed, and let $U_{1}\sim$ Uniform$$, then we define$T_{1}(x)\in\{\boldsymbol{e}_{0},\dots,\boldsymbol{e}_{k-1}\}$to have its$l$th coordinate equal to the following indicator variable

where $l=1,\dots,k-1$. It is simple to check that the first-order matching condition (3.2) holds, and so $T_{1}$ is indeed a first-order model of $Q_{1}$. Furthermore, given that $T_{1}$ is defined in terms of a single random variable $U_{1}\sim$ Uniform$$, we obtain a corresponding “first-order ensemble”$T_{1},\dots,T_{t}$via an i.i.d. sample of uniform variables$U_{1},\dots,U_{t}$, which are independent of$\mathcal{D}$. (The$l$th coordinate of the$i$th classifier$T_{i}$is given by$[T_{i}(x)]_{l}=1\{U_{i}\in I_{l}(\vartheta(x))\}$.) Hence, with regard to the representation$Q_{i}(x)=g(x,\mathcal{D},\xi_{i})$ in equation (1.2), we may make the identification

when the first-order model holds with $Q_{i}=T_{i}$.

To understand the statistical meaning of the first-order model, it is instructive to consider the simplest case of binary classification, $k=2$. In this case, $T_{1}(x)$ is a Bernoulli random variable, where $T_{1}(x)=1\{U_{1}\leq\vartheta(x)\}$. Since $\bar{Q}_{t}(x)\to\vartheta(x)$ almost surely as $t\to\infty$ (conditionally on $\mathcal{D}$), the majority vote of an infinite ensemble has a similar form, i.e. $1\{\textstyle\frac{1}{2}\leq\vartheta(x)\}$. Hence, the classifiers $\{T_{i}\}$ can be viewed as “random perturbations” of the asymptotic majority vote arising from $\{Q_{i}\}$. Furthermore, if we view the number $\vartheta(x)$ as score to be compared with a threshold, then the variable $U_{i}$ plays the role of a random threshold whose expected value is $\frac{1}{2}$. Lastly, even though the formula $T_{1}(x)=1\{U_{1}\leq\vartheta(x)\}$ might seem to yield a simplistic classifier, the complexity of $T_{1}$ is actually wrapped up in the function $\vartheta$. Indeed, the matching condition (3.2) allows for the function $\vartheta$ to be arbitrary.

2 Bootstrap consistency

We now state our main result, which asserts that the bootstrap “works” under the first-order model. To give meaning to bootstrap consistency, we first review the notion of conditional weak convergence.

If a test point $X$ is drawn from class $Y=\boldsymbol{e}_{l}$, then we denote the distribution of the random vector $\vartheta(X)$, conditionally on $\mathcal{D}$, as

For the given set $\mathcal{D}$, and each $l=0,\dots,k-1$, the distribution $\mu_{l}$ has a density $f_{l}:\Delta\to[0,\infty)$ with respect to Lebesgue measure on $\Delta$, and $f_{l}$ is continuous on $\Delta$. Also, if $\Delta^{\circ}$ denotes the interior of $\Delta$, then for each $l$, the density $f_{l}$ is $C^{1}$ on $\Delta^{\circ}$, and $\|\nabla f_{l}\|_{2}$ is bounded on $\Delta^{\circ}$.

Suppose that the first-order model $Q_{i}=T_{i}$ holds for all $i\geq 1$, and that Assumption 1 holds. Then, for the given set $\mathcal{D}$, there are numbers $\textup{err}_{\infty}=\textup{err}_{\infty}(\mathcal{D})$ and $\sigma=\sigma(\mathcal{D})$ such that as $t\to\infty$,

In a nutshell, the proof of Theorem 3.1 is composed of three pieces: showing that $\textsc{Err}_{t}$ can be represented as a functional of an empirical process (Appendix A.1), establishing the smoothness of this functional (Appendix A.2), and employing the functional delta method (Appendix A.3). With regard to theoretical techniques, there are two novel aspects of the proof. The problem of deriving this functional is solved by introducing a certain lifting operator, while the problem of showing smoothness is based on a non-smooth instance of the first-variation formula, as well as some special properties of Bernstein polynomials. Lastly, it is worth mentioning that the core technical result of the paper is Theorem A.1.

Practical considerations

In this section, we discuss some considerations that arise when the proposed method is used in practice, such as the choice of error rate, the computational cost, and the choice of a stopping criterion for algorithmic convergence.

In some applications, class-wise error rates may be of greater interest than the total error rate $\textsc{Err}_{t}$. For any $l=0,\dots,k-1$, let $\nu_{l}=\mathcal{L}(X|Y=\boldsymbol{e}_{l})$ denote the distribution of the test point $X$ given that it is drawn from class $l$. Then, the error rate on class $l$ is defined as

and the corresponding algorithmic variance is

In order to estimate $\sigma_{t,l}$, Algorithm 2 can be easily adapted using either hold-out or oob points from a particular class. Our theoretical analysis also extends immediately to the estimation of $\sigma_{t,l}$ (cf. Section A.1).

2 Computational cost and extrapolation

To explain the technique of extrapolation, the first step produces an inexpensive estimate $\widehat{\sigma}_{t_{0}}$ by applying Algorithm 2 to a small initial ensemble of size $t_{0}$. The second step then rescales $\widehat{\sigma}_{t_{0}}$ so that it approximates $\sigma_{t}$ for $t\gg t_{0}$. This rescaling relies on Theorem 3.1, which leads to the approximation, $\sigma_{t}\approx\textstyle\frac{\sigma}{\sqrt{t}}$. Consequently, we define the extrapolated estimate of $\sigma_{t}$ as

In turn, if the user desires $3\sigma_{t}\leq\epsilon$ for some $\epsilon\in(0,1)$, then $t$ should be chosen so that

which is equivalent to $t\geq(\textstyle\frac{3\sqrt{t_{0}}}{\epsilon}\cdot\widehat{\sigma}_{t_{0}})^{2}$.

In addition to applying Algorithm 2 to a smaller ensemble, a second computational benefit is that extrapolation allows the user to “look ahead” and dynamically determine how much extra computation is needed so that $\sigma_{t}$ is within a desired range. In Section 5, some examples are given showing that $\sigma_{1,000}$ can be estimated well via extrapolation when $t_{0}=200$.

Based on the reasoning just given, the cost of running Algorithm 2 does not exceed the cost of training $t$ trees, provided that

where the factor $\frac{t}{t_{0}}$ arises from the extrapolation speedup described earlier. Moreover, with regard to the selection of $B$, our numerical examples in Section 5 show that the modest choice $B=50$ allows Algorithm 2 to perform well on a variety of datasets.

Numerical Experiments

To illustrate our proposed method, we describe experiments in which the random forests method is applied to natural and synthetic datasets (6 in total). More specifically, we consider the task of estimating the parameter 3$\sigma_{t}=3\sqrt{\operatorname{var}(\textsc{Err}_{t}|\mathcal{D})}$, as well as $3\sigma_{t,l}=3\sqrt{\operatorname{var}(\textsc{Err}_{t,l}|\mathcal{D})}$. Overall, the main purpose of the experiments is to show that the bootstrap can indeed produce accurate estimates of these parameters. A second purpose is to demonstrate the value of the extrapolation technique from Section 4.2.

Each of the 6 datasets were partitioned in the following way. First, each dataset was evenly split into a training set $\mathcal{D}$ and a “ground truth” set $\mathcal{D}_{\text{ground}}$, with nearly matching class proportions in $\mathcal{D}$ and $\mathcal{D}_{\text{ground}}$. (The reason that a substantial portion of data was set aside for $\mathcal{D}_{\text{ground}}$ was to ensure that ground truth values of $\sigma_{t}$ and $\sigma_{t,l}$ could be approximated using this set.) Next, a smaller set $\mathcal{D}_{\text{hold}}\subset\mathcal{D}_{\text{ground}}$ with cardinality satisfying $|\mathcal{D}_{\text{hold}}|/(|\mathcal{D}_{\text{hold}}|+|\mathcal{D}|)\approx 1/6$ was used as the hold-out set for implementing Algorithm 2. As before, the class proportions in $\mathcal{D}_{\text{hold}}$ and $\mathcal{D}$ were nearly matching. The smaller size of $\mathcal{D}_{\text{hold}}$ was chosen to illustrate the performance of the method when hold-out points are limited.

After preparing $\mathcal{D}$, $\mathcal{D}_{\text{ground}}$, and $\mathcal{D}_{\text{hold}}$, a collection of 1,000 ensembles was trained on $\mathcal{D}$ by repeatedly running the random forests method. Each ensemble contained a total of 1,000 trees, trained under default settings from the package randomForest (Liaw and Wiener 2002). Also, we tested each ensemble on $\mathcal{D}_{\text{ground}}$ to approximate a corresponding sample path of $\textsc{Err}_{t}$ (like the ones shown in Figure 1). Next, in order to obtain “ground truth” values for $\sigma_{t}$ with $t=1,\dots,1,\!000$, we used the sample standard deviation of the 1,000 sample paths at each $t$. (Ground truth values for each $\sigma_{t,l}$ were obtained analogously.)

With regard to our methodology, we applied the hold-out and oob versions of Algorithm 2 to each of the ensembles — yielding 1,000 realizations of each type of estimate of $\sigma_{t}$. In each case, the number of bootstrap replicates was set to $B=50$, and we applied the extrapolation rule, starting from $t_{0}=200$. If we let $\widehat{\sigma}_{200,\textsc{h}}$ and $\widehat{\sigma}_{200,\textsc{o}}$ denote the initial hold-out and oob estimators, then the corresponding extrapolated estimators for $t\geq 200$ are given by

Next, as a benchmark, we considered an enhanced version of the hold-out estimator, for which the entire ground truth set $\mathcal{D}_{\text{ground}}$ was used in place of $\mathcal{D}_{\text{hold}}$. In other words, this benchmark reflects a situation where a much larger hold-out set is available, and it is referred to as the “ground estimate” in the plots. Its value based on $t_{0}=200$ is denoted $\widehat{\sigma}_{200,\textsc{g}}$, and for $t\geq 200$, we use

to refer to its extrapolated version. Lastly, class-wise versions of all extrapolated estimators were computed in an analogous way.

2 Description of datasets

The following datasets were each partitioned into $\mathcal{D}$, $\mathcal{D}_{\text{hold}}$ and $\mathcal{D}_{\text{ground}}$, as described above.

A set of census records for 48,842 people were collected with 14 socioeconomic features (continuous and discrete) (Lichman 2013). Each record was labeled as 0 or 1, corresponding to low or high income. The proportions of the classes are approximately (.76,.24). As a pre-processing step, we excluded three features corresponding to work-class, occupation, and native country, due to a high proportion of missing values.

The observations represent 67,557 board positions in the two-person game “connect-4” (Lichman 2013). For each position, a list of 42 categorical features are available, and each position is labeled as a draw $l=0$, loss $l=1$, or win $l=2$ for the first player, with the class proportions being approximately $(.10,.25,.65)$.

This dataset was prepared from a set of 12,960 applications for admission to a nursery school (Lichman 2013). Each application was associated with a list of 8 (categorical) socioeconomic features. Originally, each application was labeled as one of five classes, but in order to achieve reasonable label balance, the last three categories were combined. This led to approximate class proportions $(1/3,1/3,1/3)$.

A collection of 39,797 news articles from the website mashable.com were associated with 60 features (continuous and discrete). Each article was labeled based on the number of times it was shared: fewer than 1000 shares $(l=0)$, between 1,000 and 5,000 shares ($l=1$), and greater than 5,000 shares ($l=2$), with approximate class proportions $(.28,.59,.13)$.

3 Numerical results

For each dataset, we plot the ground truth value $3\sigma_{t}$ as a function of $t=1,\dots,1,\!000$, where the y-axis is expressed in units of %, so that a value $3\sigma_{t}=.01$ is marked as 1%. Alongside each curve for $3\sigma_{t}$, we plot the averages of $3\widehat{\sigma}_{t,\textsc{o},\text{extrap}}$ (green) , $3\widehat{\sigma}_{t,\textsc{h},\text{extrap}}$ (purple), and $3\widehat{\sigma}_{t,\textsc{g},\text{extrap}}$ (orange) over their 1,000 realizations, with error bars indicating the spread between the 10th and 90th percentiles of the estimates. Here, the error bars are only given to illustrate the variance of the estimates, conditionally on $\mathcal{D}$, and they are not proposed as confidence intervals for $\sigma_{t}$. (Indeed, our main focus is on the fluctuations of $\textsc{Err}_{t}$, rather than the fluctuations of variance estimates.) Lastly, we plot results for the class-wise parameters $\sigma_{t,l}$ in the same manner, but in order to keep the number of plots manageable, we only display the class $l$ with the highest value of $3\sigma_{t,l}$ at $t=1,\!000$. This is reflected in the plots, since the values of $3\sigma_{t,l}$ for the chosen class $l$ are generally larger than $3\sigma_{t}$.

To explain the plots from the user’s perspective, suppose the user trains an initial ensemble of $t_{0}=200$ classifiers with the ‘census income’ data. (The following considerations will apply in the same way to the other datasets in Figures 3-7.) At this stage, the user may compute either of the estimators $\widehat{\sigma}_{200,\textsc{o}}$ or $\widehat{\sigma}_{200,\textsc{h}}$. In turn, the user may follow the definitions (5.1) to plot the extrapolated estimators for all $t\geq 200$ at no additional cost. These curves will look like the purple or green curves in the left panel of Figure 2, up to a small amount of variation indicated by the error bars.

If the user wants to select $t$ so that $3\sigma_{t}$ is at most, say 0.5%, then the purple or green curves in the left panel of Figure 2 would tell the user that 200 classifiers are already sufficient, and no extra classifiers are needed (which is correct in this particular example). Alternatively, if the user happens to be interested in the class-wise error rate for $l=1$, and if the user wants $3\sigma_{t,1}$ to be at most 0.5%, then the curve for the oob estimator accurately predicts that approximately 600 total (i.e. 400 extra) classifiers are needed. By contrast, the hold-out method is conservative, and indicates that approximately 1,000 total (i.e. 800 extra) classifiers should be trained. So, in other words, the hold-out estimator would still provide the user with the desired outcome, but at a higher computational cost.

Considering all of the datasets collectively, the plots show that the extrapolated oob and ground estimators are generally quite accurate. Meanwhile, the hold-out estimator tends to be conservative, due to an upward bias. Consequently, the oob method should be viewed as preferable, since it is both more accurate, and does not require data to be held out. Nevertheless, when considering the hold-out estimator, it is worth noting that the effect of the bias actually diminishes with extrapolation, and even if the initial value $3\widehat{\sigma}_{t_{0},\textsc{h},\text{extrap}}$ has noticeable bias at $t_{0}=200$, it is possible for the extrapolated value $3\widehat{\sigma}_{t,\textsc{h},\text{extrap}}$ to have relatively small bias at $t=1,\!000$.

One last point to mention is that many of the datasets have discrete features, which may violate the theoretical conditions in Assumption 1. Nevertheless, the presence or absence of discrete features does not seem to substantially affect on the performance of the estimators. So, to this extent, the bootstrap does not seem to depend too heavily on Assumption 1. (See Appendix F.2 for further empirical assessment of that assumption.)

Conclusion

We have studied the notion of algorithmic variance $\sigma_{t}^{2}=\operatorname{var}(\textsc{Err}_{t}|\mathcal{D})$ as a criterion for deciding when a randomized ensemble will perform nearly as well as an infinite one (trained on the same data). To estimate this parameter, we have developed a new bootstrap method, which allows the user to directly measure the convergence of randomized ensembles with a guarantee that has not previously been available.

With regard to practical considerations, we have shown that our bootstrap method can be enhanced in two ways. First, the use of a hold-out set can be avoided with the oob version of Algorithm 2, and our numerical results show that the oob version is preferable when hold-out points are scarce. Second, the extrapolation technique substantially reduces the cost of bootstrapping. Furthermore, we have analyzed the cost of the method in terms of floating point operations to show that it compares favorably with the cost of training a single ensemble via random forests.

Acknowledgements

The author thanks Peter Bickel, Philip Kegelmeyer, and Debashis Paul for helpful discussions. In addition, the author thanks the editors and referees for their valuable feedback, which significantly improved the paper.

References

Outline of appendices and the proof of Theorem 3.1

Here we explain how the the main components of the proof of Theorem 3.1 fit together. First, in Appendix A.1, we show that under a first-order model, there is an explicit functional $\phi$ such that

where $U_{1}^{*},\dots,U_{t}^{*}$ are drawn with replacement from $\textup{\bf{U}}_{t}$, then the bootstrap counterpart of (1) is given by

The remainder of the appendices are organized as follows. Appendices B, C, and D contain the arguments for proving Theorem A.1 on Hadamard differentiability. Throughout these arguments, we will refer to various technical lemmas that are stated and proved in Appendix E. Also, we henceforth assume that the first-order model holds, so that $Q_{i}=T_{i}$ for $i\geq 1$, and the sets $\textup{\bf{U}}_{t}$ and $\boldsymbol{\xi}_{t}$ are synonymous. Lastly, Appendix F discusses the theoretical and empirical assessment of Assumption 1.

Appendix A Proof of Theorem 3.1

Working under the first-order model, our aim in this subsection is to construct an explicit functional $\phi$ such that

which implies that $\phi$ may be written as $\phi=\textstyle\sum_{l=0}^{k-1}\pi_{l}\phi_{l}$.

Since many of our arguments will rely on special properties of $\boldsymbol{L}$, we briefly summarize these properties below.

The fact that $\boldsymbol{L}$ respects composition of functions is the only property that takes some care to verify, but we omit the calculation for brevity. In turn, the “Inverses” property follows by combining the “Composition” and “Identities” properties.

where the set $\mathcal{S}_{l}\subset\Delta$ is defined as

where $\theta_{0}:=1-(\theta_{1}+\cdots+\theta_{k-1})$. Consequently, we have

A.2 Hadamard differentiability of ϕl\phi_{l}

Let B be a normed space. A map $\psi:\mathcal{F}\to\textup{\bf{B}}$ is Hadamard differentiable at $G_{0}\in\mathcal{F}$ tangentially to $C$, if there is a continuous linear map $\psi_{G_{0}}^{\prime}:C\to\textup{\bf{B}}$ such that as $t\to\infty$,

for all converging sequences of positive numbers $\varepsilon_{t}\to 0$ and functions $h_{t}\to h$, such that $G_{0}+\varepsilon_{t}h_{t}\in\mathcal{F}$ for every $t\geq 1$, and $h\in C$. In particular, the linear map $\psi^{\prime}_{G_{0}}$ is referred to as the Hadamard derivative of $\psi$ at $G_{0}$.

Although each functional $\phi_{l}$ is defined in terms of the particular set $\mathcal{S}_{l}$ and the particular measure $\mu_{l}$, the Hadamard differentiability of $\phi_{l}$ is only mildly dependent on their structure. For each measure $\mu_{l}$, the only property we need is that it satisfies Assumption 1. Regarding the set $\mathcal{S}_{l}$, it is simple to check that its complement in $\Delta$ is a convex set with non-empty interior. In other words, the functional $1-\phi_{l}(G)$ may be written as $\mu_{l}\big([\boldsymbol{L}(G)]^{-1}(\mathcal{S})\big)$ for some convex set $\mathcal{S}\subset\Delta$ with non-empty interior. So, given that the Hadamard derivative of $1-\phi_{l}$ is that same as that of $\phi_{l}$, up to a sign, we state the result in terms of a generic functional $\psi$ that arises from such a set $\mathcal{S}$, and such a measure $\mu$.

where $\mathbf{n}(\theta)$ is the outward normal to $\partial\mathcal{S}$ at the point $\theta$, and $d\sigma$ is Hausdorff measure on $\partial\mathcal{S}$.

A high-level proof is given in Appendix B. In the next subsection, we apply this result in conjunction with the functional delta method to complete the proof of bootstrap consistency (Theorem 3.1).

A.3 Functional delta method

With this lemma in hand, Theorem 3.1 on bootstrap consistency follows quickly from Theorem A.1. Specifically, if we consider $\phi=\sum_{l=0}^{k-1}\pi_{l}\phi_{l}$, with each $\phi_{l}$ as in equation (A.7), and define $\text{err}_{\infty}:=\phi(\text{id})$, then the relations (1) and (2) lead to

where we recall that $\boldsymbol{\xi}_{t}=\textup{\bf{U}}_{t}$ in the first-order model, as explained on p.3.1.1 of the main text. Note also that $\phi$ implicitly depends on $\mathcal{D}$ through the function $\vartheta$.

Appendix B A high-level proof of Theorem A.1

Here we give a proof of Theorem A.1 that focuses on the key ideas and delegates the technical pieces to Appendices C, D, and E. Consider a sequence of positive numbers $\varepsilon_{t}\to 0$ and functions $h_{t}\to h$ such that $\text{id}_{}+\varepsilon_{t}h_{t}\in\mathcal{F}$ for every $t\geq 1$, and $h\in C$. Define the distribution function $F_{t}:\to$ by

and define its lifted version $V_{t}:\Delta\to\Delta$ by

The fact that the range of $V_{t}$ is contained in $\Delta$ follows from the properties of $\boldsymbol{L}$ listed earlier. Our aim is to evaluate the limit of the following difference as $t\to\infty$,

Here, we have used the fact that $\psi(\text{id}_{})=\mu(\mathcal{S})$, which follows from $\boldsymbol{L}(\text{id}_{})=\text{id}_{\Delta}$. Since $V_{t}$ approaches $\text{id}_{\Delta}$ as $t\to\infty$, we may view the preimage $V_{t}^{-1}(\mathcal{S})$ as a perturbed version of the set $\mathcal{S}$. From this perspective, it is natural to interpret the right side of equation (B.1) through the lens of the first variation formula, introduced in Section 1.1, with $\mu$ playing the role of a volume, $\mathcal{S}$ playing the role of a manifold, and $V_{t}$ playing the role of the map $f_{\delta}$ with $\delta=\varepsilon_{t}$.

In its classical form, the first variation formula deals with smooth maps on smooth manifolds. However, since the map $V_{t}$ need not be smooth, our proof proceeds by constructing a smoothed version of $V_{t}$. In order to do this, it is enough to smooth the univariate function $F_{t}$ and apply the linear operator $\boldsymbol{L}$. The smoothing will be done using the linear Bernstein smoothing operator, denoted $\mathcal{B}_{s}$, where $s\geq 1$ is an integer-valued smoothing parameter (lorentz; devoreconstructive).

where $b_{j}(u;s):=\binom{s}{j}u^{j}(1-u)^{s-j}$ is the $j$th Bernstein basis polynomial where $u\in$, and $j$ ranges over $\{0,1,\dots,s\}$. Below, we will use of some special properties of this operator. First, when $G$ is a cumulative distribution function, it turns out that $\mathcal{B}_{s}(G)$ is also a cumulative distribution function. Second, when $G$ is continuous, we have the uniform limit $\mathcal{B}_{s}(G)\to G$ as $s\to\infty$. The details of all the properties of $\mathcal{B}_{s}$ we will use are summarized in Lemma E.1 of Appendix E.

When applying the operator $\mathcal{B}_{s}$, the smoothed version of $F_{t}$ will be denoted by

Likewise, the smoothed version of $V_{t}$ will be denoted by

which is a map from $\Delta$ to itself. (It is not immediately obvious that $V_{t,s}$ takes values in $\Delta$, and this follows from $F_{t,s}$ being a cumulative distribution function, by Lemma E.1, as well as the “Simplices” property of $\boldsymbol{L}$.)

The remainder of the proof involves two essential parts. First, we prove a special version of the first variation formula for the smoothed maps $V_{t,s}$. Second, we show that this smoothing leads to negligible approximation error. To quantify the approximation error from smoothing, define the remainder $R_{t,s}$ according to the following equation

Due to the smoothness of $V_{t,s}$, the difference quotient on the right may be represented with a change of variable formula

where $J(V_{t,s}^{-1})(\theta)$ is the Jacobian matrix of $V_{t,s}^{-1}$ at the point $\theta$. This step is justified by Lemmas E.4, E.6, and E.7, which also prove invertibility of $V_{t,s}$. From the above integral formula, Proposition C.1 in Appendix C provides the following limit

Next, Proposition D.1 in Appendix D shows that replacing $V_{t}^{-1}$ with its smoothed version $V_{t,s}^{-1}$ leads to negligible approximation error, i.e.

Consequently, by separately applying the operations $\liminf_{t\to\infty}$ and $\limsup_{t\to\infty}$ to equation (B.3), and then taking $\lim_{s\to\infty}$ in each case, it follows that

Lastly, it is simple to check that the right side is a continuous linear functional of $h$, as required by the definition of Hadamard differentiability. ∎

Appendix C A first variation formula

Assume the conditions of Theorem A.1. Then, in the notation of Appendix B, the following limit holds,

Using an expansion for $|\det J(V_{t,s}^{-1})(\theta)|$ given in Lemma E.5 of Appendix E, as well as the boundedness of $f$ on $\Delta$, the integral on the left side of equation (C.1) may be written as

where $\text{div}\,W_{t,s}(V_{t,s}^{-1}(\theta))$ is the divergence of $W_{t,s}$ evaluated at $V_{t,s}^{-1}(\theta)$, and $K_{s}\in[0,\infty)$ is a sequence of numbers not depending on $t$. (In addition, note that $V_{t,s}^{-1}(\theta)$ lies in $\Delta^{\circ}$ whenever $\theta\in\Delta^{\circ}$, which follows from Lemma E.4, and the invariance of domain principle (Hatcher, Theorem 2B.3).)

We now evaluate the limits of the two integrals in line (C.4) separately. Due to the multivariate mean value theorem, for each fixed $\theta\in\Delta^{\circ}$, there is a point $\zeta_{t,s}(\theta)\in\Delta^{\circ}$ on the line segment between $V_{t,s}^{-1}(\theta)$ and $\theta$ such that

Also, the points $\zeta_{t,s}(\theta)$ may be taken to satisfy $\zeta_{t,s}(\theta)\to\theta$ as $t\to\infty$, since $\lim_{t\to\infty}V_{t,s}^{-1}(\theta)=\theta$ for every fixed $\theta\in\Delta^{\circ}$ and fixed $s\geq 1$, which follows from Lemma E.9. By the Cauchy-Schwarz inequality, the inner product in equation (C.5) is dominated by a number depending only on $s$, since $\|\nabla f\|_{2}$ is bounded on $\Delta^{\circ}$ by assumption, and $\sup_{t\geq 1}\sup_{\theta\in\Delta^{\circ}}\|\widetilde{W}_{t,s}(\theta)\|_{2}$ is finite by Lemma E.9. Furthermore, Lemma E.9 gives the convergence $\lim_{t\to\infty}\widetilde{W}_{t,s}(\theta)\to-W_{s}(\theta)$ for every $\theta\in\Delta^{\circ}$, where we define

Hence, the continuity of $\nabla f$ and the dominated convergence theorem lead to

Turning our attention to the second integral in line (C.4), Lemma E.3 ensures that the divergence

converges uniformly to $\text{div}\,W_{s}(\theta)$ on $\Delta^{\circ}$ as $t\to\infty$. Combining this with the facts that $\lim_{t\to\infty}V_{t,s}^{-1}(\theta)=\theta$ for every $\theta\in\Delta^{\circ}$ (Lemma E.9) and that $\text{div}(W_{s})$ is continuous (since $W_{s}$ is smooth), it follows that

Furthermore, it is simple to check that this pointwise limit is dominated by a constant (using Lemma E.3), and so

The preceding calculations may now be combined using the basic differential identity

which shows that the quantity in line (C.4) tends to $-\int_{\mathcal{S}^{\circ}}\text{div}\big(f\,W_{s}\big)(\theta)d\theta$ as $t\to\infty$ with $s$ held fixed. In turn, Stokes’ theorem may be applied to this divergence integral. To be specific, an applicable version of the Stokes’ theorem that holds for convex domains may be found in Leoni. (When referencing that result, note that bounded convex domains have Lipschitz boundaries (Grisvard, Corollary 1.2.2.3), and also, that the regularity assumptions on $f$ imply that $f\,W_{s}$ is a Lipschitz vector field on $\Delta$.) Hence, for every $s\geq 1$,

Finally, since $h\in C$, the uniform approximation property of Bernstein polynomials for continuous functions (Lemma E.1) implies that

uniformly on $\Delta$ as $s\to\infty$. Hence, the right side of equation (C.7) tends to $\int_{\partial\mathcal{S}}\langle-\boldsymbol{L}(h)(\theta),\mathbf{n}(\theta)\rangle f(\theta)d\sigma(\theta)$ by the dominated convergence theorem, which completes the proof. ∎

Appendix D Smoothing error is negligible

Let $R_{t,s}$ be as defined in equation (B.3), and suppose the conditions of Theorem A.1 hold. Then,

It is simple to check that $R_{t,s}$ may be written as

We will argue that both terms on the right side are small. Note that the set $V_{t}^{-1}(\mathcal{S})\setminus V_{t,s}^{-1}(\mathcal{S})$ consists of the points $\theta\in\Delta$ such that $V_{t}(\theta)\in\mathcal{S}$ and $V_{t,s}(\theta)\not\in\mathcal{S}$. Now consider a particular point $\theta\in V_{t}^{-1}(\mathcal{S})\setminus V_{t,s}^{-1}(\mathcal{S})$, and note that if the Euclidean distance between the points $V_{t}(\theta)$ and $V_{t,s}(\theta)$ is written as $d_{t,s}(\theta)$, then both of the points $V_{t,s}(\theta)$ and $V_{t}(\theta)$ must be within a distance $d_{t,s}(\theta)$ from the boundary $\partial\mathcal{S}$. In other words, both of the points $V_{t}(\theta)$ and $V_{t,s}(\theta)$ lie within the tubular neighborhood of $\partial\mathcal{S}$ of radius $d_{t,s}(\theta)$. (The same reasoning applies to the other set $V_{t,s}^{-1}(\mathcal{S})\setminus V_{t}^{-1}(\mathcal{S})$.)

Consequently, applying the reasoning above to both terms on the right side of equation (D.1) gives,

Furthermore, if $\Theta\in\Delta$ is a random vector distributed according to $\mu$, then this upper bound may be written as

(This definition of $r_{t,s}$ will be used later on for handling the formal possibility that $\delta_{t,s}$ may be zero. Note that $\varepsilon_{t}/s$ is positive for all $s$ and $t$.) Clearly, the definition of $r_{t,s}$ gives

Next, we obtain an upper bound on this probability using a volume argument. Lemma E.6 in Appendix E shows that $V_{t,s}(\Theta)$ has a density, say $g_{t,s}$, with respect to $(k-1)$-dimensional Lebesgue measure. Also, Lemmas E.4 and E.7 imply that the random vector $V_{t,s}(\Theta)$ lies in $\Delta^{\circ}$ with probability 1 for all $t\geq 1$ and $s\geq 1$. Therefore,

To control the volume of the right side of the bound (D.6), it is convenient to consider the Hausdorff measure of $\partial\mathcal{S}$. Specifically, it is a fact from geometric measure theory that the $(k-2)$-dimensional Hausdorff measure of $\partial\mathcal{S}$, denoted $\mathcal{H}^{(k-2)}(\partial\mathcal{S})$, can be expressed as

In particular, $r_{t,s}$ is a sequence of positive numbers with $r_{t,s}\to 0$ as $t\to\infty$, and so when $s$ is fixed, this means

We now turn our attention to the factor $\sup_{\theta\in\Delta^{\circ}}g_{t,s}(\theta)$ in the bound (D.6). Lemma E.6 ensures there is a constant $c_{0}<\infty$, such that the following bound holds for every $s\geq 1$,

Combining lines (D.6), (D.8), and (D.9), we conclude that for every $s\geq 1$,

Finally, the proof is completed using the fact that $\kappa_{s}\to 0$ as $s\to\infty$, which is shown in Lemma E.8.∎

Appendix E Technical lemmas

To simplify the statements of the lemmas in this section, the notation in the previous appendices will be generally assumed without comment.

Recall $\bar{T}_{t}(\cdot):=\textstyle\frac{1}{t}\sum_{i=1}^{t}T_{i}(\cdot),$ and that we may express $\textsc{Err}_{t}$ and $\textsc{Err}_{t}^{\prime}$ as

for every fixed $(x,y)\in\mathcal{X}\times\mathcal{Y}$. This holds because when $x$ is fixed, we have $\mathcal{L}(\bar{Q}_{t}(x)|\mathcal{D})=\mathcal{L}(\bar{T}_{t}(x)|\mathcal{D})$, due to the first-order matching condition (3.1).∎

The Bernstein smoothing operator $\mathcal{B}_{s}$ satisfies the following properties.

For any function $h\in C$, the following limit holds

We refer to the book (devoreconstructive) for general background on these properties. The first property is given in Theorem 2.3 of (devoreconstructive, Chapter 1). The second property is proven after equation 1.7 of (devoreconstructive, Chapter 1).

Regarding the third property, if we let $u\in(0,1)$ be arbitrary, it is enough to show that the derivative $\textstyle\frac{d}{du}\mathcal{B}_{s}(h)(u)$ is strictly positive. To this end, it is shown in equation 2.2 of (devoreconstructive, Chapter 10) that $\mathcal{B}_{s}(h)$ satisfies

where we put $\underline{h}(u):=h(u+1/s)-h(u)$. If $h$ is non-decreasing and non-constant then

So, because all the terms $\underline{h}(\textstyle\frac{j}{s})$ are non-negative, at least one of them must be positive. In turn, equation (E.1) implies that $\textstyle\frac{d}{du}\mathcal{B}_{s}(h)(u)$ must be positive, since all the values $b_{j}(u;s-1)$ are positive for all $u\in(0,1)$. ∎

Let $\boldsymbol{D}$ denote the differentiation operator on univariate functions. If this operator is applied to a function $g$ with domain $$, then the domain of the derivative$\boldsymbol{D}(g)$is taken to be$(0,1)$.

Let $1_{(0,1)}$ be the indicator function of $(0,1)$. Then, for any fixed $s\geq 1$, we have the identity,

The identity (E.2) is a direct consequence of $\mathcal{B}_{s}(\text{id}_{})=\text{id}_{}$ from Lemma E.1. To prove the limit, first note that because the functions $b_{j}(\cdot;s)$ are polynomials on $$, the supremum$C_{s}:=\displaystyle\max_{0\leq j\leq s}\sup_{u\in(0,1)}|b_{j}^{\prime}(u;s)|$is finite. Hence, for fixed$s$,

For any fixed $s\geq 1$, there is a number $K_{s}\in[0,\infty)$ not depending on $t$ such that the inequality

holds for all $l=1,\dots,k-1$. Furthermore, as $t\to\infty$, we have the uniform limit

where $W_{s}=\boldsymbol{L}(\mathcal{B}_{s}(h))$.

Proof. From the definition of $W_{t,s}$ in equation (C.2), we have

The last expression is bounded in absolute value for every $l=1,\dots,k-1$, and every $t\geq 1$, by

which is finite by the uniform limit in Lemma (E.2). (Hence, the bound (E.4) is proved.) Lastly, the limit (E.5) follows from Lemma E.2 and the definition $W_{s}=\boldsymbol{L}(\mathcal{B}_{s}(h))$. ∎

For any $t\geq 1$ and $s\geq 1$, the following three statements are true:

The map $V_{t,s}:\Delta\to\Delta$ is bijective and continuous on $\Delta$, and is also $C^{1}$ on $\Delta^{\circ}$.

The Jacobian matrix $J(V_{t,s})(\theta)$ is non-singular for all $\theta\in\Delta^{\!\circ}$.

The inverse map $V_{t,s}^{-1}:\Delta\to\Delta$ is $C^{1}$ on $\Delta^{\circ}$.

It is simple to verify that $V_{t,s}$ is continuous on $\Delta$, and is $C^{1}$ on $\Delta^{\circ}$, due to the smoothness of $F_{t,s}$. To show that $V_{t,s}$ is bijective, it is enough to show that $F_{t,s}$ is bijective (due to the “inverses” property of $\boldsymbol{L}$ stated in Appendix A.1). In turn, the fact that $F_{t,s}=\mathcal{B}_{s}(F_{t})$ is bijective follows from part (c) of Lemma E.1.

To see that the Jacobian matrix $J(V_{t,s})(\theta)$ is non-singular for $\theta\in\Delta^{\circ}$, it is enough to check that $V_{t,s}^{-1}$ is $C^{1}$ on the set $\Delta^{\circ}$, because we may differentiate the identities

to establish the inverse of $J(V_{t,s})(\theta)$ via the chain rule. Using the “Inverses” property of $\boldsymbol{L}$, we have $V_{t,s}^{-1}=\boldsymbol{L}(F_{t,s}^{-1})$, and it follows that $V_{t,s}^{-1}$ will be $C^{1}$ as long as $F_{t,s}^{-1}$ is. Finally, the fact that $F_{t,s}^{-1}$ is $C^{1}$ follows from the strict monotonicity of $F_{t,s}$ and the univariate inverse function theorem.

The next lemma gives a uniform expansion for the determinant of $J(V_{t,s}^{-1})$ on $\Delta^{\circ}$.