Early stopping for kernel boosting algorithms: A general analysis with localized complexities

Introduction

While non-parametric models offer great flexibility, they can also lead to overfitting, and thus poor generalization performance. For this reason, it is well-understood that procedures for fitting non-parametric models must involve some form of regularization. When models are fit via a form of empirical risk minimization, the most classical form of regularization is based on adding some type of penalty to the objective function. An alternative form of regularization is based on the principle of early stopping, in which an iterative algorithm is run for a pre-specified number of steps, and terminated prior to convergence.

While the basic idea of early stopping is fairly old (e.g., ), recent years have witnessed renewed interests in its properties, especially in the context of boosting algorithms and neural network training (e.g., ). Over the past decade, a line of work has yielded some theoretical insight into early stopping, including works on classification error for boosting algorithms , $L^{2}$-boosting algorithms for regression , and similar gradient algorithms in reproducing kernel Hilbert spaces (e.g. ). A number of these papers establish consistency results for particular forms of early stopping, guaranteeing that the procedure outputs a function with statistical error that converges to zero as the sample size increases. On the other hand, there are relatively few results that actually establish rate optimality of an early stopping procedure, meaning that the achieved error matches known statistical minimax lower bounds. To the best of our knowledge, Bühlmann and Yu were the first to prove optimality for early stopping of $L^{2}$-boosting as applied to spline classes, albeit with a rule that was not computable from the data. Subsequent work by Raskutti et al. refined this analysis of $L^{2}$-boosting for kernel classes and first established an important connection to the localized Rademacher complexity; see also the related work with rates for particular kernel classes.

More broadly, relative to our rich and detailed understanding of regularization via penalization (e.g., see the books and papers for details), our understanding of early stopping regularization is not as well developed. Intuitively, early stopping should depend on the same bias-variance tradeoffs that control estimators based on penalization. In particular, for penalized estimators, it is now well-understood that complexity measures such as the localized Gaussian width, or its Rademacher analogue, can be used to characterize their achievable rates . Is such a general and sharp characterization also possible in the context of early stopping?

The main contribution of this paper is to answer this question in the affirmative for the early stopping of boosting algorithms for a certain class of regression and classification problems involving functions in reproducing kernel Hilbert spaces (RKHS). A standard way to obtain a good estimator or classifier is through minimizing some penalized form of loss functions of which the method of kernel ridge regression is a popular choice. Instead, we consider an iterative update involving the kernel that is derived from a greedy update. Borrowing tools from empirical process theory, we are able to characterize the “size” of the effective function space explored by taking $T$ steps, and then to connect the resulting estimation error naturally to the notion of localized Gaussian width defined with respect to this effective function space. This leads to a principled analysis for a broad class of loss functions used in practice, including the loss functions that underlie the $L^{2}$-boost, LogitBoost and AdaBoost algorithms, among other procedures.

The remainder of this paper is organized as follows. In Section 2, we provide background on boosting methods and reproducing kernel Hilbert spaces, and then introduce the updates studied in this paper. Section 3 is devoted to statements of our main results, followed by a discussion of their consequences for particular function classes in Section 4. We provide simulations that confirm the practical effectiveness of our stopping rules, and show close agreement with our theoretical predictions. In Section 5, we provide the proofs of our main results, with certain more technical aspects deferred to the appendices.

Background and problem formulation

In this section, we provide some necessary background on a gradient-type algorithm which is often referred to as boosting algorithm. We also discuss briefly about the reproducing kernel Hilbert spaces before turning to a precise formulation of the problem that is studied in this paper.

Consider a cost function $\phi:\times\rightarrow[0,\infty)$, where the non-negative scalar $\phi(y,\theta)$ denotes the cost associated with predicting $\theta$ when the true response is $y$. Some common examples of loss functions $\phi$ that we consider in later sections include:

the least-squares loss $\phi(y,\theta):\,=\frac{1}{2}(y-\theta)^{2}$ that underlies $L^{2}$-boosting ,

the logistic regression loss $\phi(y,\theta)=\ln(1+e^{-y\theta})$ that underlies the LogitBoost algorithm , and

the exponential loss $\phi(y,\theta)=\exp(-y\theta)$ that underlies the AdaBoost algorithm .

The least-squares loss is typically used for regression problems (e.g., ), whereas the latter two losses are frequently used in the setting of binary classification (e.g., ).

Given some loss function $\phi$, we define the population cost functional $f\mapsto\mathcal{L}(f)$ via

Note that with the covariates $\{x_{i}\}_{i=1}^{n}$ fixed, the functional $\mathcal{L}$ is a non-random object. Given some function space $\mathscr{F}$, the optimal functionAs clarified in the sequel, our assumptions guarantee uniqueness of $f^{*}$. minimizes the population cost functional—that is

Since we do not have access to the population distribution of the responses however, the computation of $f^{*}$ is impossible. Given our samples $\{Y_{i}\}_{i=1}^{n}$, we consider instead some procedure applied to the empirical loss

It is well-known that direct minimization of $\mathcal{L}_{n}$ over a sufficiently rich function class $\mathscr{F}$ may lead to overfitting. There are various ways to mitigate this phenomenon, among which the most classical method is to minimize the sum of the empirical loss with a penalty regularization term. Adjusting the weight on the regularization term allows for trade-off between fit to the data, and some form of regularity or smoothness in the fit. The behavior of such penalized of regularized estimation methods is now quite well understood (for instance, see the books and papers for more details).

In this paper, we study a form of algorithmic regularization, based on applying a gradient-type algorithm to $\mathcal{L}_{n}$ but then stopping it “early”—that is, after some fixed number of steps. Such methods are often referred to as boosting algorithms, since they involve “boosting” or improve the fit of a function via a sequence of additive updates (see e.g. ). Many boosting algorithms, among them AdaBoost , $L^{2}$-boosting and LogitBoost , can be understood as forms of functional gradient methods ; see the survey paper for further background on boosting. The way in which the number of steps is chosen is referred to as a stopping rule, and the overall procedure is referred to as early stopping of a boosting algorithm.

In more detail, a broad class of boosting algorithms generate a sequence $\{f^{t}\}_{t=0}^{\infty}$ via updates of the form

