Gérard Biau

Aurélie Fischer

Benjamin Guedj
email: benjamin.guedj@inria.fr

James D Malley

COBRA: A Combined Regression Strategy

Keywords: Combining estimators, Consistency, Nonlinearity, Nonparametric regression, Prediction. 2010 Mathematics Subject Classification: 62G05, 62G20

come

Introduction

Recent years have witnessed a growing interest in combined statistical procedures, supported by a considerable research and extensive empirical evidence. Indeed, the increasing number of available estimation and prediction methods (hereafter denoted machines) in a wide range of modern statistical problems naturally suggests using some efficient strategy for combining procedures and estimators. Such an approach would be a valuable research and development tool, for example when dealing with high or infinite dimensional data.

There exists an extensive literature on linear aggregation of estimators, in a wide range of statistical models: A review of these methods may be found for example in [START_REF] Giraud | Introduction to High-Dimensional Statistics[END_REF]. Our contribution relies on a nonparametric/nonlinear approach based on an original proximity criterion to combine estimators. In that sense, it is different from existing techniques. Indeed, the present article investigates a novel point of view, motivated by the sense that nonlinear, data-dependent techniques are a source of analytic flexibility. Instead of forming a linear combination of estimators, we propose an original nonlinear method for combining the outcomes over some list of candidate procedures. We call this combined scheme a regression collective over the given basic machines. We consider the problem of building a new estimator by combining M estimators of the regression function, thereby exploiting an idea proposed in the context of supervised classification by [START_REF] Mojirsheibani | Combining classifiers via discretization[END_REF]. Given a set of preliminary estimators r 1 , . . . , r M , the idea behind this combining method is a "unanimity" concept, which is based on the values predicted by r 1 , . . . , r M for the data and for a new observation x. In a nutshell, a data point is considered to be "close" to x, and consequently, reliable for contributing to the estimation of this new observation, if all estimators predict values which are close to each other for x and this data item, i.e., not more distant than a prespecified threshold ε. The predicted value corresponding to this query point x is then set to the average of the responses of the selected observations. Let us stress here that the average is over the original outcome values of the selected observations, and not over the estimates provided by the several machines for these observations.

To make the concept clear, consider the following toy example illustrated by Figure 1. Assume we are given the observations plotted in circles, and the values predicted by two known machines r 1 and r 2 (triangles pointing up and down, respectively). The goal is to predict the response for the new point x (along the dotted line). Setting a threshold ε, the black solid circles are the data points (x i , y i ) within the two dotted intervals, i.e., such that for m = 1, 2, |r m (x i ) -r m (x)| ≤ ε. Averaging the corresponding y i 's yields the prediction for x (diamond).

We stress that the central and original idea behind our approach is that the resulting regression predictor is a nonlinear, nonparametric, data-dependent function of the basic predictors r 1 , . . . , r M , where the predictors are used to determine a local distance between a new test instance and the original (a) How should we predict the response for the query point x (dotted line)? q q q q q q q q q q q q q q (b) The two primal estimators r 1 and r 2 .

q q q q q q q q q q q q q q (c) The collective operates. q q q q q q q q q q q q q q q q q q q (d) Predicted value (diamond) for the query point x.

q q q q q q q q q q q q q q q q q q q training data. To the best of our knowledge there exists no formalized procedure in the machine learning and aggregation literature that operates as ours does. In particular, note that the original nonparametric nature of our combined estimator opens up new perspectives of research.

Indeed, though we have in mind a batch setting where the data collected consists in an n-sample of i.i.d. replications of some variable (X, Y ), our pro-cedure may be linked to other situations. For example, consider the case of functional data analysis (see [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF][START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF][START_REF] Bongiorno | Contributions in infinite-dimensional statistics and related topics[END_REF], for a survey on recent developments). Even though our method is fitted for finite dimensional data, it may be naturally extended to functional data after a suitable preprocessing of the curves. For example, this can be achieved using an expansion of the curves on an appropriate functional dictionary, and/or via a variable selection approach, as in [START_REF] Aneiros | Variable selection in infinite-dimensional problems[END_REF]. Note that in a recent work, [START_REF] Cholaquidis | An nonlinear aggregation type classifier[END_REF] adapts our procedure in a classification setting, also in a functional example.

Along with this paper, we release the software COBRA [START_REF] Guedj | COBRA: COmBined Regression Alternative[END_REF] which implements the method as an additional package to the statistical software R (see [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]. COBRA is freely downloadable on the CRAN website 1 . As detailed in Section 3, we undertook a lengthy series of numerical experiments, over which COBRA proved extremely successful. These stunning results lead us to believe that regression collectives can provide valuable insights on a wide range of prediction problems. Further, these same results

demonstrate that COBRA has remarkable speed in terms of CPU timings. In the context of high-dimensional (such as genomic) data, such velocity is critical, and in fact COBRA can natively take advantage of multi-core parallel environments.

The paper is organized as follows. In Section 2, we describe the combined estimator-the regression collective-and derive a nonasymptotic risk bound.

Next we present the main result, that is, the collective is asymptotically at least as good as any functional of the basic estimators. We also provide a rate of convergence for our procedure. Section 3 is devoted to the companion R package COBRA and presents benchmarks of its excellent performance on both simulated and real data sets, including high-dimensional models.

We also show that COBRA compares favorably with two competitors, Super Learner (van der [START_REF] Van Der Laan | Super learner[END_REF] and the exponentially weighted aggregate (see for example [START_REF] Giraud | Introduction to High-Dimensional Statistics[END_REF], in that it performs similarly in most situations, much better in some, while it is consistently faster than the Super Learner in every case. Finally, for ease of exposition, proofs and additional simulation results (figures and tables with (SM) as suffix) are postponed to a Supplementary Material.

The combined estimator

Notation

Throughout the article, we assume that we are given a training sample denoted by D n = {(X 1 , Y 1 ), . . . , (X n , Y n )}. D n is composed of i.i.d. random variables taking their values in R d ×R, and distributed as an independent prototype pair (X, Y ) satisfying EY 2 < ∞ (with the notation X = (X 1 , . . . , X d )). The space R d is equipped with the standard Euclidean metric. Our goal is to consistently estimate the regression function r

(x) = E[Y |X = x], x ∈ R d , using the data D n .
To begin with, the original data set

D n is split into two data sequences D k = {(X 1 , Y 1 ), . . . , (X k , Y k )} and D = {(X k+1 , Y k+1 ), . . . , (X n , Y n )}, with = n -k ≥ 1.
For ease of notation, the elements of D are renamed {(X 1 , Y 1 ), . . . , (X , Y )}.

There is a slight abuse of notation here, as the same letter is used for both subsets D k and D -however, this should not cause any trouble since the context is clear.

Now, suppose that we are given a collection of M ≥ 1 competing candidates r k,1 , . . . , r k,M to estimate r . These basic estimators-basic machines-are assumed to be generated using only the first subsample D k . These machines can be any among the researcher's favorite toolkit, such as linear regression, kernel smoother, SVM, Lasso, neural networks, naive Bayes, or random forests. They could equally well be any ad hoc regression rules suggested by the experimental context. The essential idea is that these basic machines can be parametric, nonparametric, or semi-parametric, with possible tuning rules. All that is asked for is that each of the r k,m (x), m = 1, . . . , M, is able to provide an estimation of r (x) on the basis of D k alone. Thus, any collection of model-based or model-free machines are allowed, and our way of combining such a collection is here called the regression collective. Let us emphasize that the number of basic machines M is considered as fixed throughout this paper. Hence, the number of machines is not expected to grow and is typically of a reasonable size (M is chosen on the order of 10 in Section 3).

Given the collection of basic machines r k = (r k,1 , . . . , r k,M ), we define the collective estimator T n to be

T n (r k (x)) = i=1 W n,i (x)Y i , x ∈ R d ,
where the random weights W n,i (x) take the form

W n,i (x) = 1 M m=1 {|r k,m (x)-r k,m (X i )|≤ε } j=1 1 M m=1 {|r k,m (x)-r k,m (X j )|≤ε } . (2.1)
In this definition, ε is some positive parameter and, by convention, 0/0 = 0.

The weighting scheme used in our regression collective is distinctive but not obvious. Starting from [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF] and [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF], we see that T n is a local averaging estimator in the following sense: The predicted value for r (x), that is, the estimated outcome at the query point x, is the unweighted average over those Y i 's such that X i is "close" to the query point. More precisely, for each X i in the sample D , "close" means that the output at the query point, generated from each basic machine, is within an ε -distance of the output generated by the same basic machine at X i . If a basic machine evaluated at X i is close to the basic machine evaluated at the query point x, then the corresponding outcome Y i is included in the average, and not otherwise. Also, as a further note of clarification: "Closeness" of the X i 's is not here to be understood in the Euclidean sense. It refers to closeness of the primal estimators outputs at the query point as compared to the outputs over all points in the training data. Training points X i that are close, in this sense, to the corresponding outputs at the query point contribute to the indicator function for the corresponding outcome Y i . This alternative approach is motivated by the fact that a major issue in learning problems consists of devising a metric that is suited to the data (see, e.g., the monograph by [START_REF] Pekalska | The Dissimilarity Representation for Pattern Recognition: Foundations and Applications[END_REF].

In this context, ε plays the role of a smoothing parameter: Put differently, in order to retain Y i , all basic estimators r k,1 , . . . , r k,M have to deliver predictions for the query point x which are in a ε -neighborhood of the predictions r k,1 (X i ), . . . , r k,M (X i ). Note that the greater ε , the more tolerant the process. It turns out that the practical performance of T n strongly relies on an appropriate choice of ε . This important question will be discussed in Section 3, where we devise an automatic (i.e., data-dependent) selection strategy of ε .

Next, we note that the subscript n in T n may be a little confusing, since T n is a weighted average of the Y i 's in D only. However, T n depends on the entire data set D n , as the rest of the data is used to set up the original machines r k,1 , . . . , r k,M . Most importantly, it should be noticed that the combined estimator T n is nonlinear with respect to the basic estimators r k,m . As such, it is inspired by the preliminary work of [START_REF] Mojirsheibani | Combining classifiers via discretization[END_REF] in the supervised classification context.

In addition, let us mention that, in the definition of the weights (2.1), all original estimators are invited to have the same, equally valued opinion on the importance of the observation X i (within the range of ε ) for the corresponding Y i to be integrated in the combination T n . However, this unanim-ity constraint may be relaxed by imposing, for example, that a fixed fraction α ∈ {1/M, 2/M, . . . , 1} of the machines agrees on the importance of X i . In that case, the weights take the more sophisticated form

W n,i (x) = 1 { M m=1 1 {|r k,m (x)-r k,m (X i )|≤ε } ≥Mα} j=1 1 { M m=1 1 {|r k,m (x)-r k,m (X j )|≤ε } ≥Mα}
.

It turns out that adding the parameter α does not change the asymptotic properties of T n , provided α → 1. Thus, to keep a sufficient degree of clarity in the mathematical statements and subsequent proofs, we have decided to consider only the case α = 1 (i.e., unanimity). Extension of the results to more general values of α is left for future work. On the other hand, as highligthed by Section 3, α has a nonnegligible impact on the performance of the combined estimator. Accordingly, we will discuss in Section 3 an automatic procedure to select this extra parameter.

Theoretical performance

This section is devoted to the study of some asymptotic and nonasymptotic properties of the combined estimator T n , whose quality will be assessed by the quadratic risk

E T n (r k (X)) -r (X) 2 .
Here and later, E denotes the expectation with respect to both X and the sample D n . Everywhere in the document, it is assumed that E|r k,m (X)| 2 < ∞ for all m = 1, . . . , M.

For any m = 1, . . . , M, let r -1 k,m denote the inverse image of machine r k,m . Assume that for any m = 1, . . . , M, r -1 k,m ((t, +∞))

t↑+∞ and r -1 k,m ((-∞, t)) t↓-∞ . (2.2)
It is stressed that this is a mild assumption which is met, for example, whenever the machines are bounded. Throughout, we let

T (r k (X)) = E [Y |r k (X)]
and note that, by the very definition of the L 2 conditional expectation,

E |T(r k (X)) -Y | 2 ≤ inf f E | f (r k (X)) -Y | 2 , (2.3)
where the infimum is taken over all square integrable functions of r k (X).

Our first result is a nonasymptotic inequality, which states that the combined estimator behaves as well as the best one in the original list, within a term measuring how far T n is from T.

Proposition 2.1. Let r k = (r k,1 , . . . , r k,M ) be the collection of basic estimators, and let T n (r k (x)) be the combined estimator. Then, for all distributions of (X, Y ) with EY 2 < ∞,

E|T n (r k (X)) -r (X)| 2 ≤ E|T n (r k (X)) -T(r k (X))| 2 + inf f E| f (r k (X)) -r (X)| 2 ,
where the infimum is taken over all square integrable functions of r k (X). In particular,

E|T n (r k (X)) -r (X)| 2 ≤ min m=1,...,M E|r k,m (X) -r (X)| 2 + E|T n (r k (X)) -T(r k (X))| 2 .
Proposition 2.1 guarantees the performance of T n with respect to the basic machines, whatever the distribution of (X, Y ) is and regardless of which initial estimator is actually the best. The term min m=1,...,M E|r k,m (X) -r (X)| 2 may be regarded as a bias term, whereas the term E|T n (r k (X)) -T(r k (X))| 2 is a variance-type term, which can be asymptotically neglected, as shown by the following result.

Proposition 2.2. Assume that ε → 0 and ε M → ∞ as → ∞. Then

E |T n (r k (X)) -T (r k (X))| 2 → 0 as → ∞,
for all distributions of (X, Y ) with EY 2 < ∞. Thus,

lim sup →∞ E T n (r k (X)) -r (X) 2 ≤ inf f E f (r k,m (X)) -r (X) 2 .
In particular,

lim sup →∞ E T n (r k (X)) -r (X) 2 ≤ min m=1,...,M E r k,m (X) -r (X) 2 .
This result is remarkable, for two reasons. Firstly, it shows that, in terms of predictive quadratic risk, the combined estimator does asymptotically at least as well as the best primitive machine. Secondly, the result is nearly universal, in the sense that it is true for all distributions of (X, Y ) such that EY 2 < ∞.

This is especially interesting because the performance of any estimation procedure eventually depends upon some model and smoothness assumptions on the observations. For example, a linear regression fit performs well if the distribution is truly linear, but may behave poorly otherwise. Similarly, the Lasso procedure is known to do a good job for non-correlated designs, with no clear guarantee however in adversarial situations. Likewise, performance of nonparametric procedures such as the k-nearest neighbor method, kernel estimators and random forests dramatically deteriorate as the ambient dimension increases, but may be significantly improved if the true underlying dimension is reasonable. Note that this phenomenon is thoroughly analyzed for the random forests algorithm in [START_REF] Biau | Analysis of a random forests model[END_REF].

The result exhibited in Proposition 2.2 holds under a minimal regularity assumption on the basic machines. However, this universality comes at a price since we have no guarantee on the rate of convergence of the variance term. Nevertheless, assuming some light additional smoothness conditions, one has the following result, which is the central statement of the paper.

Theorem 2.1. Assume that Y and the basic machines r k are bounded by some constant R. Assume moreover that there exists a constant L ≥ 0 such that, for every k ≥ 1,

|T(r k (x)) -T(r k (y))| ≤ L|r k (x) -r k (y)|, x, y ∈ R d .
Then, with the choice ε ∝ -1 M+2 , one has

E T n (r k (X)) -r (X) 2 ≤ min m=1,...,M E r k,m (X) -r (X) 2 + C -2 M+2 ,
for some positive constant C = C(R, L), independent of k.

Theorem 2.1 offers an oracle-type inequality with leading constant 1 (i.e., sharp oracle inequality), stating that the risk of the regression collective is bounded by the lowest risk among those of the basic machines, i.e., our procedure mimics the performance of the oracle over the set {r k,m : m = 1, . . . , M}, plus a remainder term of the order of -2/(M+2) which is the price to pay for combining M estimators. In our setting, it is important to observe that this term has a limited impact. As a matter of fact, since the number of basic machines M is assumed to be fixed and not too large (the implementation presented in Section 3 considers M at most 6), the remainder term is negligible compared to the standard nonparametric rate -2/(d+2) in dimension d. While the rate -2/(d+2) is affected by the curse of dimensionality when d is large, this is not the case for the term -2/(M+2) . That way, our procedure appears well armed to face high dimensional problems. When d n, many methods deteriorate and suffer from the curse of dimensionality. However, it is important to note here that even if some of the basic machines r k,1 , . . . , r k,M might be less performant in that context, this does not affect in any way our combining procedure. Indeed, forming the regression collective T n does not require any additional effort if d grows. Obviously, when d is large, the best choice would be to include as basic machines methods and models which are adapted to the high dimensional setting. This is an interesting track for future research, which is connected to functional data analysis and dimension-reduction models (see [START_REF] Goia | A partitioned single functional index model[END_REF].

Obviously, under the assumption that the distribution of (X, Y ) might be described parametrically and that one of the initial estimators is adapted to this distribution, faster rates of the order of 1/ could emerge in the bias term. Nonetheless, the regression collective is designed for much more adversarial regression problems, hence the rate exhibited in Theorem 2.1 appears satisfactory. We stress that our approach carries no assumption on the random design and mild ones over the primal estimators, in line with our attempt to design a procedure which is as model-free as possible.

The central motivation for our method is that model and smoothness assumptions are usually unverifiable, especially in modern high-dimensional and large scale data sets. To circumvent this difficulty, researchers often try many different methods and retain the one exhibiting the best empirical (e.g., cross-validated) results. Our combining strategy offers a nice alternative, in the sense that if one of the initial estimators is consistent for a given class M of distributions, then, under light smoothness assumptions, T n inherits the same property. To be more precise, assume that the initial pool of estimators includes a consistent estimator, i.e., that one of the original estimators, say r k,m 0 , satisfies

E r k,m 0 (X) -r (X) 2 → 0 as k → ∞,
for all distributions of (X, Y ) in some class M . Then, under the assumptions of Theorem 2.1, with the choice ε ∝

-1 M+2 , one has lim k, →∞ E T n (r k (X)) -r (X) 2 = 0.

Implementation and numerical studies

This section is devoted to the implementation of the described method. Its excellent performance is then assessed in a series of experiments. The companion R package COBRA (standing for COmBined Regression Alternative) is available on the CRAN website 2 , for Linux, Mac and Windows platforms (see [START_REF] Guedj | COBRA: COmBined Regression Alternative[END_REF]. Note that in a will to favor its execution speed, COBRA includes a parallel option, allowing for improved performance on multi-core computers (from [START_REF] Knaus | snowfall: Easier cluster computing (based on snow)[END_REF].

As raised in the previous section, a precise calibration of the smoothing parameter ε is crucial. Clearly, a value that is too small will discard many machines and most weights will be zero. Conversely, a large value sets all weights to 1/Σ with

Σ = j=1 1 M m=1 {|r k,m (x)-r k,m (X j )|≤ε } ,
giving the naive predictor that does not account for any new data point and predicts the mean over the sample D . We also consider a relaxed version of the unanimity constraint: Instead of requiring global agreement over the implemented machines, consider some α ∈ (0, 1] and keep observation Y i in the construction of T n if and only if at least a proportion α of the machines agrees on the importance of X i . This parameter requires some calibration.

To understand this better, consider the following toy example: On some data set, assume most machines but one have nice predictive performance. For any new data point, requiring global agreement will fail since the pool of machines is heterogeneous. In this regard, α should be seen as a measure of homogeneity: If a small value is selected, it may be an indicator that some machines perform (possibly much) better than some others. Conversely, a large value indicates that the predictive abilities of the machines are close.

A natural measure of the risk in the prediction context is the empirical quadratic loss, namely

R( Ŷ) = 1 p p j=1 ( Ŷj -Y j ) 2 ,
where Ŷ = ( Ŷ1 , . . . , Ŷp ) is the vector of predicted values for the responses Y 1 , . . . , Y p and {(X j , Y j )} p j=1 is a testing sample. We adopted the following protocol: Using a simple data-splitting device, ε and α are chosen by minimizing the empirical risk R over the set {ε ,min , . . . , ε ,max } × {1/M, . . . , 1}, where ε ,min = 10 -300 and ε ,max is proportional to the largest absolute difference between two predictions of the pool of machines.

In the package, the number #{ε ,min , . . . , ε ,max } of evaluated values may be modified by the user, otherwise the default value 200 is chosen. It is also possible to choose either a linear or a logistic scale. Figure 2 (SM) illustrates the discussion about the choice of ε and α.

By default, COBRA includes the following classical packages dealing with regression estimation and prediction. However, note that the user has the choice to modify this list to her/his own convenience:

• Lasso (R package lars, see [START_REF] Hastie | lars: Least Angle Regression, Lasso and Forward Stagewise[END_REF].

• Ridge regression (R package ridge, see [START_REF] Cule | ridge: Ridge Regression with automatic selection of the penalty parameter[END_REF].

• k-nearest neighbors (R package FNN, see [START_REF] Li | FNN: Fast Nearest Neighbor search algorithms and applications[END_REF].

• CART algorithm (R package tree, see Ripley, 2012).

• Random Forests algorithm (R package randomForest, see [START_REF] Liaw | Classification and regression by randomforest[END_REF].

First, COBRA is benchmarked on synthetic data. For each of the following eight models, two designs are considered: Uniform over (-1, 1) d (referred to as "Uncorrelated" in Table 1, Table 2 andTable 3), and Gaussian with mean 0 and covariance matrix Σ with Σ i j = 2 -|i-j| ("Correlated"). Models considered cover a wide spectrum of contemporary regression problems. Indeed, Model 1 is a toy example, Model 2 comes from van der Laan et al. ( 2007), Model 3 and Model 4 appear in [START_REF] Meier | High-dimensional additive modeling[END_REF]. Model 5 is somewhat a classic setting. Model 6 is about predicting labels, Model 7 is inspired by high-dimensional sparse regression problems. Finally, Model 8 deals with probability estimation, forming a link with nonparametric model-free approaches such as in [START_REF] Malley | Probability machines: Consistent probability estimation using nonparametric learning machines[END_REF]. In the sequel, we let N (µ, σ 2 ) denote a Gaussian random variable with mean µ and variance σ 2 . In the simulations, the training data set was usually set to 80% of the whole sample, then split into two equal parts corresponding to D k and D .

Model 1. n = 800, d = 50, Y = X 2 1 + exp(-X 2 2 ). Model 2. n = 600, d = 100, Y = X 1 X 2 + X 2 3 -X 4 X 7 + X 8 X 10 -X 2 6 + N (0, 0.5). Model 3. n = 600, d = 100, Y = -sin(2X 1 ) + X 2 2 + X 3 -exp(-X 4 ) + N (0, 0.5). Model 4. n = 600, d = 100, Y = X 1 + (2X 2 -1) 2 + sin(2πX 3 )/(2 -sin(2πX 3 )) + sin(2πX 4 ) + 2 cos(2πX 4 ) + 3 sin 2 (2πX 4 ) + 4 cos 2 (2πX 4 ) + N (0, 0.5). Model 5. n = 700, d = 20, Y = 1 {X 1 >0} +X 3 2 +1 {X 4 +X 6 -X 8 -X 9 >1+X 14 } +exp(-X 2 2 )+ N (0, 0.5). Model 6. n = 500, d = 30, Y = 10 k=1 1 {X 3 k <0} -1 {N (0,1)>1.25} . Model 7. n = 600, d = 300, Y = X 2 1 + X 2 2 X 3 exp(-|X 4 |) + X 6 -X 8 + N (0, 0.5).
4 +X 9 +sin(X 12 X 18 )+N (0,0.1)>0.38} . Table 1 presents the empirical mean quadratic error and standard deviation over 100 independent replications, for each model and design. Bold numbers identify the lowest error, i.e., the apparent best competitor. Boxplots of errors are presented in Figure 3 (SM) and Figure 4 (SM). Further, Figure 5 (SM) and Figure 6 (SM) show the predictive capacities of COBRA, and Figure 7 (SM) depicts its ability to reconstruct the functional dependence over the covariates in the context of additive regression, assessing the striking performance of our approach in a wide spectrum of statistical settings. A persistent and notable fact is that COBRA performs at least as well as the best machine, especially so in Model 3, Model 5 and Model 6.

Next, since more and more problems in contemporary statistics involve highdimensional data, we have tested the abilities of COBRA in that context. As highlighted by Table 4 (SM) and Figure 8 (SM), the main message is that COBRA is perfectly able to deal with high-dimensional data, provided that it is generated over machines, at least some of which are known to perform well in such situations (possibly at the price of a sparsity assumption). In that context, we conducted 200 independent replications for the three following models:

Model 9. n = 500, d = 1000, Y = X 1 + 3X 2
3 -2 exp(-X 5 ) + X 6 . Uncorrelated design.

Model 10. n = 500, d = 1000, Y = X 1 + 3X 2 3 -2 exp(-X 5 ) + X 6 . Correlated design.

Model 11. n = 500, d = 1500, Y = exp(-X 1 ) + exp(X 1 ) + d j=2 X j/100 j . Uncorrelated design.

A legitimate question that arises is where one should cut the initial sample D n ? In other words, for a given data set of size n, what is the optimal value for k? A naive approach is to cut the initial sample in two halfs (i.e., k = n/2): This appears to be satisfactory provided that n is large enough, which may be too much of an unrealistic assumption in numerous experimental settings. A more involved choice is to adopt a random cut scheme, where k is chosen uniformly in {1, . . . , n}. Figure 9 (SM) presents the boxplots of errors of the five default machines and COBRA with that random cutting strategy, and also shows the risk of COBRA with respect to k. To illustrate this phenomenon, we tested a thousand random cuts on the following Model 12. As showed in Figure 9 (SM), for that particular model, the best value seems to be near 3n/4.

The average risk of COBRA on a thousand replications of Model 12 is 0.3124.

Since this delivered a thousand prediction vectors, a natural idea is to take their mean or median. The risk of the mean is 0.2306, and the median has an even better risk (0.2184). Since a random cut scheme may generate some unstability, we advise practitioners to compute a few COBRA estimators, then compute the mean or median vector of their predictions.

Next, we compare COBRA to the Super Learner algorithm [START_REF] Polley | SuperLearner: Super Learner Prediction[END_REF]. This widely used algorithm was first described in van der Laan et al. ( 2007) and extended in Polley and van der Laan (2010). Super Learner is used in this section as the key competitor to our method. In a nutshell, the Super Learner trains basic machines r 1 , . . . , r M on the whole sample D n . Then, following a V -fold cross-validation procedure, Super Learner adopts a V -blocks partition of the set {1, . . . , n} and computes the matrix

H = (H i j ) 1≤ j≤M 1≤i≤n ,
where H i j is the prediction for the query point X i made by machine j trained on all remaining V -1 blocks, i.e., excluding the block containing X i . The Super Learner estimator is then

SL = M j=1 α j r j , where α ∈ arg inf α∈Λ M n i=1 |Y i -(Hα) i | 2 , with Λ M denoting the simplex Λ M = α ∈ R M : M j=1 α j = 1, α j ≥ 0 for any j = 1, . . . , M .
This convex aggregation scheme is significantly different from our collective approach. Yet, we feel close to the philosophy carried by the SuperLearner package, in that both methods allow the user to aggregate as many machines as desired, then combining them to deliver predictive outcomes. For that reason, it is reasonable to deploy Super Learner as a benchmark in our study of our collective approach.

Table 2 summarizes the performance of COBRA and SuperLearner (used with SL.randomForest, SL.ridge and SL.glmnet, for the fairness of the comparison) through the described protocol. Both methods compete on similar terms in most models, although COBRA proves much more efficient on correlated design in Model 2 and Model 4. This already remarkable result is to be stressed by the flexibility and velocity showed by COBRA. Indeed, as emphasized in Table 3 , without even using the parallel option, COBRA obtains similar or better results than SuperLearner roughly five times faster. Note also that COBRA suffers from a disadvantage: SuperLearner is built on the whole sample D n whereas COBRA only uses < n data points. Finally, observe that the algorithmic cost of computing the random weights on n test query points is × M × n test operations. In the package, those calculations are handled in C language for optimal speed performance.

Super Learner is a natural competitor on the implementation side. However, on the theoretical side, we do not assume that it should be the only benchmark. Thus, we compared COBRA to the popular exponentially weighted aggregate estimator (EWA, see [START_REF] Giraud | Introduction to High-Dimensional Statistics[END_REF]. We implemented the following version of the EWA: For all preliminary estimators r k,1 , . . . , r k,M , their empirical risks R1 , . . . , RM are computed on a subsample of D and the EWA is

EWA β : x → M j=1 ŵ j r k, j (x), x ∈ R d , where ŵ j = exp(-β R j ) M i=1 exp(-β Ri ) , j = 1, . . . , M.
The temperature parameter β > 0 is selected by minimizing the empirical risk of EWA β over a data-based grid, in the same spirit as the selection of ε and α. We conducted 200 independent replications, on Models 9 to 12. The conclusion is that COBRA outperforms the EWA estimator in some models, and delivers similar performance in others, as shown in Figure 10 (SM) and Table 5 (SM).

Finally, COBRA is used to process the following real-life data sets:

• Concrete Slump Test3 (see [START_REF] Yeh | Modeling slump flow of concrete using second-order regressions and artificial neural networks[END_REF].

• Concrete Compressive Strength4 (see [START_REF] Yeh | Modeling of strength of high performance concrete using artificial neural networks[END_REF]).

• Wine Quality5 (see [START_REF] Cortez | Modeling wine preferences by data mining from physicochemical properties[END_REF]. We point out that the Wine Quality data set involves supervised classification and leads naturally to a line of future research using COBRA over probability machines (see [START_REF] Malley | Probability machines: Consistent probability estimation using nonparametric learning machines[END_REF].

The good predictive performance of COBRA is summarized in Figure 11 (SM)

and errors are presented in Figure 12 (SM). For every data set, the sample is divided into a training set (90%) and a testing set (10%) on which the predictive performance is evaluated. Boxplots are obtained by randomly shuffling the data points a hundred times.

As a conclusion to this thorough experimental protocol, it is our belief that COBRA sets a new high standard of reference, a benchmark procedure, both in terms of performance and velocity, for prediction-oriented problems in the context of regression, including high-dimensional problems. where the infimum is taken over all square integrable functions of r k (X). In particular,

E|T n (r k (X)) -r (X)| 2 ≤ min m=1,...,M E|r k,m (X) -r (X)| 2 + E|T n (r k (X)) -T(r k (X))| 2 ,
as desired.

A.2. Proof of Proposition 2.2 Note that the second statement is an immediate consequence of the first statement and Proposition 2.1, therefore we only have to prove that

E |T n (r k (X)) -T (r k (X))| 2 → 0 as → ∞.
We start with a technical lemma, whose proof can be found in the monograph by [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF].

Lemma A.1. Let B(n, p) be a binomial random variable with parameters n ≥ 1 and p > 0. Then

E 1 1 + B(n, p) ≤ 1 p(n + 1)
and

E 1 {B(n,p)>0} B(n, p) ≤ 2 p(n + 1)
.

For all distributions of (X, Y ), using the elementary inequality (a

+ b + c) 2 ≤ 3(a 2 + b 2 + c 2 ), note that E|T n (r k (X)) -T(r k (X))| 2 = E i=1 W n,i (X) (Y i -T(r k (X i )) + T(r k (X i )) -T(r k (X)) + T(r k (X))) -T(r k (X)) 2 ≤ 3E i=1 W n,i (X)(T(r k (X i )) -T(r k (X))) 2 (A.1) + 3E i=1 W n,i (X)(Y i -T(r k (X i ))) 2 (A.2) + 3E i=1 W n,i (X) -1 T(r k (X)) 2 . (A.3)
Consequently, to prove the proposition, it suffices to establish that (A.1), (A.2) and (A.3) tend to 0 as tends to infinity. This is done, respectively, in Proposition A.1, Proposition A.2 and Proposition A.3 below.

Proposition A.1. Under the assumptions of Proposition 2.2,

lim →∞ E i=1 W n,i (X)(T(r k (X i )) -T(r k (X))) 2 = 0.
Proof of Proposition A.1. By the Cauchy-Schwarz inequality,

E i=1 W n,i (X)(T(r k (X i )) -T(r k (X))) 2 = E i=1 W n,i (X) W n,i (X) (T(r k (X i )) -T(r k (X))) 2 ≤ E j=1 W n, j (X) i=1 W n,i (X) |T(r k (X i )) -T(r k (X))| 2 = E i=1 W n,i (X) |T(r k (X i )) -T(r k (X))| 2 := A n .
The function T is such that E[T 2 (r k (X))] < ∞. Therefore, it can be approximated in an L 2 sense by a continuous function with compact support, say T (see, e.g., Theorem A.1 in [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF]. More precisely, for any η > 0, there exists a function T such that

E T(r k (X)) -T(r k (X)) 2 < η.
Consequently, we obtain

A n = E i=1 W n,i (X)|T(r k (X i )) -T(r k (X))| 2 ≤ 3E i=1 W n,i (X)|T(r k (X i )) -T(r k (X i ))| 2 + 3E i=1 W n,i (X)| T(r k (X i )) -T(r k (X))| 2 + 3E i=1 W n,i (X)| T(r k (X)) -T(r k (X))| 2 := 3A n1 + 3A n2 + 3A n3 .
Computation of A n3 . Thanks to the approximation of T by T,

A n3 = E i=1 W n,i (X)|T(r k (X)) -T(r k (X))| 2 ≤ E T(r k (X)) -T(r k (X)) 2 < η.
Computation of A n1 . Denote by µ the distribution of X. Then, du) .

A n1 = E i=1 W n,i (X)| T(r k (X i )) -T(r k (X i ))| 2 = E   1 M m=1 {|r k,m (X)-r k,m (X 1 )|≤ε } j=1 1 M m=1 {|r k,m (X)-r k,m (X j )|≤ε } | T(r k (X 1 )) -T(r k (X 1 ))| 2   . = E R d | T(r k (u)) -T(r k (u))| 2 ×E R d 1 M m=1 {|r k,m (x)-r k,m (u)|≤ε } 1 M m=1 {|r k,m (x)-r k,m (u)|≤ε } + j=2 1 M m=1 {|r k,m (x)-r k,m (X j )|≤ε } µ(dx) D k µ(
Letting

A n1 = E R d 1 M m=1 {|r k,m (x)-r k,m (u)|≤ε } 1 M m=1 {|r k,m (x)-r k,m (u)|≤ε } + j=2 1 M m=1 {|r k,m (x)-r k,m (X j )|≤ε } µ(dx) D k ,
let us prove that A n1 ≤ 2 M . To this aim, observe that

A n1 = E   R d 1 {x∈ M m=1 r -1 k,m ([r k,m (u)-ε ,r k,m (u)+ε ])} 1 + j=2 1 {X j ∈ M m=1 r -1 k,m ([r k,m (x)-ε ,r k,m (x)+ε ])} µ(dx) D k   = E   R d 1 {x∈ (a 1 ,...,a M )∈{1,2} M r -1 k,1 (I a 1 n,1 (u))∩•••∩r -1 k,M (I a M n,M (u))} 1 + j=2 1 {X j ∈ M m=1 r -1 k,m ([r k,m (x)-ε ,r k,m (x)+ε ])} µ(dx) D k   ≤ 2 M p=1 E   R d 1 {x∈R p n (u)} 1 + j=2 1 {X j ∈ M m=1 r -1 k,m ([r k,m (x)-ε ,r k,m (x)+ε ])} µ(dx) D k   .
Here,

I 1 n,m (u) = [r k,m (u) -ε , r k,m (u)], I 2 n,m (u) = [r k,m (u), r k,m (u) + ε ], and R p n (u) is the p-th set of the form r -1 k,1 (I a 1 n,1 (u)) ∩ • • • ∩ r -1 k,M (I a M
n,M (u)) assuming that they have been ordered using the lexicographic order of (a 1 , . . . , a M ).

Next, note that

x ∈ R p n (u) ⇒ R p n (u) ⊂ M m=1 r -1 k,m ([r k,m (x) -ε , r k,m (x) + ε ]).
To see this, just observe that, for all m = 1, . . . ,

M, if r k,m (z) ∈ [r k,m (u) - ε , r k,m (u)], i.e., r k,m (u)-ε ≤ r k,m (z) ≤ r k,m (u), then, as r k,m (u)-ε ≤ r k,m (x) ≤ r k,m (u), one has r k,m (x) -ε ≤ r k,m (z) ≤ r k,m (x) + ε . Similarly, if r k,m (u) ≤ r k,m (z) ≤ r k,m (u) + ε , then r k,m (u) ≤ r k,m (x) ≤ r k,m (u) + ε implies r k,m (x) - ε ≤ r k,m (z) ≤ r k,m (x) + ε . Consequently, A n1 ≤ 2 M p=1 E R d 1 {x∈R p n (u)} 1 + j=2 1 {X j ∈R p n (u)} µ(dx) D k = 2 M p=1 E µ{R p n (u)} 1 + j=2 1 {X j ∈R p n (u)} D k ≤ 2 M p=1 E µ{R p n (u)} µ{R p n (u)} D k ≤ 2 M
(by the first statement of Lemma A.1). Thus, returning to A n1 , we obtain

A n1 ≤ 2 M E T(r k (X) -T(r k (X))) 2 < 2 M η.
Computation of A n2 . For any δ > 0, write

A n2 = E i=1 W n,i (X)| T(r k (X i )) -T(r k (X))| 2 = E i=1 W n,i (X)| T(r k (X i )) -T(r k (X))| 2 1 M m=1 {|r k,m (X)-r k,m (X i )|>δ} + E i=1 W n,i (X)| T(r k (X i )) -T(r k (X))| 2 1 M m=1 {|r k,m (X)-r k,m (X i )|≤δ}
from which we get that

A n2 ≤ 4 sup u∈R d | T(r k (u))| 2 E i=1 W n,i (X)1 M m=1 {|r k,m (X)-r k,m (X i )|>δ} (A.4) + sup u,v∈R d , M m=1 {|r k,m (u)-r k,m (v)|≤δ} | T(r k (v)) -T(r k (u))| 2 . (A.5) With respect to the term (A.4), if δ > ε , then i=1 W n,i (X)1 M m=1 {|r k,m (X)-r k,m (X i )|>δ} = i=1 1 M m=1 {|r k,m (X)-r k,m (X i )|≤ε } 1 M m=1 {|r k,m (X)-r k,m (X i )|>δ} j=1 1 M m=1 {|r k,m (X)-r k,m (X j )|≤ε } = 0.
It follows that, for all δ > 0, this term converges to 0 as tends to infinity. On the other hand, letting δ → 0, we see that the term (A.5) tends to 0 as well, by uniform continuity of T. Hence, A n2 tends to 0 as tends to infinity. Letting finally η go to 0, we conclude that A n vanishes as tends to infinity.

Proposition A.2. Under the assumptions of Proposition 2.2,

lim →∞ E i=1 W n,i (X)(Y i -T(r k (X i ))) 2 = 0. Proof of Proposition A.2. E i=1 W n,i (X)(Y i -T(r k (X i ))) 2 = i=1 j=1 E[W n,i (X)W n, j (X)(Y i -T(r k (X i )))(Y j -T(r k (X j )))] = E i=1 W 2 n,i (X)|Y i -T(r k (X i ))| 2 = E i=1 W 2 n,i (X)σ 2 (r k (X i )) , where σ 2 (r k (x)) = E[|Y -T(r k (X))| 2 |r k (x)].
For any η > 0, σ 2 can be approximated in an L 1 sense by a continuous function with compact support σ2 , i.e.,

E| σ2 (r k (X)) -σ 2 (r k (X))| < η. Thus E i=1 W 2 n,i (X)σ 2 (r k (X i )) ≤ E i=1 W 2 n,i (X) σ2 (r k (X i )) + E i=1 W 2 n,i (X)|σ 2 (r k (X i )) -σ2 (r k (X i ))| ≤ sup u∈R d | σ2 (r k (u))|E i=1 W 2 n,i (X) + E i=1 W n,i (X)|σ 2 (r k (X i )) -σ2 (r k (X i ))| .
With the same argument as for A n1 , we obtain

E i=1 W n,i (X)|σ 2 (r k (X i )) -σ2 (r k (X i ))| ≤ 2 M η.
Therefore, it remains to prove that E i=1 W 2 n,i (X) → 0 as → ∞. To this aim, fix δ > 0, and note that

i=1 W 2 n,i (X) = i=1 1 M m=1 {|r k,m (X)-r k,m (X i )|≤ε } j=1 1 M m=1 {|r k,m (X)-r k,m (X j )|≤ε } 2 ≤ min    δ, 1 i=1 1 M m=1 {|r k,m (X)-r k,m (X i )|≤ε }    ≤ δ + 1 i=1 1 M m=1 {|r k,m (X)-r k,m (X i )|≤ε } >0 i=1 1 M m=1 {|r k,m (X)-r k,m (X i )|≤ε } .
To complete the proof, we have to establish that the expectation of the righthand term tends to 0. Denoting by I a bounded interval on the real line, we have

E        1    i=1 1 X i ∈ M m=1 r -1 k,m ([r k,m (X)-ε ,r k,m (X)+ε ]) >0    i=1 1 X i ∈ M m=1 r -1 k,m ([r k,m (X)-ε ,r k,m (X)+ε ])        ≤ E        1    i=1 1 X i ∈ M m=1 r -1 k,m ([r k,m (X)-ε ,r k,m (X)+ε ]) >0    1 X∈ M m=1 r -1 k,m (I) i=1 1 X i ∈ M m=1 r -1 k,m ([r k,m (X)-ε ,r k,m (X)+ε ])        + µ M m=1 r -1 k,m (I c ) = E        E        1    i=1 1 X i ∈ M m=1 r -1 k,m ([r k,m (X)-ε ,r k,m (X)+ε ]) >0    1 X∈ M m=1 r -1 k,m (I) i=1 1 X i ∈ M m=1 r -1 k,m ([r k,m (X)-ε ,r k,m (X)+ε ]) D k , X + µ M m=1 r -1 k,m (I c ) ≤ 2 ( + 1) E   1 X∈ M m=1 r -1 k,m (I) µ( M m=1 r -1 k,m ([r k,m (X) -ε , r k,m (X) + ε ]))   + µ M m=1 r -1 k,m (I c ) .
The last inequality arises from the second statement of Lemma A.1. By an appropriate choice of I, according to the technical statement (2.2), the second term on the right-hand side can be made as small as desired. Regarding the first term, there exists a finite number N of points z 1 , . . . , z N such that

M m=1 r -1 k,m (I) ⊂ ( j 1 ,..., j M )∈{1,...,N } M r -1 k,1 (I n,1 (z j 1 )) ∩ • • • ∩ r -1 k,M (I n,M (z j M )),
where I n,m (z j ) = [z j -ε /2, z j +ε /2]. Suppose, without loss of generality, that the sets Here |I| denotes the length of the interval I and x denotes the smallest integer greater than x. For all p,

r -1 k,1 (I n,1 (z j 1 )) ∩ • • • ∩ r -1 k,M (I n,M (z j M ))
x ∈ R p n ⇒ R p n ⊂ M m=1 r -1 k,m ([r k,m (x) -ε , r k,m (x) + ε ]). Indeed, if v ∈ R p n
, then, for all m = 1, . . . , M, there exists j ∈ {1, . . . ,

N } such that r k,m (v) ∈ [z j -ε /2, z j + ε /2], that is z j -ε /2 ≤ r k,m (v) ≤ z j + ε /2. Since we also have z j -ε /2 ≤ r k,m (X) ≤ z j + ε /2, we obtain r k,m (X) -ε ≤ r k,m (v) ≤ r k,m (X) + ε . In conclusion, E   1 X∈ M m=1 r -1 k,m (I) µ( M m=1 r -1 k,m ([r k,m (X) -ε , r k,m (X) + ε ]))   ≤ N M p=1 E 1 { X∈R p n } µ( M m=1 r -1 k,m ([r k,m (X) -ε , r k,m (X) + ε ])) ≤ N M p=1 E 1 { X∈R p n } µ(R p n ) = N M = |I| ε M .
The result follows from the assumption lim →∞ ε M = ∞.

Proposition A.3. Under the assumptions of Proposition 2.2,

lim →∞ E i=1 W n,i (X) -1 T(r k (X)) 2 = 0. Proof of Proposition A.3. Since | i=1 W n,i (X) -1| ≤ 1, one has i=1 W n,i (X) -1 T(r k (X)) 2 ≤ T 2 (r k (X)).
Consequently, by Lebesgue's dominated convergence theorem, to prove the proposition, it suffices to show that W n,i (X) tends to 1 almost surely. Now,

P i=1 W n,i (X) = 1 = P i=1 1 M m=1 {|r k,m (X)-r k,m (X i ))|≤ε } = 0 = P i=1 1 X i ∈ M m=1 r -1 k,m ( [r k,m (X)-ε ,r k,m (X)+ε ] ) = 0 = R d P ∀i = 1, . . . , , 1 X i ∈ M m=1 r -1 k,m ( [r k,m (x)-ε ,r k,m (x)+ε ] ) = 0 µ(dx) = R d 1 -µ(∩ M m=1 r -1 k,m [r k,m (x) -ε , r k,m (x) + ε ] ) µ(dx).
Denote by I a bounded interval. Then,

P i=1 W n,i (X) = 1 ≤ R d exp -µ(∩ M m=1 r -1 k,m [r k,m (x) -ε , r k,m (x) + ε ] ) × 1 {x∈ M m=1 r -1 k,m (I)} µ(dx) + µ M m=1 r -1 k,m (I c ) ≤ max u ue -u R d 1 {x∈ M m=1 r -1 k,m (I)} µ(∩ M m=1 r -1 k,m [r k,m (x) -ε , r k,m (x) + ε ] ) µ(dx) + µ M m=1 r -1 k,m (I c ) .
Using the same arguments as in the proof of Proposition A.2, the probability

P i=1 W n,i (X) = 1 is bounded by e -1 |I| ε M
. This bound vanishes as n tends to infinity since, by assumption, lim →∞ ε M = ∞.

A.3. Proof of Theorem 2.1 Choose x ∈ R d . An easy calculation yields that

E[|T n (r k (x)) -T(r k (x))| 2 r k (X 1 ), . . . , r k (X ), D k ] = E T n (r k (x)) -E[T n (r k (x)) r k (X 1 ), . . . , r k (X ), D k ] 2 (A.6) r k (X 1 ), . . . , r k (X ), D k + E[T n (r k (x)) r k (X 1 ), . . . , r k (X ), D k ] -T(r k (x)) 2 := E 1 + E 2 . (A.7)
On the one hand, we have

E 1 = E T n (r k (x)) -E[T n (r k (x))|r k (X 1 ), . . . , r k (X ), D k ] 2 r k (X 1 ), . . . , r k (X ), D k = E i=1 W n,i (x)(Y i -E[Y i |r k (X i )]) 2 |r k (X 1 ), . . . , r k (X ), D k .
Developing the square and noticing that E Y j |Y i , r k (X 1 ), . . . , r k (X ),

D k = E[Y j |r k (X j )],
since Y j is independent of Y i and of the X j 's with j = i, we have

E 1 = E i=1 1 M m=1 {|r k,m (x)-r k,m (X i )|≤ε } |Y i -E[Y i |r k (X i )]| 2 i=1 1 M m=1 {|r k,m (x)-r k,m (X i )|≤ε } 2 (A.8) r k (X 1 ), . . . , r k (X ), D k = i=1 V(Y i |r k (X i )) 1 M m=1 {|r k,m (x)-r k,m (X i )|≤ε } i=1 1 M m=1 {|r k,m (x)-r k,m (X i )|≤ε } 2 .
Thus,

E 1 ≤ 4R 2 1 i=1 1 M m=1 {|r k,m (x)-r k,m (X i )|≤ε } >0 i=1 1 M m=1 {|r k,m (x)-r k,m (X i )|≤ε } , (A.9)
where V(Z) denotes the variance of a random variable Z. On the other hand, recalling the notation Σ introduced in Section 3, we obtain for the second term E 2 :

E 2 = E[T n (r k (x))|r k (X 1 ), . . . , r k (X ), D k ] -T(r k (x)) 2 = i=1 W n,i (x)E[Y i |r k (X i )] -T(r k (x)) 2 1 {Σ>0} + T 2 (r k (x))1 {Σ=0} ≤ i=1 1 M m=1 {|r k,m (x)-r k,m (X i )|≤ε } |E[Y i |r k (X i )] -T(r k (x))| 2 j=1 1 M m=1 {|r k,m (x)-r k,m (X j )|≤ε } 1 {Σ>0} (A.10) + T 2 (r k (x))1 {Σ=0}
(by Jensen's inequality)

= i=1 1 M m=1 {|r k,m (x)-r k,m (X i )|≤ε } |T(r k (X i )) -T(r k (x))| 2 j=1 1 M m=1 {|r k,m (x)-r k,m (X j )|≤ε } 1 {Σ>0} (A.11) + T 2 (r k (x))1 {Σ=0} ≤ L 2 ε 2 + T 2 (r k (x))1 {Σ=0} . (A.12) Now, E|T n (r k (X)) -T(r k (X))| 2 ≤ R d E|(T n (r k (x)) -T(r k (x))| 2 µ(dx).
Then, using the decomposition (A.7) and the upper bounds (A.9) and (A.12),

E|T n (r k (X)) -T(r k (X))| 2 ≤ R d E 4R 2 1 {Σ>0} B µ(dx) + L 2 ε 2 + R d E T 2 (r k (x))1 {Σ=0} µ(dx) ≤ R d E E 4R 2 1 {Σ>0} B D k µ(dx) + L 2 ε 2 + R d E E T 2 (r k (x))1 {Σ=0} |D k µ(dx).
Thus, thanks to Lemma A.1,

E|T n (r k (X)) -T(r k (X))| 2 ≤ 8R 2 ( + 1) R d 1 µ( M m=1 {|r k,m (x) -r k,m (X)| ≤ ε }) µ(dx) + L 2 ε 2 + R d T 2 (r k (x)) 1 -µ( M m=1 {|r k,m (x) -r k,m (X)| ≤ ε }) µ(dx).
Consequently, dx) .

E|T n (r k (X)) -T(r k (X))| 2 ≤ 8R 2 ( + 1) R d 1 µ( M m=1 {|r k,m (x) -r k,m (X)| ≤ ε }) µ(dx) + L 2 ε 2 + R d T 2 (r k (x)) exp -µ( M m=1 {|r k,m (x) -r k,m (X)| ≤ ε }) µ(dx) ≤ 8R 2 ( + 1) R d 1 µ( M m=1 {|r k,m (x) -r k,m (X)| ≤ ε }) µ(dx) + L 2 ε 2 + sup x∈R d T 2 (r k (x)) max u∈R + ue -u × R d 1 µ( M m=1 {|r k,m (x) -r k,m (X)| ≤ ε }) µ(
Introducing a bounded interval I as in the proof of Proposition 2.2, we observe that the boundedness of the r k yields that

µ M m=1 r -1 k,m (I c ) = 0,
as soon as I is sufficiently large, independently of k. Then, proceeding as in the proof of Proposition 2.2, we obtain

E|T n (r k (X)) -T(r k (X))| 2 ≤ 8R 2 |I| ε M 1 + 1 + L 2 ε 2 + R 2 max u∈R + ue -u |I| ε M 1 ≤ C 1 R 2 ε M + L 2 ε 2 ,
for some positive constant C 1 , independent of k. Hence, for the choice ε ∝ -1

M+2 , we obtain

E|T n (r k (X)) -T(r k (X))| 2 ≤ C -2 M+2 ,
for some positive constant C depending on L, R and independent of k, as desired. q q q q q q q LR Lasso Ridge KNN CART RF CR 0.0 0.5 1.0 1.5

B. Numerical results

(d) Model 4.

q q q q q q q q q q q LR Lasso Ridge KNN CART RF CR 0 2 4 6 8 10 12 (e) Model 5.

q q q q q q q q q q q q q q q q q q q q q LR Lasso Ridge (f) Model 6.

q q q q q q q q q q q q q q q q q q q q q q LR Lasso Ridge KNN CART RF COBRA 0.0 0.5 1.0 1.5 2.0 2.5 3.0 (g) Model 7.

q q q q q q q q q q Lasso Ridge KNN CART RF COBRA 0.0 0.5 1.0 1.5

(h) Model 8. q q q q q q q Lasso Ridge KNN CART RF COBRA q q q q q q q q q q q q q q q q q q q q q q q q q q q LR Lasso Ridge q q q q q q q q q q q q q q q q q q LR Lasso Ridge KNN CART RF COBRA 0 20 40 60 80 (e) Model 5. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q LR Lasso Ridge KNN CART RF COBRA 0 5 10 15 20

(f) Model 6.

q q q q q q q q q q q q q q q q q q LR Lasso Ridge KNN CART RF COBRA 0 1 2 3 4 5 (g) Model 7.

q q q q q q q q q q q q q q q q q q q q q q q q q Lasso Ridge KNN CART RF COBRA 0 10 20 30 40 50 (h) Model 8. q q q q q q q q Lasso Ridge KNN CART RF COBRA (a) Model 9 q q q q q q q Lasso Ridge FNN CART RF COBRA 0 1 2 3 4 5 (b) Model 10 q q q q q q q q q q q q Lasso Ridge FNN CART RF COBRA 0 2 4 6 8 10 (c) Model 11 q q q q q q q q q q q q q q q q q Lasso Ridge FNN CART RF COBRA (a) Boxplot of errors: Initial sample is randomly cut (1000 replications of Model 12). q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Lasso Ridge FNN q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q EWA COBRA 0 q q q q q q q q q q q q q q q q q q EWA COBRA 0.00 0.05 0.10 0.15 (d) Model 12. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q EWA COBRA 0.0

Figure 1 :

 1 Figure 1: A toy example: Combining two primal estimators.

  are ordered, and denote by R p n the p-th among the N M = ( |I|/ε ) M sets.

Figure 2 (Figure 3 (

 23 Figure 2 (SM): Examples of calibration of parameters ε and α. The bold point is the minimum. (a) Model 5, uncorrelated design.

Figure 4 (

 4 Figure 4 (SM): Boxplots of quadratic errors, correlated design. From left to right: lars, ridge, fnn, tree, randomForest, COBRA. (a) Model 1.

Figure 5 (Figure 6 (

 56 Figure 5 (SM): Prediction over the testing set, uncorrelated design. The more points on the first bissectrix, the better the prediction. (a) Model 1.

Figure 8 (

 8 Figure 8 (SM): Boxplot of errors, high-dimensional models.

Figure 9 (

 9 Figure 9 (SM): How stable is COBRA?

Figure 10 (

 10 Figure 10 (SM): Boxplot of errors: EWA vs COBRA

Figure 11 (

 11 Figure 11 (SM): Prediction over the testing set, real-life data sets. (a) Concrete Slump Test.

  Figure 12 (SM): Boxplot of quadratic errors, real-life data sets. (a) Concrete Slump Test.

Table 1 :

 1 Quadratic errors of the implemented machines and COBRA. Means and standard deviations over 100 independent replications.

	Uncorr.	lars	ridge	fnn	tree	rf	COBRA
	Model 1	m. 0.1561 sd. 0.0123	0.1324 0.0094	0.1585 0.0123	0.0281 0.0043	0.0330 0.0033	0.0259 0.0036
	Model 2	m. 0.4880 sd. 0.0676	0.2462 0.0233	0.3070 0.0303	0.1746 0.0270	0.1366 0.0161	0.1645 0.0207
	Model 3	m. 0.2536 sd. 0.0271	0.5347 0.4469	1.1603 0.1227	0.4954 0.0772	0.4027 0.0558	0.2332 0.0272
	Model 4	m. 7.6056 sd. 0.9419	6.3271 10.5890 3.7358 1.0800 0.9404 0.8067	3.5262 0.3223	3.3640 0.5178
	Model 5	m. 0.2943 sd. 0.0214	0.3311 0.1012	0.5169 0.0439	0.2918 0.0279	0.2234 0.0216	0.2060 0.0210
	Model 6	m. 0.8438 sd. 0.0916	1.0303 0.4840	2.0702 0.2240	2.3476 0.2814	1.3354 0.1590	0.8345 0.1004
	Model 7	m. 1.0920 sd. 0.2265	0.5452 0.0920	0.9459 0.0833	0.3638 0.0456	0.3110 0.0325	0.3052 0.0298
	Model 8	m. 0.1308 sd. 0.0120	0.1279 0.0161	0.2243 0.0189	0.1715 0.0270	0.1236 0.0100	0.1021 0.0155
	Corr.	lars	ridge	fnn	tree	rf	COBRA
	Model 1	m. 2.3736 sd. 0.4108	1.9785 0.3538	2.0958 0.3414	0.3312 0.1285	0.5766 0.1914	0.3301 0.1239
	Model 2	m. 8.1710 sd. 1.5532	4.0071 0.6840	4.3892 0.7190	1.3609 0.4647	1.4768 0.4415	1.3612 0.4654
	Model 3	m. 6.1448 sd. 11.9450 12.0861 13.3121 11.7386 12.4160 11.1806 6.0185 8.2154 4.3175 4.0177 3.7917
	Model 4	m. 60.5795 42.2117 51.7293 9.6810 14.7731 9.6906 sd. 11.1303 9.8207 10.9351 3.9807 5.9508 3.9872
	Model 5	m. 6.2325 sd. 2.4320	7.1762 10.1254 3.1525 3.5448 3.1190 2.1468	4.2289 2.4826	2.1743 1.6640
	Model 6	m. 1.2765 sd. 0.1381	1.5307 0.9593	2.5230 0.2762	2.6185 0.3445	1.2027 0.1600	0.9925 0.1210
	Model 7	m. 20.8575 4.4367 sd. 7.1821 1.0770	5.8893 1.2226	3.6865 1.0139	2.7318 0.8945	2.9127 0.9072
	Model 8	m. 0.1366 sd. 0.0127	0.1308 0.0143	0.2267 0.0179	0.1701 0.0302	0.1226 0.0102	0.0984 0.0144

Table 3 :

 3 Average CPU-times in seconds.

	No parallelization. Means and standard
	deviations over 10 independent replica-
	tions.
	Uncorr.	SL	COBRA
	Model 1	m. 53.92 10.92 sd. 1.42 0.29
	Model 2	m. 57.96 11.90 sd. 0.95 0.31
	Model 3	m. 53.70 10.66 sd. 0.55 0.11
	Model 4	m. 55.00 11.15 sd. 0.74 0.18
	Model 5	m. 28.46 sd. 0.73	5.01 0.06
	Model 6	m. 22.97 sd. 0.27	3.99 0.05
	Model 7	m. 127.80 35.67 sd. 5.69 1.91
	Model 8	m. 32.98 sd. 1.33	6.46 0.33
	Corr.	SL	COBRA
	Model 1	m. 61.92 11.96 sd. 1.85 0.27
	Model 2	m. 70.90 14.16 sd. 2.47 0.57
	Model 3	m. 59.91 11.92 sd. 2.06 0.41
	Model 4	m. 63.58 13.11 sd. 1.21 0.34
	Model 5	m. 31.24 sd. 0.86	5.02 0.07
	Model 6	m. 24.29 sd. 0.82	4.12 0.15
	Model 7	m. 145.18 41.28 sd. 8.97 2.84
	Model 8	m. 31.31 sd. 0.73	6.24 0.11

Table 4 (

 4 SM): Quadratic errors of the implemented machines and COBRA in high-dimensional situations. Means and standard deviations over 200 independent replications.

		lars	ridge	fnn	tree	rf	COBRA
	Model 9	m. 1.5698 2.9752 3.9285 1.8646 1.5001 0.9996 sd. 0.2357 0.4171 0.5356 0.3751 0.2491 0.1733
	Model 10	m. 5.2356 5.1748 6.1395 6.1585 4.8667 2.7076 sd. 0.6885 0.7139 0.9192 0.9298 0.6634 0.3810
	Model 11	m. 0.1584 0.1055 0.1363 0.0058 0.0327 0.0049 sd. 0.0199 0.0119 0.0176 0.0010 0.0052 0.0009

Table 5 (

 5 SM): Quadratic errors of exponentially weighted aggregate (EWA) and COBRA. 200

	independent replications.
		EWA	COBRA
	Model 9	m. 1.1712 1.1360 sd. 0.2090 0.2468
	Model 10	m. 9.4789 12.4353 sd. 5.6275 9.1267
	Model 11	m. 0.0244 0.0128 sd. 0.0042 0.0237
	Model 12	m. 0.4175 0.3124 sd. 0.0513 0.0884

http://cran.r-project.org/web/packages/COBRA/index.html

http://cran.r-project.org/web/packages/COBRA/index.html

-2 exp(-X 5 ) + X 6 . Uncorrelated design.

http://archive.ics.uci.edu/ml/datasets/Concrete+Slump+Test.

http://archive.ics.uci.edu/ml/datasets/Concrete+Compressive+Strength.

http://archive.ics.uci.edu/ml/datasets/Wine+Quality.

Acknowledgements

The authors thank the Editor and two anonymous referees for providing constructive and helpful remarks, thus greatly improving the paper.

SUPPLEMENTARY MATERIAL

We have

As for the double product, notice that

Consequently,

Thus, by definition of the conditional expectation, and using the fact that T