where the scalar $\{\alpha^{t}\}_{t=0}^{\infty}$ is a sequence of step sizes chosen by the user, the constraint $\|d\|_{\mathscr{F}}\leq 1$ defines the unit ball in a given function class $\mathscr{F}$, $\nabla\mathcal{L}_{n}(f)\in^{n}$ denotes the gradient taken at the vector \big{(}f(x_{1}),\ldots,f(x_{n})), and $\langle h,\,g\rangle$ is the usual inner product between vectors $h,g\in^{n}$. For non-decaying step sizes and a convex objective $\mathcal{L}_{n}$, running this procedure for an infinite number of iterations will lead to a minimizer of the empirical loss, thus causing overfitting. In order to illustrate this phenomenon, Figure 1 provides plots of the squared error \|f^{t}-f^{*}\|_{n}^{2}:\,=\frac{1}{n}\sum_{i=1}^{n}\big{(}f^{t}(x_{i})-f^{*}(x_{i})\big{)}^{2} versus the iteration number, for LogitBoost in panel (a) and AdaBoost in panel (b). See Section 4.2 for more details on exactly how these experiments were conducted.

In the plots in Figure 1, the dotted line indicates the minimum mean-squared error $\rho_{n}^{2}$ over all iterates of that particular run of the algorithm. Both plots are qualitatively similar, illustrating the existence of a “good” number of iterations to take, after which the MSE greatly increases. Hence a natural problem is to decide at what iteration $T$ to stop such that the iterate $f^{T}$ satisfies bounds of the form

with high probability. Here $f(n)\precsim g(n)$ indicates that $f(n)\leq cg(n)$ for some universal constant $c\in(0,\infty)$. The main results of this paper provide a stopping rule $T$ for which bounds of the form (5) do in fact hold with high probability over the randomness in the observed responses.

2 Reproducing Kernel Hilbert Spaces

for some coefficient vector $\omega\in^{n}$. Among those functions which achieve the infimum in expression (1), let us define $f^{*}$ as the one with the minimum Hilbert norm. This definition is equivalent to restricting $f^{*}$ to be in the linear subspace $\mathscr{H}_{n}$.

3 Boosting in kernel spaces

With a change of variable $d(x_{1}^{n})=\sqrt{n}\sqrt{K}z$ we then have

In this paper, we study the choice $g^{t}=\langle\nabla\mathcal{L}_{n}(f^{t}),\,d^{t}(x_{1}^{n})\rangle d^{t}$ in the boosting update (4), so that the function value iterates take the form

where $\alpha>0$ is a constant stepsize choice. Choosing $f^{0}(x_{1}^{n})=0$ ensures that all iterates $f^{t}(x_{1}^{n})$ remain in the range space of $K$.

In this paper we consider the following three error measures for an estimator $\widehat{f}$:

Main results

We now turn to the statement of our main results, beginning with the introduction of some regularity assumptions.

Recall from our earlier set-up that we differentiate between the empirical loss function $\mathcal{L}_{n}$ in expression (3), and the population loss $\mathcal{L}$ in expression (1). Apart from assuming differentiability of both functions, all of our remaining conditions are imposed on the population loss. Such conditions at the population level are weaker than their analogues at the empirical level.

For a given radius $r>0$, let us define the Hilbert ball around the optimal function $f^{*}$ as

Our analysis makes particular use of this ball defined for the radius $C_{\mathscr{H}}^{2}:\,=2\max\{\|f^{*}\|_{\mathscr{H}}^{2},~{}32,\sigma^{2}\}$ where the effective noise level $\sigma$ is defined in the sequel.

In addition to the least-squares cost, our theory also applies to losses $\mathcal{L}$ induced by scalar functions $\phi$ that satisfy the $\phi^{\prime}$-boundedness:

This condition holds with $B=1$ for the logistic loss for all $\mathcal{Y}$, and $B=\exp(2.5C_{\mathscr{H}})$ for the exponential loss for binary classification with $\mathcal{Y}=\{-1,1\}$, using our kernel boundedness condition. Note that whenever this condition holds with some finite $B$, we can always rescale the scalar loss $\phi$ by $1/B$ so that it holds with $B=1$, and we do so in order to simplify the statement of our results.

2 Upper bound in terms of localized Gaussian width

Our upper bounds involve a complexity measure known as the localized Gaussian width. In general, Gaussian widths are widely used to obtain risk bounds for least-squares and other types of $M$-estimators. In our case, we consider Gaussian complexities for “localized” sets of the form

with $f,g\in\mathscr{H}$. The Gaussian complexity localized at scale $\delta$ is given by

where $(w_{1},\ldots,w_{n})$ denotes an i.i.d. sequence of standard Gaussian variables.

An essential quantity in our theory is specified by a certain fixed point equation that is now standard in empirical process theory . Let us define the effective noise level

The critical radius $\delta_{n}$ is the smallest positive scalar such that

We note that past work on localized Rademacher and Gaussian complexity guarantees that there exists a unique $\delta_{n}>0$ that satisfies this condition, so that our definition is sensible.

Suppose that the sample size $n$ large enough such that $\delta_{n}\leq\frac{M}{m}$, and we compute the sequence $\{f^{t}\}_{t=0}^{\infty}$ using the update (9) with initialization $f^{0}=0$ and any step size $\alpha\in(0,\min\{\frac{1}{M},M\}]$. Then for any iteration T\in\big{\{}0,1,\ldots\lfloor\frac{m}{8M\delta_{n}^{2}}\rfloor\big{\}}, the averaged function estimate $\bar{f}^{T}$ satisfies the bounds

where both inequalities hold with probability at least $1-c_{1}\exp(-C_{2}\frac{m^{2}n\delta_{n}^{2}}{\sigma^{2}})$.

A few comments about the constants in our statement: in all cases, constants of the form $c_{j}$ are universal, whereas the capital $C_{j}$ may depend on parameters of the joint distribution and population loss $\mathcal{L}$. In Theorem 1, we have the explicit value $C_{2}=\{\frac{m^{2}}{\sigma^{2}},1\}$ and $C^{2}$ is proportional to the quantity $2\max\{\|f^{*}\|_{\mathscr{H}}^{2},~{}32,~{}\sigma^{2}\}$. While inequalities (15a) and (15b) are stated as high probability results, similar bounds for expected loss (over the response $y_{i}$, with the design fixed) can be obtained by a simple integration argument.

In order to gain intuition for the claims in the theorem, note that apart from factors depending on $(m,M)$, the first term $\frac{1}{\alpha mT}$ dominates the second term $\frac{\delta_{n}^{2}}{m^{2}}$ whenever $T\lesssim 1/\delta_{n}^{2}$. Consequently, up to this point, taking further iterations reduces the upper bound on the error. This reduction continues until we have taken of the order $1/\delta_{n}^{2}$ many steps, at which point the upper bound is of the order $\delta_{n}^{2}$.

More precisely, suppose that we perform the updates with step size $\alpha=\frac{m}{M}$; then, after a total number of $\tau:\,=\frac{1}{\delta_{n}^{2}\max\{8,M\}}$ many iterations, the extension of Theorem 1 to expectations guarantees that the mean squared error is bounded as

where $C^{\prime}$ is another constant depending on $C_{\mathscr{H}}$. Here we have used the fact that $M\geq m$ in simplifying the expression. It is worth noting that guarantee (16) matches the best known upper bounds for kernel ridge regression (KRR)—indeed, this must be the case, since a sharp analysis of KRR is based on the same notion of localized Gaussian complexity (e.g. ) . Thus, our results establish a strong parallel between the algorithmic regularization of early stopping, and the penalized regularization of kernel ridge regression. Moreover, as will be clarified in Section 3.3, under suitable regularity conditions on the RKHS, the critical squared radius $\delta_{n}^{2}$ also acts as a lower bound for the expected risk, meaning that our upper bounds are not improvable in general.

Note that the critical radius $\delta_{n}^{2}$ only depends on our observations $\{(x_{i},y_{i})\}_{i=1}^{n}$ through the solution of inequality (14). In many cases, it is possible to compute and/or upper bound this critical radius, so that a concrete and valid stopping rule can indeed by calculated in advance. In Section 4, we provide a number of settings in which this can be done in terms of the eigenvalues $\{\mu_{j}\}_{j=1}^{n}$ of the normalized kernel matrix.

2.2 Consequences for random design regression

In order to state an upper bound on this error, we introduce a population analogue of the critical radius $\delta_{n}$, which we denote by $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{n}$. Consider the set

It is analogous to the previously defined set $\mathcal{E}(\delta,1)$, except that the empirical norm $\|\cdot\|_{n}$ has been replaced by the population version. The population Gaussian complexity localized at scale $\delta$ is given by

with probability at least $1-c_{1}\exp(-C_{2}\frac{m^{2}n\delta_{n}^{2}}{\sigma^{2}})$ over the random samples.

The proof of Corollary 1 follows directly from standard empirical process theory bounds on the difference between empirical risk $\|\bar{f}^{T}-f^{*}\|_{n}^{2}$ and population risk $\|\bar{f}^{T}-f^{*}\|_{2}^{2}$. In particular, it can be shown that $\|\cdot\|_{2}$ and $\|\cdot\|_{n}$ norms differ only by a factor proportion to $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{n}$. Furthermore, one can show that the empirical critical quantity $\delta_{n}$ is bounded by the population $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{n}$. By combining both arguments the corollary follows. We refer the reader to the papers for further details on such equivalences.

It is worth comparing this guarantee with the past work of Raskutti et al. , who analyzed the kernel boosting iterates of the form (9), but with attention restricted to the special case of the least-squares loss. Their analysis was based on first decomposing the squared error into bias and variance terms, then carefully relating the combination of these terms to a particular bound on the localized Gaussian complexity (see equation (19) below). In contrast, our theory more directly analyzes the effective function class that is explored by taking $T$ steps, so that the localized Gaussian width (18) appears more naturally. In addition, our analysis applies to a broader class of loss functions.

In the case of reproducing kernel Hilbert spaces, it is possible to sandwich the localized Gaussian complexity by a function of the eigenvalues of the kernel matrix. Mendelson provides this argument in the case of the localized Rademacher complexity, but similar arguments apply to the localized Gaussian complexity. Letting $\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{n}\geq 0$ denote the ordered eigenvalues of the normalized kernel matrix $K$, define the function

Up to a universal constant, this function is an upper bound on the Gaussian width \mathcal{G}_{n}\big{(}\mathcal{E}(\delta,1)\big{)} for all $\delta\geq 0$, and up to another universal constant, it is also a lower bound for all $\delta\geq\frac{1}{\sqrt{n}}$.

3 Achieving minimax lower bounds

In this section, we show that the upper bound (16) matches known minimax lower bounds on the error, so that our results are unimprovable in general. We establish this result for the class of regular kernels, as previously defined by Yang et al. , which includes the Gaussian and Sobolev kernels as special cases.

The class of regular kernels is defined as follows. Let $\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{n}\geq 0$ denote the ordered eigenvalues of the normalized kernel matrix $K$, and define the quantity $d_{n}:\,=\operatorname{argmin}_{j=1,\dots,n}\{\mu_{j}\leq\delta_{n}^{2}\}$. A kernel is called regular whenever there is a universal constant $c$ such that the tail sum satisfies $\sum_{j=d_{n}+1}^{n}\mu_{j}\leq c\,d_{n}\delta_{n}^{2}$. In words, the tail sum of the eigenvalues for regular kernels is roughly on the same or smaller scale as the sum of the eigenvalues bigger than $\delta_{n}^{2}$.

For such kernels and under the Gaussian observation model ($Y_{i}\sim N(f^{*}(x_{i}),\sigma^{2})$), Yang et al. prove a minimax lower bound involving $\delta_{n}$. In particular, they show that the minimax risk over the unit ball of the Hilbert space is lower bounded as

Comparing the lower bound (20) with upper bound (16) for our estimator $\bar{f}^{T}$ stopped after $O(1/\delta_{n}^{2})$ many steps, it follows that the bounds proven in Theorem 1 are unimprovable apart from constant factors.

We now state a generalization of this minimax lower bound, one which applies to a sub-class of generalized linear models, or GLM for short. In these models, the conditional distribution of the observed vector $Y=(Y_{1},\ldots,Y_{n})$ given \big{(}f^{*}(x_{1}),\ldots,f^{*}(x_{n})\big{)} takes the form

where $s(\sigma)$ is a known scale factor and $\Phi:\to$ is the cumulant function of the generalized linear model. As some concrete examples:

The linear Gaussian model is recovered by setting $s(\sigma)=\sigma^{2}$ and $\Phi(t)=t^{2}/2$.

The logistic model for binary responses $y\in\{-1,1\}$ is recovered by setting $s(\sigma)=1$ and $\Phi(t)=\log(1+\exp(t))$.

Our minimax lower bound applies to the class of GLMs for which the cumulant function $\Phi$ is differentiable and has uniformly bounded second derivative $|\Phi^{\prime\prime}|\leq L$. This class includes the linear, logistic, multinomial families, among others, but excludes (for instance) the Poisson family. Under this condition, we have the following:

Suppose that we are given i.i.d. samples $\{y_{i}\}_{i=1}^{n}$ from a GLM (21) for some function $f^{*}$ in a regular kernel class with $\|f^{*}\|_{\mathscr{H}}\leq 1$. Then running $T:\,=\lfloor\frac{1}{\delta_{n}^{2}\max\{8,M\}}\rfloor$ iterations with step size $\alpha\in(0,\min\{\frac{1}{M},M\}]$ and $f^{0}=0$ yields an estimate $\bar{f}^{T}$ such that

At a high level, the statement in Corollary 2 shows that early stopping prevents us from overfitting to the data; in particular, using the stopping time $T$ yields an estimate that attains the optimal balance between bias and variance.

Consequences for various kernel classes

In this section, let us consider two broad types of eigen-decay:

$\gamma$-exponential decay: For some $\gamma>0$, the kernel matrix eigenvalues satisfy a decay condition of the form $\mu_{j}\leq c_{1}\exp(-c_{2}j^{\gamma})$, where $c_{1},c_{2}$ are universal constants. Examples of kernels in this class include the Gaussian kernel, which for the Lebesgue measure satisfies such a bound with $\gamma=2$ (real line) or $\gamma=1$ (compact domain).

$\beta$-polynomial decay: For some $\beta>1/2$, the kernel matrix eigenvalues satisfy a decay condition of the form $\mu_{j}\leq c_{1}j^{-2\beta}$, where $c_{1}$ is a universal constant. Examples of kernels in this class include the $k^{th}$-order Sobolev spaces for some fixed integer $k\geq 1$ with Lebesgue measure on a bounded domain. We consider Sobolev spaces that consist of functions that have $k^{th}$-order weak derivatives $f^{(k)}$ being Lebesgue integrable and $f(0)=f^{(1)}(0)=\dots=f^{(k-1)}(0)=0$. For such classes, the $\beta$-polynomial decay condition holds with $\beta=k$.

Given eigendecay conditions of these types, it is possible to compute an upper bound on the critical radius $\delta_{n}$. In particular, using the fact that the function $\mathcal{R}$ from equation (19) is an upper bound on the function \mathcal{G}_{n}\big{(}\mathcal{E}(\delta,1)\big{)}, we can show that for $\gamma$-exponentially decaying kernels, we have $\delta_{n}^{2}\precsim\frac{(\log n)^{1/\gamma}}{n}$, whereas for $\beta$-polynomial kernels, we have $\delta_{n}^{2}\precsim n^{-\frac{2\beta}{2\beta+1}}$ up to universal constants. Combining with our Theorem 1, we obtain the following result:

For kernels with $\gamma$-exponential eigen-decay, we have

For kernels with $\beta$-polynomial eigen-decay, we have

See Section 5.3 for the proof of Corollary 3.

To the best of our knowledge, this result is the first to show non-asymptotic and optimal statistical rates for the $\|\cdot\|_{n}^{2}$-error when early stopping LogitBoost or AdaBoost with an explicit dependence of the stopping rule on $n$. Our results also yield similar guarantees for $L^{2}$-boosting, as has been established in past work . Note that we can observe a similar trade-off between computational efficiency and statistical accuracy as in the case of kernel least-squares regression : although larger kernel classes (e.g. Sobolev classes) yield higher estimation errors, boosting updates reach the optimum faster than for a smaller kernel class (e.g. Gaussian kernels).

2 Numerical experiments

We now describe some numerical experiments that provide illustrative confirmations of our theoretical predictions. While we have applied our methods to various kernel classes, in this section, we present numerical results for the first-order Sobolev kernel as two typical examples for exponential and polynomial eigen-decay kernel classes.

Figure 2 shows plots of the mean-squared error $\|\bar{f}^{T}-f^{*}\|_{n}^{2}$ over the sample size $n$ averaged over $40$ trials, for the gold standard $T=G$ and stopping rules based on $T=(7n)^{\kappa}$ for different choices of $\kappa$. Error bars correspond to the standard errors computed from our simulations. Panel (a) shows the behavior for $L^{2}$-boosting, whereas panel (b) shows the behavior for LogitBoost.

Note that both plots are qualitatively similar and that the theoretically derived stopping rule $T=(7n)^{\kappa}$ with $\kappa^{*}=2/3=0.67$, while slightly worse than the Gold standard, tracks its performance closely. We also performed simulations for some “bad” stopping rules, in particular for an exponent $\kappa$ not equal to $\kappa^{*}=2/3$, indicated by the green and black curves. In the log scale plots in Figure 3 we can clearly see that for $\kappa\in\{0.33,1\}$ the performance is indeed much worse, with the difference in slope even suggesting a different scaling of the error with the number of observations $n$. Recalling our discussion for Figure 1, this phenomenon likely occurs due to underfitting and overfitting effects. These qualitative shifts are consistent with our theory.

Proof of main results

In this section, we present the proofs of our main results. The technical details are deferred to Appendix A.

The proof of our main theorem is based on a sequence of lemmas, all of which are stated with the assumptions of Theorem 1 in force. The first lemma establishes a bound on the empirical norm $\|\cdot\|_{n}$ of the error $\Delta^{t+1}:\,=\theta^{t+1}-\theta^{*}$, provided that its Hilbert norm is suitably controlled.

For any stepsize $\alpha\in(0,~{}\frac{1}{M}]$ and any iteration $t$ we have

See Section A.1 for the proof of this claim.

The second term on the right-hand side of the bound (1) involves the difference between the population and empirical gradient operators. Since this difference is being evaluated at the random points $\Delta^{t}$ and $\Delta^{t+1}$, the following lemma establishes a form of uniform control on this term.

See Section A.2 for the proof of Lemma 2.

Note that Lemma 1 applies only to error iterates with a bounded Hilbert norm. Our last lemma provides this control for some number of iterations:

There are constants $(C_{1},C_{2})$ independent of $n$ such that for any step size \alpha\in\big{(}0,\min\{M,\frac{1}{M}\}\big{]}, we have

with probability at least $1-C_{1}\exp(-C_{2}n\delta_{n}^{2})$, where $C_{2}=\max\{\frac{m^{2}}{\sigma^{2}},1\}$.

See Section A.3 for the proof of this lemma which also uses Lemma 2.

Taking these lemmas as given, we now complete the proof of the theorem. We first condition on the event $\mathcal{E}$ from Lemma 2, so that we may apply the bound (5.1). We then fix some iterate $t$ such that $t<\frac{m}{8M\delta_{n}^{2}}-1$, and condition on the event that the bound (26) in Lemma 3 holds, so that we are guaranteed that $\|\Delta^{t+1}\|_{\mathscr{H}}\leq C_{\mathscr{H}}$. We then split the analysis into two cases:

First, suppose that $\|\Delta^{t+1}\|_{n}\leq\delta_{n}C_{\mathscr{H}}$. In this case, inequality (15b) holds directly.

Otherwise, we may assume that $\|\Delta^{t+1}\|_{n}>\delta_{n}\|\Delta^{t+1}\|_{\mathscr{H}}$. Applying the bound (5.1) with the choice $(\widetilde{\Delta},\Delta)=(\Delta^{t},\Delta^{t+1})$ yields

Substituting inequality (27) back into equation (1) yields

where we have introduced the shorthand notation D^{t}:\,=\frac{1}{2\alpha}\Big{\{}\|\Delta^{t}\|_{\mathscr{H}}^{2}-\|\Delta^{t+1}\|_{\mathscr{H}}^{2}\Big{\}}, as well as $\gamma=\frac{1}{2}-\frac{1}{c_{3}}$

Equation (28) defines a quadratic inequality with respect to $\|\Delta^{t+1}\|_{n}$; solving it and making use of the inequality $(a+b)^{2}\leq 2a^{2}+2b^{2}$ yields the bound

for some universal constant $c$. By telescoping inequality (29), we find that

so that inequality (15b) follows from the bound (30).

On the other hand, by the smoothness assumption, we have

2 Proof of Corollary 2

Similar to the proof of Theorem 1 in Yang et al. , a generalization can be shown using a standard argument of Fano’s inequality. By definition of the transformed parameter $\theta=DU\alpha$ with $K=U^{T}DU$, we have for any estimator $\widehat{f}=\sqrt{n}U^{T}\theta$ that $\|\widehat{f}-f^{*}\|_{n}^{2}=\|\theta-{\theta^{*}}\|_{2}^{2}$. Therefore our goal is to lower bound the Euclidean error $\|\theta-{\theta^{*}}\|_{2}$ of any estimator of ${\theta^{*}}.$ Borrowing Lemma 4 in Yang et al. , there exists $\delta/2$-packing of the set $B=\{\theta\in^{n}\mid\|D^{-1/2}\theta\|_{2}\leq 1\}$ of cardinality $M=e^{d_{n}/64}$ with $d_{n}:\,=\arg\min_{j=1,\dots,n}\{\mu_{j}\leq\delta_{n}^{2}\}.$ This is done through packing the following subset of $B$

Let us denote the packing set by $\{\theta^{1},\ldots,\theta^{M}\}.$ Since $\theta\in\mathcal{E}(\delta)$, by simple calculation, we have $\|\theta^{i}\|_{2}\leq\delta$.

where $I(y_{1}^{n};Z)$ is the mutual information between the samples $Y$ and the random index $Z$.

Since each $\|\theta^{i}\|_{2}\leq\delta$, triangle inequality yields $\|\theta_{i}-\theta_{j}\|_{2}\leq 2\delta$ for all $i\neq j$. It is therefore guaranteed that

Therefore, similar to Yang et al. , following by the fact that the kernel is regular and hence $s(\sigma)d_{n}\geq cn\delta_{n}^{2}$, any estimator $\widehat{f}$ has prediction error lower bounded as

Corollary 2 thus follows using the upper bound in Theorem 1.

Direct calculations of the KL-divergence yield

which follows by assumption on $\Phi$ having a uniformly bounded second derivative. Putting the above inequality with inequality (5.2) establishes our claim (32).

3 Proof of Corollary 3

With $D:\,=C_{\mathscr{H}}+\|\theta^{*}\|_{\mathscr{H}}$, the logistic regression cost function satisfies the $m$-$M$-condition with parameters

The AdaBoost cost function satisfies the $m$-$M$-condition with parameters

See Section A.4 for the proof of Lemma 4.

If the kernel eigenvalues satisfy a decay condition of the form $\mu_{j}\leq c_{1}\exp(-c_{2}j^{\gamma})$, where $c_{1},c_{2}$ are universal constants, the function $\mathcal{R}$ from equation (19) can be upper bounded as

where $k$ is the smallest integer such that $c_{1}\exp(-c_{2}k^{\gamma})<\delta^{2}$. Since the localized Gaussian width \mathcal{G}_{n}\big{(}\mathcal{E}_{n}(\delta,1)\big{)} can be sandwiched above and below by multiples of $\mathcal{R}(\delta)$, some algebra shows that the critical radius scales as $\delta_{n}^{2}\asymp\frac{n}{\log(n)^{1/\gamma}\sigma^{2}}$.

Consequently, if we take $T\asymp\frac{\log(n)^{1/\gamma}\sigma^{2}}{n}$ steps, then Theorem 1 guarantees that the averaged estimator $\bar{\theta}^{T}$ satisfies the bound

with probability $1-c_{1}\text{exp}(-c_{2}m^{2}\log^{1/\gamma}n)$.

Now suppose that the kernel eigenvalues satisfy a decay condition of the form $\mu_{j}\leq c_{1}j^{-2\beta}$ for some $\beta>1/2$ and constant $c_{1}$. In this case, a direct calculation yields the bound

where $k$ is the smallest integer such that $c_{2}k^{-2}<\delta^{2}$. Combined with upper bound $c_{2}\sum_{j=k+1}^{n}j^{-2}\leq c_{2}\int_{k+1}j^{-2}\leq k\delta^{2}$, we find that the critical radius scales as $\delta_{n}^{2}\asymp n^{-2\beta/(1+2\beta)}$.

Consequently, if we take $T\asymp n^{-2\beta/(1+2\beta)}$ many steps, then Theorem 1 guarantees that the averaged estimator $\bar{\theta}^{T}$ satisfies the bound

with probability at least $1-c_{1}\text{exp}(-c_{2}m^{2}(\frac{n}{\sigma^{2}})^{1/(2\beta+1)})$.

Discussion

In this paper, we have proven non-asymptotic bounds for early stopping of kernel boosting for a relatively broad class of loss functions. These bounds allowed us to propose simple stopping rules which, for the class of regular kernel functions , yield minimax optimal rates of estimation. Although the connection between early stopping and regularization has long been studied and explored in the theoretical literature and applications alike, to the best of our knowledge, this paper is the first one to establish a general relationship between the statistical optimality of stopped iterates and the localized Gaussian complexity. This connection is important, because this localized Gaussian complexity measure, as well as its Rademacher analogue, are now well-understood to play a central role in controlling the behavior of estimators based on regularization .

There are various open questions suggested by our results. The stopping rules in this paper depend on the eigenvalues of the empirical kernel matrix; for this reason, they are data-dependent and computable given the data. However, in practice, it would be desirable to avoid the cost of computing all the empirical eigenvalues. Can fast approximation techniques for kernels be used to approximately compute our optimal stopping rules? Second, our current theoretical results apply to the averaged estimator $\bar{f}^{T}$. We strongly suspect that the same results apply to the stopped estimator $f^{T}$, but some new ingredients are required to extend our proofs.

This work was partially supported by DOD Advanced Research Projects Agency W911NF-16-1-0552, National Science Foundation grant NSF-DMS-1612948, and Office of Naval Research Grant DOD-ONR-N00014.

References

Appendix A Proof of technical lemmas

Recalling that $K^{\dagger}$ denotes the pseudoinverse of $K$, our proof is based on the linear transformation

If we transform this update on $z$ back to an equivalent one on $\theta$ by multiplying both sides by $\sqrt{n}\sqrt{K}$, we see that ordinary gradient descent on $\mathcal{J}_{n}$ is equivalent to the kernel boosting update $\theta^{t+1}=\theta^{t}-\alpha nK\nabla\mathcal{L}_{n}(\theta^{t})$.

Our goal is to analyze the behavior of the update (36) in terms of the population cost $\mathcal{J}(z^{t})$. Thus, our problem is one of analyzing a noisy form of gradient descent on the function $\mathcal{J}$, where the noise is induced by the difference between the empirical gradient operator $\nabla\mathcal{J}_{n}$ and the population gradient operator $\nabla\mathcal{J}$.

Recall that the $\mathcal{L}$ is $M$-smooth by assumption. Since the kernel matrix $K$ has been normalized to have largest eigenvalue at most one, the function $\mathcal{J}$ is also $M$-smooth, whence

Morever, since the function $\mathcal{J}$ is convex, we have $\mathcal{J}(z^{*})\geq\mathcal{J}(z^{t})+\langle\nabla\mathcal{J}(z^{t}),\,z^{*}-z^{t}\rangle$, whence

Now define the difference of the squared errors V^{t}:\,=\frac{1}{2}\Big{\{}\|z^{t}-z^{*}\|_{2}^{2}-\|z^{t+1}-z^{*}\|_{2}^{2}\Big{\}}. By some simple algebra, we have

Substituting back into equation (37) yields

where we have used the fact that $\frac{1}{\alpha}\geq M$ by our choice of stepsize $\alpha$.

Finally, we transform back to the original variables $\theta=\sqrt{n}\sqrt{K}z$, using the relation $\nabla\mathcal{J}(z)=\sqrt{n}\sqrt{K}\nabla\mathcal{L}(\theta)$, so as to obtain the bound

Note that the optimality of $\theta^{*}$ implies that $\nabla\mathcal{L}(\theta^{*})=0$. Combined with $m$-strong convexity, we are guaranteed that $\frac{m}{2}\|\Delta^{t+1}\|_{n}^{2}\leq\mathcal{L}(\theta^{t+1})-\mathcal{L}(\theta^{*})$, and hence

A.2 Proof of Lemma 2

We split our proof into two cases, depending on whether we are dealing with the least-squares loss $\phi(y,\theta)=\frac{1}{2}(y-\theta)^{2}$, or a classification loss with uniformly bounded gradient ($\|\phi^{\prime}\|_{\infty}\leq 1$).

The least-squares loss is $m$-strongly convex with $m=M=1$. Moreover, the difference between the population and empirical gradients can be written as $\nabla\mathcal{L}(\theta^{*}+\widetilde{\Delta})-\nabla\mathcal{L}_{n}(\theta^{*}+\widetilde{\Delta})=\frac{\sigma}{n}(w_{1},\ldots,w_{n})$, where the random variables $\{w_{i}\}_{i=1}^{n}$ are i.i.d. and sub-Gaussian with parameter $1$. Consequently, we have

Under these conditions, one can show (see for reference) that

which implies that Lemma 2 holds with $c_{3}=16$.

A.2.2 Gradient-bounded ϕitalic-ϕ\phi-functions

Applying the bound (5.1) to $g$ and then multiplying both sides by $\|\Delta\|_{\mathscr{H}}$, we obtain

where the second inequality uses the fact that $\|\Delta\|_{\mathscr{H}}>1$ by assumption.

In order to establish the bound (5.1) for functions with $\|g\|_{\mathscr{H}}=1$, we first prove it uniformly over the set $\{g\mid\|g\|_{\mathscr{H}}=1,\quad\|g\|_{n}\leq t\}$, where $t>1$ is a fixed radius (of course, we restrict our attention to those radii $t$ for which this set is non-empty.) We then extend the argument to one that is also uniform over the choice of $t$ by a “peeling” argument.

The following two lemmas, respectively, bound the mean of this random variable, and its deviations above the mean:

For any $t>0$, the mean is upper bounded as

There are universal constants $(c_{1},c_{2})$ such that

See Appendices A.2.3 and A.2.4 for the proofs of these two claims.

Equipped with Lemmas 5 and 6, we now prove inequality (5.1). We divide our argument into two cases:

We first prove inequality (5.1) for $t=\delta_{n}$. From Lemma 5, we have

where inequality (i) follows from the definition of $\delta_{n}$ in inequality (14). Setting $\alpha=\delta_{n}^{2}$ in expression (41) yields

which establishes the claim for $t=\delta_{n}$.

On the other hand, for any $t>\delta_{n}$, we have

Note that the precise values of the universal constants $c_{2}$ may change from line to line throughout this section.

Equipped with the tail bounds (43) and (44), we are now ready to complete the peeling argument. Let $\mathcal{A}$ denote the event that the bound (5.1) is violated for some function $g\in\partial\mathscr{H}$ with $\|g\|_{\mathscr{H}}=1$. For real numbers $0\leq a<b$, let $\mathcal{A}(a,b)$ denote the event that it is violated for some function such that $\|g\|_{n}\in[a,b]$, and $\|g\|_{\mathscr{H}}=1$. For $k=0,1,2,\ldots$, define $t_{k}=2^{k}\delta_{n}$. We then have the decomposition $\mathcal{E}=(0,t_{0})\cup(\bigcup_{k=0}^{\infty}\mathcal{A}(t_{k},t_{k+1}))$ and hence by union bound,

where step (i) uses the $\|g\|_{n}\geq t_{k}$ and step (ii) uses the fact that $t_{k+1}=2t_{k}.$ This lower bound implies that $\mathcal{Z}_{n}(t_{k+1})>t_{k+1}\delta_{n}+\frac{t_{k+1}^{2}m}{4c_{3}}$ and applying the tail bound (44) yields

Substituting this inequality and our earlier bound (43) into equation (45) yields

where the reader should recall that the precise values of universal constants may change from line-to-line. This concludes the proof of Lemma 2.

A.2.3 Proof of Lemma 5

Recalling the definitions (1) and (3) of $\mathcal{L}$ and $\mathcal{L}_{n}$, we can write

Thus, we can restrict our attention to vectors $\Delta,\widetilde{\Delta}$ with $\|\Delta\|_{\infty},\|\widetilde{\Delta}\|_{\infty}\leq 2C_{\mathscr{H}}$ from hereonwards.

Letting $\{\varepsilon_{i}\}_{i=1}^{n}$ denote an i.i.d. sequence of Rademacher variables, define the symmetrized variable

where the second inequality follows by applying the Rademacher contraction inequality , using the fact that $\|\phi^{\prime}\|_{\infty}\leq 1$ for the first term, and $\|\Delta\|_{\infty}\leq 2C_{\mathscr{H}}$ for the second term.

where step (i) follows since each function $\varphi_{i}$ is $M$-Lipschitz by assumption; and step (ii) follows since the Gaussian complexity upper bounds the Rademacher complexity up to a factor of $\sqrt{\frac{\pi}{2}}$. Similarly, we have

and putting together the pieces yields the claim.

A.2.4 Proof of Lemma 6

where we have used the facts that $\|\phi^{\prime}\|_{\infty}\leq 1$ and $\|\Delta\|_{n}\leq t$. By Ledoux’s concentration for convex and Lipschitz functions , we have

Since the right-hand side does not involve $\{y_{i}\}_{i=1}^{n}$, the same bound holds unconditionally over the randomness in both the Rademacher variables and the sequence $\{y_{i}\}_{i=1}^{n}$. Consequently, the claimed bound (41) follows, with suitable redefinitions of the universal constants.

A.3 Proof of Lemma 3

We first require an auxiliary lemma, which we state and prove in the following section. We then prove Lemma 3 in Section A.3.2.

The following result relates the Hilbert norm of the error to the difference between the empirical and population gradients:

For any convex and differentiable loss function $\mathcal{L}$, the kernel boosting error $\Delta^{t+1}:\,=\theta^{t+1}-\theta^{*}$ satisfies the bound

Recall that $\|\Delta^{t}\|_{\mathscr{H}}^{2}=\|\theta^{t}-\theta^{*}\|_{\mathscr{H}}^{2}=\|z^{t}-z^{*}\|_{2}^{2}$ by definition of the Hilbert norm. Let us define the population update operator $G$ on the population function $\mathcal{J}$ and the empirical update operator $G_{n}$ on $\mathcal{J}_{n}$ as

Since $\mathcal{J}$ is convex and smooth, it follows from standard arguments in convex optimization that $G$ is a non-expansive operator—viz.

In addition, we note that the vector $z^{*}$ is a fixed point of $G$—that is, $G(z^{*})=z^{*}$. From these ingredients, we have

where step (i) follows by applying the Cauchy-Schwarz to control the inner product, and step (ii) follows since $\Delta^{t+1}=\sqrt{n}\sqrt{K}(z^{t+1}-z^{*})$, and the square root kernel matrix $\sqrt{K}$ is symmetric. ∎

A.3.2 Proof of Lemma 3

We now prove Lemma 3. The argument makes use of Lemmas 1 and 2 combined with Lemma 7.

In order to prove inequality (26), we follow an inductive argument. Instead of proving (26) directly, we prove a slightly stronger relation which implies it, namely

Here $\widetilde{\gamma}$ and $c_{3}$ are constants linked by the relation

We claim that it suffices to prove that the error iterates $\Delta^{t+1}$ satisfy the inequality (50). Indeed, if we take inequality (50) as given, then we have

where we used the definition $C_{\mathscr{H}}^{2}=2\max\{\|\theta^{*}\|_{\mathscr{H}}^{2},~{}32\}$. Thus, it suffices to focus our attention on proving inequality (50).

For $t=0$, it is trivially true. Now let us assume inequality (50) holds for some $t\leq\frac{m}{8M\delta_{n}^{2}}$, and then prove that it also holds for step $t+1$.

If $\|\Delta^{t+1}\|_{\mathscr{H}}<1$, then inequality (50) follows directly. Therefore, we can assume without loss of generality that $\|\Delta^{t+1}\|_{\mathscr{H}}\geq 1$.

We break down the proof of this induction into two steps:

First, we show that $\|\Delta^{t+1}\|_{\mathscr{H}}\leq 2C_{\mathscr{H}}$ so that Lemma 2 is applicable.

Second, we show that the bound (50) holds and thus in fact $\|\Delta^{t+1}\|_{\mathscr{H}}\leq C_{\mathscr{H}}$.

with $C_{2}=\max\{\frac{m^{2}}{\sigma^{2}},1\}$.

Applying the triangle inequality yields the bound

where the population update operator $G$ was previously defined (A.3.1), and observed to be non-expansive (49). From this non-expansiveness, we find that

Observe that because of uniform boundedness of the kernel by one, the quantity $T$ can be bounded as

where we have used the fact that $\alpha\leq m/M<1\leq\frac{C_{\mathscr{H}}}{2}.$ For least-squares $\phi$ we instead have

Putting together the pieces yields that $\|\Delta^{t+1}\|_{\mathscr{H}}\leq 2C_{\mathscr{H}}$, as claimed.

We are now ready to complete the induction step for proving inequality (50) using Lemma 1 and Lemma 2 since $\|\Delta^{t+1}\|_{\mathscr{H}}\geq 1$. We split the argument into two cases separately depending on whether or not $\|\Delta^{t+1}\|_{\mathscr{H}}\delta_{n}\geq\|\Delta^{t+1}\|_{n}$. In general we can assume that $\|\Delta^{t+1}\|_{\mathscr{H}}>\|\Delta^{t}\|_{\mathscr{H}}$, otherwise the induction inequality (50) satisfies trivially.

When $\|\Delta^{t+1}\|_{\mathscr{H}}\delta_{n}\geq\|\Delta^{t+1}\|_{n}$, inequality (5.1) implies that

Combining Lemma 7 and inequality (A.3.2), we obtain

where the last inequality uses the fact that $\|\Delta^{t+1}\|_{n}\leq\delta_{n}\|\Delta^{t+1}\|_{\mathscr{H}}.$

When $\|\Delta^{t+1}\|_{\mathscr{H}}\delta_{n}<\|\Delta^{t+1}\|_{n}$, we use our assumption $\|\Delta^{t+1}\|_{\mathscr{H}}\geq\|\Delta^{t}\|_{\mathscr{H}}$ together with Lemma 7 and inequality (5.1) which guarantee that

Using the elementary inequality $2ab\leq a^{2}+b^{2}$, we find that

where in the final step, we plug in the constants $\widetilde{\gamma},c_{3}$ which satisfy equation (51).

where step (i) again uses $2ab\leq a^{2}+b^{2}$. Thus, we have $\frac{m}{4}\|\Delta^{t+1}\|_{n}^{2}\leq D^{t}+\frac{1}{\widetilde{\gamma}m}\delta_{n}^{2}$. Together with expression (A.3.2), we find that

By combining the two previous cases, we arrive at the bound

where $\kappa:\,=\frac{1}{(1-\alpha\delta_{n}^{2}\frac{m}{c_{3}})}$ and we used that $\alpha\leq\min\{\frac{1}{M},M\}$.

Now it is only left for us to show that with the constant $c_{3}$ chosen such that $\widetilde{\gamma}=\frac{1}{32}-\frac{1}{4c_{3}}=1/C_{\mathscr{H}}^{2}$, we have

Define the function $f:(0,C_{\mathscr{H}}]\rightarrow$ via $f(\xi):\,=\kappa^{2}(\xi+4\alpha\delta_{n}^{2})^{2}-\xi^{2}-\frac{4M}{\widetilde{\gamma}m}\delta_{n}^{2}$. Since $\kappa\geq 1$, in order to conclude that $f(\xi)<0$ for all $\xi\in(0,C_{\mathscr{H}}]$, it suffices to show that $\operatorname{argmin}_{x\in}f(x)<0$ and $f(C_{\mathscr{H}})<0$. The former is obtained by basic algebra and follows directly from $\kappa\geq 1$. For the latter, since $\widetilde{\gamma}=\frac{1}{32}-\frac{1}{4c_{3}}=1/C_{\mathscr{H}}^{2}$, $\alpha<\frac{1}{M}$ and $\delta_{n}^{2}\leq\frac{M^{2}}{m^{2}}$ it thus suffices to show

Since $(4x+1)(1-\frac{x}{8})^{2}\geq 1$ for all $x\leq 1$ and $\frac{m}{M}\leq 1$, we conclude that $f(C_{\mathscr{H}})<0$.

Now that we have established $\max\{1,\|\Delta^{t+1}\|_{\mathscr{H}}^{2}\}\leq\max\{1,\|\Delta^{t}\|_{\mathscr{H}}^{2}\}+\frac{4M}{\widetilde{\gamma}m}\delta_{n}^{2}$, the induction step (50) follows. which completes the proof of Lemma 3.

A.4 Proof of Lemma 4

Recall that the LogitBoost algorithm is based on logistic loss $\phi(y,\theta)=\ln(1+e^{-y\theta})$, whereas the AdaBoost algorithm is based on the exponential loss $\phi(y,\theta)=\exp(-y\theta)$. We now verify the $m$-$M$-condition for these two losses with the corresponding parameters specified in Lemma 4.

The first and second derivatives are given by

It is easy to check that $|\frac{\partial\phi(y,\theta)}{\partial\theta}|$ is uniformly bounded by $B=1$.

Turning to the second derivative, recalling that $y\in\{-1,+1\}$, it is straightforward to show that

which implies that $\frac{\partial\phi(y,\theta)}{\partial\theta}$ is a $1/4$-Lipschitz function of $\theta$, i.e. with $M=1/4$.

which completes the proof for the logistic loss.

A.4.2 m𝑚m-M𝑀M-condition for AdaBoost

The AdaBoost algorithm is based on the cost function $\phi(y,\theta)=e^{-y\theta}$, which has first and second derivatives (with respect to its second argument) given by

As in the preceding argument for logistic loss, we have the bound $|y|\leq 1$ and $|\theta|\leq D$. By inspection, the absolute value of the first derivative is uniformly bounded $B:\,=e^{D}$, whereas the second derivative always lies in the interval $[m,M]$ with $M:\,=e^{D}$ and $m:\,=e^{-D}$, as claimed.