Anas Mifrani

A Multi-objective Model Search Algorithm for Linear Regression

Keywords: Model selection, Multi-objective optimization, Ant colony optimization, Pareto front, Tradeoffs

Inherent in model selection is the problem of simultaneously optimizing multiple performance metrics. Some of these metrics express potentially conflicting criteria, like accuracy and simplicity. Pareto optimization is a branch of mathematical optimization that deals precisely with problems involving conflicting objective functions. In this article, an algorithm was developed that searches automatically for Pareto optimal linear regression models given a dataset and a set of performance metrics. The optimization task was framed as one of sequential variable selection on a graph. A search strategy was proposed that draws on ant colony optimization, a probabilistic technique well suited for graph-based problems. Experiments were run in which the metrics to be minimized were the root-mean-square error, expressing accuracy, and the number of coefficients, expressing simplicity. To substantiate the usefulness of our algorithm, cases were presented in which it outperformed AIC-based stepwise regression. Results suggested that our algorithm copes well with small datasets and correlated predictors, that it is efficient and that it informs model selection. Key properties of our algorithm were discussed and areas of improvement highlighted.

Introduction and Motivation

Part of the complexity of model selection is that it is usually motivated by multiple and potentially competing objectives. A typical example is when a model has to conform as much as possible to available data, but also has to be sufficiently simple for it to generalize and be interpretable [START_REF] Breiman | Statistical modeling: The two cultures[END_REF][START_REF] Gregorutti | Correlation and variable importance in random forests[END_REF]. Some sort of compromise between accuracy and simplicity is often necessary because complex models tend to approximate data better than simple ones. The same observation seems warranted for any pair of conflicting criteria for model assessment.

Optimal model selection has been studied primarily within the automated machine learning community [START_REF]Automatic Machine Learning: Methods, Systems, Challenges[END_REF][START_REF] He | Automl: A survey of the state-of-the-art[END_REF][START_REF] Shubhra | Automl to date and beyond: Challenges and opportunities[END_REF], with a common point of departure being that model selection can be framed as a sequential decision making problem. Decisions include the choice of algorithms (e.g., linear regression), methods (e.g., imputation of missing data) and hyperparameters [START_REF] Thornton | Auto-weka: Combined selection and hyperparameter optimization of classification algorithms[END_REF]. The goal is to find sequences, and therefore models, that minimize a loss function. Hutter et al. [START_REF]Automatic Machine Learning: Methods, Systems, Challenges[END_REF] point out the difficulty in applying classical optimization to such a problem, citing as their main reasons (a) the complexity of the search space and (b) the fact that little is known about the loss's analytical properties, say convexity, when viewed as a function of the decisions. Accordingly, Feurer et al. [START_REF] Feurer | Efficient and robust automated machine learning[END_REF] explore the search space using Bayesian optimization (BO) [START_REF] Mockus | The application of Bayesian methods for seeking the extremum[END_REF], which involves repetitive sampling of sequences around the current best sequence, based on performance estimates provided by a surrogate.

Rokatoarison et al. [START_REF] Herilalaina Rakotoarison | Automated machine learning with Monte-Carlo Tree Search[END_REF] use Monte Carlo tree search [START_REF] Maurice | Monte-carlo tree search[END_REF] for structures, i.e algorithms and methods, and BO for hyperparameters. Structures are built incrementally by "walking down" a tree, adding one algorithm or method at a time, whereas hyperparameters are sampled as previously. Other works have focused on particular areas of model selection.

For instance, Khurana et al. [START_REF] Khurana | Cognito: Automated feature engineering for supervised learning[END_REF][START_REF] Khurana | Feature engineering for predictive modeling using reinforcement learning[END_REF] propose automatic feature engineering algorithms for regression and classification. Their goal is to determine, based on a predefined set of functions, what sequence of variable transformations yields the most predictively accurate model when applied on a dataset. For this they define a "transformation tree" in which nodes represent transformed datasets and edges represent functions, then they compare ad hoc search heuristics [START_REF] Khurana | Cognito: Automated feature engineering for supervised learning[END_REF] with strategies obtained through reinforcement learning [START_REF] Khurana | Feature engineering for predictive modeling using reinforcement learning[END_REF]. a convolutional neural network's hyperparameters. Their method involves training a recurrent neural network to select hyperparameters in such a way as to maximize expected model accuracy.

A drawback to the general approach represented by the aforementioned works is that it is single-objective, whereas model selection is essentially multi-objective. Here we propose a different approach, in which more than one metric is considered. We assume that the metrics conflict; that is, there exists no solution that optimizes all of them simultaneously, in which case the expert would be prepared to make tradeoffs. Thus we are primarily interested in identifying the set of models that cannot be improved in any metric without causing the others to deteriorate; this is an adaptation of a core Pareto (or multi-objective) optimization concept known as the Pareto front [START_REF] Cao | A hybrid elitist Paretobased coordinate exchange algorithm for constructing multi-criteria optimal experimental designs[END_REF][START_REF] Wei | The fairness-accuracy Pareto front[END_REF][START_REF] Bierbrauer | The Pareto-frontier in a simple Mirrleesian model of income taxation[END_REF].

Albeit equivalent, models which lie on the Pareto front can be seen as representing different tradeoffs between the objectives. From a decision making perspective, therefore, the principal appeal of our approach is that it focuses the expert's attention on a set of equally good alternatives that they can ultimately choose from based on their preferences and the value placed on each objective [START_REF] Chapman | Impact of response variability on Pareto front optimization[END_REF][START_REF] Lu | Optimal designed experiments using a Pareto front search for focused preference of multiple objectives[END_REF]. Another appeal is that by examining the shape of the Pareto front, the expert is able to gauge how much compromise is needed from some metric to improve others [START_REF] Cao | A hybrid elitist Paretobased coordinate exchange algorithm for constructing multi-criteria optimal experimental designs[END_REF], and is hence able to locate regions in the model space where small concessions in some metrics permit considerable improvements in others. From a statistical perspective, our approach may assist in improving our understanding of the relationships among metrics and in informing theoretical investigations in that direction.

In this article, we present a model selection algorithm which searches for Pareto optimal models of a continuous response variable, given a dataset and a set of performance metrics. It does this by efficiently exploring a model search graph. We devise and implement a search heuristic which draws on ant colony optimization (ACO), a combinatorial optimization technique modeled after ant colonies [START_REF] Dorigo | Ant Colony Optimization[END_REF]. ACO's first application was the traveling salesperson problem [START_REF] Dorigo | Ant system: Optimization by a colony of cooperating agents. IEEE transactions on systems, man, and cybernet-ics[END_REF], but its scope has expanded tremendously in the last two decades, with applications in experiment design [START_REF] Borrotti | Naïve bayes ant colony optimization for designing high dimensional experiments[END_REF], project scheduling [START_REF] Zhu | A discrete oppositional multi-verse optimization algorithm for multi-skill resource constrained project scheduling problem[END_REF], spatial clustering [START_REF] İnkaya | Ant colony optimization based clustering methodology[END_REF], developed for Pareto optimization [START_REF] García-Martínez | A taxonomy and an empirical analysis of multiple objective ant colony optimization algorithms for the bi-criteria traveling salesperson problem[END_REF][START_REF] Angus | Multiple objective ant colony optimisation[END_REF], one of which was used in our work.

The remainder of the article is organized as follows. Section 2 states the combinatorial optimization problem underpinning the article, provides some terminology and introduces the model search graph. Section 3 describes a generic ACO algorithm and shows how it can be adapted to Pareto optimization. Section 4 presents our algorithm. On the basis of experiments, Section 5 compares the algorithm with AIC-based stepwise regression [START_REF] Kabe | Stepwise multivariate linear regression[END_REF][START_REF] Akaike | Information Theory and an Extension of the Maximum Likelihood Principle[END_REF], and sheds light on some of its basic properties.

Model Search Graph

We concentrate on linear regression models of a response Y that may feature powers of variables X 1 , ..., X p , p ≥ 1, together with two-way interactions. For p = 3 and a maximum exponent e max = 2, these models have as an upper bound

Y = β 0 +β 1 X 1 +β 2 X 2 +β 3 X 3 +β 4 X 2 1 +β 5 X 2 2 +β 6 X 2 3 +β 7 X 1 X 2 +β 8 X 1 X 3 +β 9 X 2 X 3 +ϵ (1)
This work is an attempt to construct an efficient procedure for building satisfactory models of this form, in a way similar to forward selection but one which is stochastic and multiobjective. We first show that building such models is equivalent to traversing a graph.

Models as Paths on a Graph

Given p and e max , M (p, e max ) will denote the set of models described above. This set is finite. Let again p = 3 and e max = 2. There are 8 ways of choosing two-way interactions:

(1)

X 1 X 2 ; (2) X 1 X 3 ; (3) X 2 X 3 ; (4) X 1 X 2 and X 1 X 3 ; (5) X 1 X 2 and X 2 X 3 ; (6) X 1 X 3 and X 2 X 3 ; (7) X 1 X 2 , X 2 X 3 and
= β 0 + β 1 A + β 2 B + +β 3 A 2 + β 4 AB + ϵ.
and e max by following the same construction principles as above; that is, by enumerating the possible combinations of interaction terms, of power-of-one terms, and so on. Some comments on the graph are supplied in the online Appendix C.

Problem Statement

Let F 1 , ..., F Q be Q performance metrics to be minimized with respect to a dataset D. For fixed p and e max , a combinatorial Pareto optimization problem follows:

min x ∈ M (p, emax) (F 1 (x, D), ..., F Q (x, D)) (2)
where F q (x, D), q ∈ {1, . . . , Q}, denotes the value of F q for a model x ∈ M (p, e max ) with respect to D.

In terms of the model search graph, and assuming as in the Introduction that the metrics conflict, solving (2) amounts to finding the full paths corresponding to models that are Pareto optimal with respect to F 1 , . . . , F Q (Definition 1 and 2). The Pareto front [START_REF] Lu | Optimal designed experiments using a Pareto front search for focused preference of multiple objectives[END_REF] is the set of such models.

Definition 1 Let x and y be models in M (p, e max ). x is said to dominate y on D if F q (x, D) ≤ F q (y, D) for all q ∈ {1, . . . , Q} and F k (x, D) < F k (y, D) for some k ∈ {1, . . . , Q}.

Definition 2 Let x in M (p, e max ). x is Pareto optimal (or nondominated, or efficient) if, for all y ∈ M (p, e max ), y does not dominate x.

We would like to point out, in passing, that the definitions above make no reference to the character of the relationships among the metrics. Indeed, even with nonconflicting metrics, one can still talk of a Pareto front, namely the singleton containing the model that minimizes all the metrics at once. In this case, however, the front's interest is clearly reduced, and Problem (2) could in principle be solved using single-objective optimization techniques.

Computing the Pareto front may be intractable [START_REF] Chapman | Impact of response variability on Pareto front optimization[END_REF]. Only an exact algorithm could determine with certainty whether a model is Pareto optimal, but such an algorithm typically proceeds by implicit enumeration of the search space [START_REF] Bierlaire | Optimization: Principles and Algorithms[END_REF], which makes it unsuitable for large graphs. This motivates the interest in heuristics, which are computationally efficient techniques for generating good solutions [START_REF] García-Ródenas | A comparison of general-purpose optimization algorithms for finding optimal approximate experimental designs[END_REF]; that is, solutions that are "feasible and significantly better than a solution that would have been designed by a human expert" [START_REF] Bierlaire | Optimization: Principles and Algorithms[END_REF].

Constructive heuristics, meaning those that build solutions stepwise, are known to be the fastest [START_REF] Dorigo | Ant Colony Optimization[END_REF][START_REF] Lai | Self-adaptive randomized constructive heuristics for the multi-item capacitated lot sizing problem[END_REF], but they also seem to be the most relevant to the objective we outlined in introducing this section. Even more relevant are algorithms based on ant colony optimization [START_REF] Dorigo | Ant Colony Optimization[END_REF], which are designed specifically for problems that can be reduced to finding desirable paths in a graph.

It is important to note that most applications of Pareto optimization aim not at computing the Pareto front but rather at approximating it [START_REF] Sambo | A coordinate-exchange two-phase local search algorithm for the d-and i-optimal design of split-plot experiments[END_REF]. Typically, a good Pareto front approximation will be a set of solutions in which (a) no solution dominates another; (b) the solutions are close to those on the true front; and (c) the solutions are diverse [START_REF] Wei | The fairness-accuracy Pareto front[END_REF][START_REF] Cao | On using the hypervolume indicator to compare Pareto fronts: Applications to multi-criteria optimal experimental design[END_REF][START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF].

Perhaps the best known measure of the quality of a Pareto front approximation that accounts for (b) and (c) is the hypervolume, which is the size of the objective space enclosed by the points on the front and a dominated reference point [START_REF] Zitzler | Multiobjective optimization using evolutionary algorithms-a comparative case study[END_REF]. The hypervolume has the property of being strictly Pareto compliant [START_REF] Cao | On using the hypervolume indicator to compare Pareto fronts: Applications to multi-criteria optimal experimental design[END_REF]; that is, if A and B are Pareto front approximations such that the points in B are dominated by those in A, then A has a higher hypervolume than B. This observation will be useful in Section 5, where we compare different algorithm settings based on the Pareto front approximations obtained in each one.

Multi-objective Ant Colony Optimization

We now turn to ant colony optimization, an umbrella term for a group of algorithms in which artificial agents imitate the behavior of foraging ants to solve combinatorial optimization problems [START_REF] Dorigo | Ant Colony Optimization[END_REF]. The basis for ACO is the assumption that simulation of stigmergy,

i.e the discharging and following of pheromones by insects, can solve problems analogous to finding the shortest path in a graph. ACO deals traditionally with problems involving a function to be minimized with respect to a set of discrete decision variables subject to certain feasibility constraints. We will see, however, that the underlying ideas can easily be extended to the case where multiple functions are to be minimized.

ACO ants construct solutions incrementally and independently by traversing a graph whose edges represent value assignments for decision variables. An edge (i, j) is described at time t by the amount of pheromones present on it, τ ij (t), and a so-called "heuristic value" (HV), η ij [START_REF] Dorigo | Ant Colony Optimization[END_REF]. Pheromones are "discharged" during travel, with paths closer to optimal solutions receiving greater amounts of pheromones. Hence, τ ij (t) is a cumulative indicator of the quality of solutions which had been built up to time t using (i, j). η ij expresses an a priori belief about whether (i, j) belongs to an optimal solution. When approaching the TSP, for example, it is customary to let η ij = 1 d ij , d ij being the distance between cities-nodes i and j [START_REF] Branke | Solving the probabilistic TSP with ant colony optimization[END_REF]; the intuition is that a short tour can be obtained by repeatedly going to the closest yet-to-be-visited city until no city is left.

An ant chooses its destinations as a function of the pheromones and HVs of the edges available. Regions with relatively high concentrations of pheromones tend to be preferred [START_REF] Borrotti | Naïve bayes ant colony optimization for designing high dimensional experiments[END_REF]. Once a complete solution has been built, it is evaluated and pheromones are deposited on its constituent edges. The amount of pheromones deposited is commensurate with the quality of the solution. Additionally, pheromones "evaporate" with time. Ideally this prevents a too rapid convergence to suboptimal regions, and encourages exploration [START_REF] Dorigo | Ant Colony Optimization[END_REF][START_REF] Wan | Acomcd: A multiple cluster detection algorithm based on the spatial scan statistic and ant colony optimization[END_REF].

A multi-objective version of ACO was proposed by [START_REF] García-Martínez | A taxonomy and an empirical analysis of multiple objective ant colony optimization algorithms for the bi-criteria traveling salesperson problem[END_REF] in which τ ij (t) and η ij are Qsized vectors, where Q is the number of objectives, such that τ ij (t) = (τ ij,1 (t), ..., τ ij,Q (t)) and η ij = (η ij,1 , ..., η ij,Q ). Ants evaluate their solutions and deposit a separate amount of pheromones for each objective. For all (i, j) and q ∈ {1, . . . , Q}, the q-th component of

τ ij (t) is updated as τ ij,q (t + 1) = (1 -ρ) • τ ij,q (t) + m k=1 ∆τ k ij,q (3)
where ρ ∈ [0, 1] is the evaporation rate; m the number of ants in the colony; ∆τ k ij,q the amount of pheromones deposited by ant k with respect to F q , which measures solution quality with respect to F q . The Q updates occur simultaneously.

For concreteness, suppose Q = 2. Assuming ant k occupies node i in the construction graph, it will move to a successor node j with a probability

p k ij =      (τ λ 1 ij,1 •τ λ 2 ij,2 ) α (η λ 1 ij,1 •η λ 2 ij,2 ) β p∈S(i) (τ λ 1 ip,1 •τ λ 2 ip,2 ) α (η λ 1 ip,1 •η λ 2 ip,2 ) β (i, j) ∈ N (s k ) 0 (i, j) / ∈ N (s k ) (4)
where s k denotes the path that k has traveled so far; N (s k ) the set of edges that can be added to s k without violating the problem's constraints; S(i) the successors of i; (α, β)

parameters controlling the importance of an edge's pheromones relative to its HVs; and where λ q ∈ [0, 1] is the weight of F q relative to the other objective, such that λ 1 + λ 2 = 1.

Time was dropped from the expression for brevity.

Proposed Algorithm

We now describe an algorithm for solving Problem (2). We will call it MOMSACO, which is shorthand for Multi-objective Model Selection based on Ant Colony Optimization. We treat the model search graph as a construction graph for a colony of size m. For simplicity, a Pareto front approximation will henceforth be referred to as a Pareto front. To avoid ambiguity, the Pareto front will be referred to as the true Pareto front. We use the notation of Section 2 and 3.

Initialization

The algorithm requires a dataset in which the columns representing Y and X 1 , ..., X p were designated. The dataset is randomly partitioned into a training dataset, D train , and a test dataset, D test . A value for e max is required for generating the model search graph.

HVs are initialized after the graph has been generated. Informally, η Fq ij should be high if selection of j is likely to bring an immediate improvement in F q . In the TSP, for example, the biggest immediate improvement in traveled distance comes from visiting the closest city.

It is hard to draw an analogy here, because model assessment metrics are not necessarily separable in the way distance is. Instead, if the current node is i, we infer the likely improvement in F q resulting from the selection of successor node j based on a sample of models containing both i and j. Specifically, we let, for q ∈ {1, . . . , Q},

η Fq ij = Φ Fq i F q (., D test ) ij (5)
where:

-F q (., D test ) ij is the mean of F q , with respect to D test , in a random sample of n init models containing (i, j). n init ≥ 1 is a user-defined parameter.

-Φ Fq i is the maximum of F q (., D test ) il over l ∈ S(i).

Model Search

The m ants construct models independently by walking from Start to End, selecting at each stage a node from the distribution defined by Equation [START_REF] Thornton | Auto-weka: Combined selection and hyperparameter optimization of classification algorithms[END_REF]. Let us describe this process in detail for some ant k, k ∈ {1, . . . , m}. The ant is located initially in Start. A successor node is chosen after calculating p k ij , where i = Start, for all successors j. The ant joins this successor then continues in the same manner until it reaches End. Note that if an interaction node other than N one was selected, the ant is required to choose N one for its next moves until the next node is a power-of-one node (for an explanation, see the online Appendix C). Upon reaching End, a model will have been fully constructed; call it m k . Its coefficients are estimated on D train then its metrics are evaluated on D test . The ant subsequently returns to Start by taking the opposite path. Along the way, it deposits a separate amount of pheromones for each metric, ∆τ k ij,q (q ∈ {1, . . . , Q}), on every visited edge. These are identical for all the edges taken. A percentage ρ of the old pheromones evaporate. These steps occur at the same time colony-wide, and constitute a full iteration.

Conceptually, ∆τ k ij,q should be a nonincreasing function of F q . Hence, we take ∆τ k ij,q to be the inverse of the value obtained for F q by ant k, i.e

∆τ k ij,q = δ ij,k F q (m k , D test ) ( 6
)
where δ ij,k = 1 if (i, j) was visited by k, and 0 otherwise. At the end of an iteration, τ ij,q

is updated for all q ∈ {1, . . . , Q} according to

τ ij,q (t + 1) = (1 -ρ) • τ ij,q (t) + m k=1 δ ij,k F q (m k , D test ) (7)

Populating a Pareto Front

A Pareto front is formed as iterations pass. To form the initial front, we compare the first iteration's models, then we drop the dominated ones. The remaining models could be Pareto optimal, so they constitute the initial front. From the second iteration onward, new models are compared with each other then with the present front. If a new model is not dominated by any of the others with which it was compared, then it is added to the front.

However, if a model on the front turns out to be dominated by some new model, then it is removed. This procedure ensures that out of all the models found up to the present, only those that have not been proven to be dominated compose the front.

Results and Analysis

We begin this section with an experimental comparison of AIC-based stepwise regression (henceforth AIC-SR) [START_REF] Kabe | Stepwise multivariate linear regression[END_REF] and MOMSACO, the two being sequential, multi-objective procedures for model selection (Section 5.1). Section 5.2 through 5.4 deal with some salient properties of MOMSACO.

In all the experiments reported, the metrics to be minimized were the root-mean-square error (RMSE) and the number of coefficients (NoCoeff). RMSE indicates accuracy because it captures the differences between a model's predictions and the values observed for the response [START_REF] Gilmore Pontius | Components of information for multiple resolution comparison between maps that share a real variable[END_REF]. NoCoeff describes a model's simplicity. The relevance of a Pareto front-based approach stems from the observation that it is usually impossible to minimize RMSE and NoCoeff simultaneously: when NoCoeff is minimal, i.e when there are no variables in the model, RMSE is generally suboptimal. One exception is when the response is constant, in which case the intercept-only model would score an RMSE of zero.

To initialize MOMSACO, eight models per edge were drawn uniformly (n init = 8), and Equation ( 7) was used. Initial pheromones were of one per metric per edge. Experiments spanned 100 iterations, and were run 50 times due to the randomness of navigation. For each dataset, 70% of the observations were used for training, and 30% for calculating the metrics. There were 15 ants. In Section 5.1, ρ was set at 0.5, λ 1 and λ 2 at 0.5, α and β at 1. A Pareto front's hypervolume was the total area of the rectangles enclosed by the points on the front and a reference point. This point's location varied with each dataset.

Details on how the reference points and hypervolumes were determined are supplied in the online Appendix A. In Section 5.4 and 5.5, the term "transition probability" refers to the quantity defined by Equation [START_REF] Thornton | Auto-weka: Combined selection and hyperparameter optimization of classification algorithms[END_REF]. The experiments were run in R, on a laptop with an Intel i3 processor running at 2.10GHz and using 4GB of RAM.

Examples Where AIC-SR Fails But MOMSACO Does Not

The motivation for comparing MOMSACO with AIC-SR is threefold: (1) they are both multi-objective procedures; (2) they differ in fundamental ways; and (3) if cases exist in which MOMSACO outperforms AIC-SR, then positive conclusions can be drawn about its usefulness. Perhaps the most conspicuous difference between the two lies in the nature of the procedure. AIC-SR is deterministic, subject to a fixed set of rules based on statistical tests;

MOMSACO is stochastic, iterative and one that "learns" over time. Secondly, consider the way in which candidate models are compared. AIC-SR compares models on the basis of AIC, with lower values preferred [START_REF] Akaike | Information Theory and an Extension of the Maximum Likelihood Principle[END_REF]. In contrast, MOMSACO does not score models; rather, it compares their metrics then discards those that are dominated. The former mode of comparison produces a ranking of models; the latter supplies a range of models, each of which represents a particular compromise among the metrics being studied. The latter also allows some awareness of the relationships among the metrics; in particular, it helps gauge the extent to which they conflict, and identify opportunities for substantial improvements in one metric that require small concessions in another. Awareness of such relationships, which is important to decision making, is precluded by the first mode of comparison.

We now give examples of synthetic datasets where (a) MOMSACO outperforms AIC-SR in the sense of being able to find data-generating models, and (b) comparing models according to Pareto dominance is preferable to ranking them according to AIC. Such datasets abound, but for the sake of exposition we focus on two. Both were comprised of variables X 1 , X 2 , X 3 and X 4 and a response Y. In D 1 , the four variables were drawn independently from normal distributions; in D 2 , they were drawn from a multivariate normal distribution with X 1 and X 2 correlated. Y was generated in each case according to

Y = 2 + X 1 + 1.5X 2 + X 1 X 2 + X 2 4 + ϵ (8)
where ϵ was drawn from the standard normal distribution. We will henceforth refer to [START_REF] Maurice | Monte-carlo tree search[END_REF] as the true model. X 1 , X 2 , X 3 and X 4 were all allowed to feature in the models, so p was equal to four. e max was taken to be two. This produced a search graph with 16 384 models.

Details on the simulation of D 1 and D 2 (and on a third, unreported example) are provided in the online Appendix B.

To substantiate (a) and (b), we generated N = 1000 observations in D 1 and 10 000 in D 2 , then we ran two experiments per dataset. In the first, we computed the AIC value of every model in the graph, then we ran a full search of M (4, 2) to find the true Pareto front with respect to RMSE and NoCoeff. In the second, we ran forward, backward and bidirectional AIC-SR, then we initialized MOMSACO and ran 50 trials of it. An aggregate Pareto front was then obtained by combining the fronts of the MOMSACO trials and discarding the dominated models. Figure 3 shows both aggregate fronts. To ensure a fair comparison, we required AIC-SR's scope to be the same as MOMSACO's, i.e Table 1: Effect of sample size and correlation on the ability of both algorithms to find the true models. Y denotes success in finding these models, and N denotes failure.

{X 1 , X 2 , X 3 , X 4 , X 2 1 , . . . , X 2 4 , X 1 X 2 , . . . , X 3 X 4 }.
Sample size 1000 10 000 100 000

D 1 (MOMSACO) Y Y Y D 1 (AIC-SR) N Y Y D 2 (MOMSACO) Y Y Y D 2 (AIC-SR) N N N
In both datasets, the true model was truly Pareto optimal with respect to RMSE and

NoCoeff, but suboptimal with respect to AIC. Several models had better AIC values than the true one. This suggests that Pareto dominance is a sound basis for comparing models when RMSE and NoCoeff are its metrics. It also shows that there are situations in which a ranking of models based on AIC can be misleading in that it places data-generating models below others. We also found that whereas AIC-SR did not succeed in finding the true model, MOMSACO did. This was true for both datasets (see the triangles in Figure 3), which proves that MOMSACO can outperform AIC-SR.

Illustrations of how Pareto dominance can inform model selection are provided by both fronts' bottom two models (Figure 3), which show that a small concession in NoCoeff (from two variables to three) allows for significant gains in RMSE (from approximately 34.4 to 1.8 in D 1 ). Both Pareto fronts also give clear insight into the tradeoff represented by each model, an insight not provided by AIC values.

Additional experiments with D 1 and D 2 showed that the relative performance of AIC-SR was affected by the sample size and the presence of correlated predictors (Table 1). For N = 10 000, AIC-SR was finally able to retrieve the true model in D 1 . This was not observed in D 2 , which differs with D 1 only in that it was generated by a model with correlated predictors. This suggests the following: (a) in the absence of correlated predictors, AIC-SR requires a larger sample than MOMSACO to find the data-generating model; (b) when correlations exist, AIC-SR may require a much larger sample, which however was not confirmed by the experiments, because AIC-SR failed even with N = 10 6 . In any event, it seems that in the presence of correlations, MOMSACO can recover the data-generating model regardless of the sample size.

Pareto Fronts

Now we concentrate solely on MOMSACO. In each of the previous experiments, the Pareto front's size and hypervolume kept increasing until they leveled out. Figure 4 depicts this process for D 1 (N = 1000) for various evaporation rates, omitting the fronts' sizes for clarity. In general, the hypervolume improved most rapidly within the first few iterations.

This was when the fronts were being populated at the fastest rate (Table 2).

While work is needed to pinpoint the effect of each parameter on the Pareto front, experiments suggested that evaporation is an important factor. Notably, the hypervolume This means that for smaller values of ρ, the discovery of new nondominated models spanned a longer period of time, and higher quality fronts were formed.

It is noteworthy that, on average, MOMSACO was able to recover 26.54% of the true Pareto front (see the rates of contribution in Table 3) by evaluating as few as 219 models,

i.e 1.33% of the total size of the search graph.

Ant Navigation

Figure 4 describes the evolution of Pareto front quality over time, using the hypervolume as an indicator. It may be interesting, from an optimization point of view, to ask how the mean of each metric evolves while the front is being populated. Figure 5 shows that as the search was proceeding, there was a gradual decline, and therefore improvement, in the average values of RMSE and NoCoeff among models on the Pareto front. This indicates that both metrics were being optimized.

The improvement observed in Figure 4 and 5 can be attributed to the movement of the pheromones. Initially, these do not contribute to navigation because they are equally distributed among the edges. Rather, it is the HVs, which differ from an edge to another, that determine the initial transition probabilities. Therefore, the ants will most likely start by taking the paths with the highest HVs. Later, the pheromones begin to vary as a result of their deposition and evaporation. This affects navigation by altering the probabilities of the traversed edges, so that those leading to lower metrics become more probable.

On Initialization

Since the HVs determine the initial transition probabilities, they significantly alter the behavior of the colony in the first iteration. MOMSACO's performance in subsequent iterations is also, by that very fact, dependent on the HVs. The method by which these are computed is a critical component of ACO algorithms in general, and of MOMSACO in particular.

To evaluate the initialization scheme described in Section 4.1, further experiments were run (1) with an arbitrary assignment of HVs (50 trials), then [START_REF] Gregorutti | Correlation and variable importance in random forests[END_REF] with no HVs at all, i.e β = 0 in Equation ( 6) (50 trials). Figure 6 juxtaposes the mean hypervolumes of both sets of experiments with that of the original D 1 experiment (N = 1000). We can see that in the original experiment, i.e when the HVs had been set using Equation ( 7), the Pareto front was of a consistently superior quality, and was thus closer to the true front. When the HVs were initialized arbitrarily, the resulting Pareto front was the poorest. Two conclusions

Conclusion

A model selection algorithm was developed that searches automatically for Pareto optimal linear regression models, given a dataset and a set of performance metrics. Candidate models combined two-way interactions with powers of the original explanatory variables.

The task of finding Pareto optimal models was framed as that of finding the corresponding paths on the model search graph. A search strategy was proposed that draws on ACO, a probabilistic technique well suited for problems involving graphs.

Experiments with simulated data showed that MOMSACO could recover the underlying models where AIC-SR could not. They also showed that these were actually suboptimal with respect to AIC, but Pareto optimal with respect to RMSE and NoCoeff. Only by increasing the sample size or reducing the correlations in the data-generating model was AIC-SR able to find it. These results suggest that MOMSACO is better suited for correlated predictors and small datasets than AIC-SR, and that comparing models according to Pareto dominance can inform model selection.

By juxtaposing the Pareto front's rate of contribution with the total number of models explored, we were able to conclude that MOMSACO can efficiently navigate the model search graph. We also found our initialization method to be a desirable component of MOMSACO because it enabled the formation of a higher quality Pareto front.

MOMSACO can be used as a variable selection procedure for potentially any type of regression model. Because it does not presuppose the nature of the models it is searching for, and because these are trained independently of it, the procedure can apply equally to non-linear models. If hyperparameters are required, however, then grid search or Bayesian optimization can supplement MOMSACO.

A number of questions remain. One is how the parameters should be set, given some reference point for hypervolumes, in order to obtain the best possible Pareto fronts. This is in itself an optimization problem, and a complex one, primarily because the data is a factor in the optimality of a configuration of parameters. We intend to study the matter in detail, hoping to enunciate some guidelines on how the parameters should be chosen in broad enough situations. A second question is how to better the representation of M (p, e max ) so that computational issues associated with large p and e max could be alleviated. For large p, the model search graph may not fit into the memory of a computer. Immediate solutions are to consider only interactions and power-of-one terms, or to subject the variables to preliminary dimensionality reduction via techniques like principal component analysis [START_REF] Johnstone | On consistency and sparsity for principal components analysis in high dimensions[END_REF].

A third question is how the current initialization procedure could be improved. We think incorporating some form of past knowledge may be beneficial. If, for instance, datasets similar to the one being used are available, and models that performed well on them are known, then information about these can be leveraged -for example, by awarding relatively high values to the constituent edges of similar models in the model search graph.

A database could be set up that stores key statistical characteristics of a large number of datasets and corresponding models, which characteristics could be used to measure the similarity between a dataset of interest and those in the database, and between the graph models and the database models. Datasets and models could be extracted from academic papers, open machine learning competitions, experts and other relevant sources.

SUPPLEMENTARY MATERIAL

Online Appendix A supplies an analytic expression for the hypervolume of a twodimensional Pareto front and explains how reference points were chosen in Section

Figure 1 :

 1 Figure 1: Model search graph for p = 3 and e max = 2, with A = X 1 , B = X 2 and C = X 3 . The red path corresponds to the model Y = β 0 + β 1 A + β 2 B + +β 3 A 2 + β 4 AB + ϵ.

Figure 2 :

 2 Figure 2: The aggregate Pareto fronts obtained for D 1 (N = 1000; left) and D 2 (N = 10 000; right). Triangles represent the true models.

Figure 3 :

 3 Figure 3: Evolution of the hypervolume for various evaporation rates (averaged in each case over the 50 trials). Dataset: D 1 (N = 1000). Should be viewed in color.

Figure 4 :

 4 Figure 4: Evolution of the Pareto front's average RMSE and NoCoeff values (averaged over 50 trials). Dataset: D 1 (N = 1000).

Figure 5 :

 5 Figure 5: Evolution of the hypervolume when (a) the HVs were initialized under the scheme of Section 4.1, when (b) they were initialized randomly, and when (c) they were not present. Hypervolume was averaged over the 50 trials. Dataset: D 1 (N = 1000). Should be viewed in color.

5 .

 5 Appendix B gives details about data generation in Section 5. Appendix C offers comments about the model search graph. An R script and a package demonstrating the use of MOMSACO are also provided. Comments on the two files are in Appendix D.

Table 2 :

 2 Pareto front size per iteration (averaged over the 50 trials). Dataset: D 1 (N = 1000).

	Iteration Avg. no. of models on the Pareto front
	1	3.6
	3	4.6
	12	5.8
	15	6.0
	24	6.3
	100	6.3
	plateaued earlier for larger evaporation rates, but attained higher values for smaller ones.

Table 3 :

 3 Pareto front data after 100 iterations for different evaporation rates. Figures indicated are averages (standard deviations) over the 50 trials. The rate of contribution was defined as the percentage of models on the true Pareto front that belonged to the MOMSACO front for which the rate was computed. Dataset: D 1 (N = 1000).

	Evaporation	0.1	0.3	0.5	0.7	0.9
	rate
	Hypervolume	804.68 (18.44)	784.11 (23.56)	766.81 (27.75)	751.32 (28.68)	744.63 (28.87)
	Pareto front size	6.30 (0.76)	5.76 (1.06)	5.48 (1.01)	5.04 (1.03)	4.14 (1.14)
	Proportion	of	46.85% (16.25%) 29.61% (18.65%) 16.27% (19.45%) 11.93% (17.26%) 16.20% (21.30%)
	Pareto	optimal
	models
	Rate of contribu-	26.54% (8.58%) 15.09% (9.65%) 8.18% (10.28%)	5.09% (7.15%)	4.91% (5.57%)
	tion to the true
	Pareto front

COMPETING INTERESTS

The authors declare that they have no conflict of interest.


Aymen Bouaziz
email: bouazizaymen18@yahoo.com

LOCAL LIMITS OF GALTON-WATSON TREES CONDITIONED ON LARGE WIDTH

Keywords: Mathematics Subject Classification. 60J80, 60B10, 05C05 Galton-Watson tree, random trees, local limits, width

We study the local convergence of critical Galton-Watson trees under various conditionings. We give a sufficient condition, which serves to cover all the previous cases, for the convergence in distribution of a conditioned Galton-Watson tree to Kesten's tree. We also propose an other proof to give the limit in distribution of a critical Galton-Watson tree, with bounded support, conditioned on having a large width.

Introduction

In [START_REF] Kesten | Subdiffusive behavior of random walk on a random cluster[END_REF], Kesten proved that a critical or sub-critical Galton-Watson (GW) tree conditioned on reaching at least height h converges in distribution (for the local topology on trees) as h goes to infinity toward the so-called sized-biased tree (that we call here Kesten's tree and whose distribution is described in Section 2.3). Since then, several different conditioning have been studied, in particular, the conditioning of extinction after large time, the conditioning of large total population size and the conditioning of large number of leaves. In [START_REF] Abraham | Local limits of conditioned Galton-Watson trees : the infinite spine case[END_REF], Abraham and Delmas provided a criterion for local convergence of finite random trees to Kestens tree, then gave short and elementary proofs of essentially all previous related results and some new ones. Later, in [START_REF] Abraham | Local limits of Galton-Watson trees conditioned of the number of the protected nodes[END_REF], we prove that a critical Galton-Watson tree conditioned on having a large number of marked vertices converges in distribution to Kesten's tree and we then apply this result to give the limit in distribution of a critical Galton-Watson tree conditioned on having a large number of protected nodes.

Let L(t) be the width of the tree t. Remark that this functional L is clearly monotone in the sense of [START_REF] He | Local Convergence of Critical Random Trees and Continuous-State Branching Processes[END_REF]; therefore, using Theorem 2.1 of [START_REF] He | Local Convergence of Critical Random Trees and Continuous-State Branching Processes[END_REF], we immediately get that a critical GW tree τ conditioned on {L(τ ) > n} converges in distribution toward Kesten's tree as n goes to infinity. In this paper, we propose another proof by generalizing somewhat the monotonicity property. Notice that the functional L does not satisfy additivity in [START_REF] Abraham | Local limits of conditioned Galton-Watson trees : the infinite spine case[END_REF]. Thus, considering the conditioning event {L(τ ) = n} is still an open problem.

We say that a probability distribution p = (p(0), p(1), . . .) on nonnegative integers is bounded, if the set {n ∈ N, p(n) > 0} is bounded. For any critical and bounded distribution p satisfying assumption (1), Xin He [START_REF] He | Local Convergence of Critical Random Trees and Continuous-State Branching Processes[END_REF] provides a positive answer to this question, more precisely he proves that a critical GW tree τ conditioned on {L(τ ) = n} converges in distribution toward Kesten's tree as n goes to infinity.

Provide an other proof of this result is the main objective of this paper. On the technical level, our proofs are extremely short and elementary, thanks in particular to the convenient framework in [START_REF] Abraham | Local limits of conditioned Galton-Watson trees : the infinite spine case[END_REF].

The paper is then organized as follows: In section 2, we recall briefly the framework we use for discrete trees and define the Galton-Watson tree τ and Kesten's tree τ * associated with offspring distribution p. In section 3, we state and prove our general result on local convergence of conditioned critical and sub-critical Galton-Watson trees and we apply it to the conditioning of large width, in the critical case, in Corollary 3.6, which is one of our main objective of this paper.

In section 4, we study the conditioning on having large width of a critical Galton-Watson tree, with bounded support, where we give an elementary and short proof. Finally, we generalized this result by considering a type of conditioning never treated.

Technical background on GW trees

2.1. The set of discrete trees.

We denote by N = {0, 1, 2, . . .} the set of non-negative integers and by N * = {1, 2, . . .} the set of positive integers. We recall Neveu's formalism [START_REF] Neveu | Arbres et processus de Galton-Watson[END_REF] for ordered rooted trees. Let

U = n≥0 (N * ) n
be the set of finite sequences of positive integers with the convention (N * ) 0 = {∅}. For u ∈ U, its length or generation |u| ∈ N is defined by u ∈ (N * ) |u| . If u and v are two sequences of U, we denote by uv the concatenation of the two sequences, with the convention that uv

= u if v = ∅ and uv = v if u = ∅. The set of ancestors of u is the set An(u) = {v ∈ U; ∃w ∈ U such that u = vw}.
The most recent common ancestor of a subset s of U, denoted by M (s), is the unique element u of ∩ u∈s An(u) with maximal length |u|. For two distinct elements u and v of U, we denote by u < v the lexicographic order on U i.e. u < v if u ∈ An(v) and u = v or if u = wiu and v = wjv for some i, j ∈ N * with i < j. We write

u ≤ v if u = v or u < v.
A tree t is a subset of U that satisfies:

• ∅ ∈ t.

• If u ∈ t, then An(u) ⊂ t.

• For every u ∈ t, there exists k u (t) ∈ N such that, for every i ∈ N * , ui ∈ t iff 1 ≤ i ≤ k u (t). The vertex ∅ is called the root of t. The integer k u (t) represents the number of offsprings of the vertex u ∈ t, and we call it the outdegree of the node u in the tree t. The maximal outdegree M (t) of a tree t is defined by M (t) = sup{k u (t), u ∈ t} The set of children of a vertex u ∈ t is given by:

C u (t) = {ui; 1 ≤ i ≤ k u (t)}. By convention, we set k u (t) = -1 if u ∈ t. A vertex u ∈ t is called a leaf if k u (t) = 0. We denote by L 0 (t) the set of leaves of t. A vertex u ∈ t is called a protected node if C u (t) = ∅ and C u (t) L 0 (t) = ∅,
that is u is not a leaf and none of its children is a leaf.

For a tree t, we denote by

Z(t) = Card({u ∈ t, |u| = n})
the size of the n-th generation of t, the height of t is defined by

H(t) = sup{|u|, u ∈ t}
and can be infinite. We define the width L(t) of t as:

L(t) = sup k≥0 Z k (t)•
We denote by T the set of trees, by T 0 = {t ∈ T; Card(t) < +∞} the subset of finite trees and by

T 1 = {t ∈ T; lim n-→+∞ |M ({u ∈ t; |u| = n})| = +∞}
the subset of trees with a unique infinite spine.

We say that a sequence of trees (t n , n ∈ N) converges locally to a tree t if and only if lim n→∞ k u (t n ) = k u (t) for all u ∈ U.

Let (T n , n ∈ N) and T be T-valued random variables. We denote by dist(T ) the distribution of the random variable T and write

lim n-→+∞ dist(T n ) = dist(T )
for the convergence in distribution of the sequence (T n , n ∈ N) to T with respect to the local topology. If t, t ∈ T and x ∈ L 0 (t) we denote by

t x t = {u ∈ t} ∪ {xv; v ∈ t }
the tree obtained by grafting the tree t on the leaf x of the tree t. For every t ∈ T and every x ∈ L 0 (t), we shall consider the set of trees obtained by grafting a tree on the leaf x of t:

T(t, x) = {t x t ; t ∈ T}.
For convergence in distribution in the set T 0 ∪ T 1 , we recall the following key characterization in [START_REF] Abraham | An introduction to Galton-Watson trees and their locals limits[END_REF].

Lemma 2.1. Let (T n ) n∈N and T be random trees taking values in the set T 0 ∪ T 1 . Then the sequence (T n ) n∈N converges in distribution to T if and only if:

(1) for every finite tree t ∈ T 0 , lim n-→+∞ P(T n = t) = P(T = t);

(2) for every finite tree t ∈ T 0 and every leaf x of t, lim inf n-→+∞ P(T n ∈ T(t, x)) ≥ P(T ∈ T(t, x)).

Galton Watson trees.

Let p = (p(n), n ∈ N) be a probability distribution on N. We assume that [START_REF] Abraham | Local limits of Galton-Watson trees conditioned of the number of the protected nodes[END_REF] p(0) > 0, p(0) + p(1) < 1, and µ :=

+∞ n=0 np(n) < +∞.
A T-valued random variable τ is a GW tree with offspring distribution p if the distribution of k ∅ (τ ) is p and it enjoys the branching property: for n ∈ N * , conditionally on {k ∅ (τ ) = n}, the subtrees (S 1 (τ ), . . . , S n (τ )) are independent and distributed as the original tree τ .

The GW tree and the offspring distribution are called critical (resp. sub-critical, super-critical) if µ = 1 (resp. µ < 1, µ > 1).

In the critical and sub-critical case, we have that a.s τ belongs to T 0 .

Kesten's tree.

Let p be an offspring distribution satisfying Assumption (1) with µ ≤ 1 (i.e. the associated GW process is critical or sub-critical). We denote by p * = (p * (n) = np(n)/µ, n ∈ N) the corresponding sizebiased distribution.

We define an infinite random tree τ * (the size-biased tree that we call Kesten's tree in this paper) whose distribution is described as follows:

There exists a unique infinite sequence (v k , k ∈ N * ) of positive integers such that, for every h

∈ N, v 1 • • • v h ∈ τ * , with the convention that v 1 • • • v h = ∅ if h = 0. The joint distribution of (v k , k ∈ N *
) and τ * is determined recursively as follows. For each h ∈ N, conditionally given (v 1 , . . . , v h ) and {u ∈ τ * ; |u| ≤ h} the tree τ * up to level h, we have:

This yields a short and an elementary proof of theorem 4.1 in [START_REF] He | Conditioning Galton-Watson trees on large maximal out-degree[END_REF]. In general, assume that p is critical, A satisfies the identity property in [START_REF] Abraham | An introduction to Galton-Watson trees and their locals limits[END_REF] and P(A(τ ) = n) > 0 for any n, Then by theorem 3.1, as n -→ +∞,

dist(τ |A(τ ) = n) -→ dist(τ * )• and as n -→ +∞, dist(τ |A(τ ) ≥ n) -→ dist(τ * )• Remark 3.4.
Assume that p is critical. If A satisfies the monotonicity property in [START_REF] He | Local Convergence of Critical Random Trees and Continuous-State Branching Processes[END_REF], which is equivalent to A(t x t ) ≥ A(t ) ∀t ∈ T, and P(A(τ ) ≥ n) > 0 for any integer n (consider n 0 = +∞), Then, using Theorem 3.1, we have as n -→ +∞:

dist(τ |A(τ ) ≥ n) -→ dist(τ * )•
Thus, the monotonicity property is a particular case of our result, which provide another and a short proof of Theorem 2.1 in [START_REF] He | Local Convergence of Critical Random Trees and Continuous-State Branching Processes[END_REF].

Remark 3.5. As a conclusion, Theorem 3.1 serves to cover all the different cases (identity, additivity and monotonicity) which allows to give a simple and short proof of theorem 2.2.1 and 2.2.4 in [START_REF] Abraham | An introduction to Galton-Watson trees and their locals limits[END_REF] As direct application, we can recover several specific conditionings in the critical case:

(1) Conditioning on extinction after large time, see [START_REF] Kesten | Subdiffusive behavior of random walk on a random cluster[END_REF][START_REF] Abraham | Local limits of conditioned Galton-Watson trees : the infinite spine case[END_REF], (A = H),

as n -→ +∞ dist(τ |H(τ ) ≥ n) -→ dist(τ * )•
(2) Conditioning on the total population size, (A = Card),

As n -→ +∞ dist(τ |Card(τ ) ≥ n) -→ dist(τ * )•
(3) Conditioning on the number of individuals having a given number of children: Let A be a non-empty subset of N. For a tree t, we denote by L A (t) the total number of nodes in the tree t with outdegree in A.

Assume that k∈A p(k) > 0. Then as

n -→ +∞ dist(τ |L A (τ ) ≥ n) -→ dist(τ * )• (4)
Conditioning on the number of protected nodes, see [START_REF] Abraham | Local limits of Galton-Watson trees conditioned of the number of the protected nodes[END_REF], Let A(t) be the number of protected nodes in the tree t. Then as

n -→ +∞ dist(τ |A(τ ) ≥ n) -→ dist(τ * )•
One of our original motivations for this result is the local convergence under the conditioning of large width. So right now we will apply Theorem 3.1 to this specific conditioning.

Since p 0 + p 1 < 1, we have P(L(τ ) ≥ n) > 0 for any n. Then Theorem 3.1 immediately gives the local convergence of critical GW trees to Kesten's tree, under the conditioning of large width.

Recall that the width L(t) of a tree t is defined to be

L(t) = sup k≥0 Z k (t)• Corollary 3.6. Assume that p is critical. Then as n -→ +∞, dist(τ |L(τ ) ≥ n) -→ dist(τ * )• 4.
The conditioning on the largest generation, critical case Theorem 4.1. Let τ be a critical GW tree with offspring distribution p satisfying (1) with bounded support and let τ * be a Kesten's tree associated to p. Let τ n be a random tree distributed according to τ conditionally on {L(τ

) = n}. Then as n -→ +∞, dist(τ n ) -→ dist(τ * )•
where the limit is understood along the infinite sub-sequence {n ∈ N : P(L(τ ) = n) > 0}.

Proof. We consider A(t) = L(t) with n 0 = 1 that is

A n = {t ∈ T; L(t) = n}• Since p(0) + p(1) < 1, the set {n ∈ N : P(L(τ ) = n) > 0} is infinite.
Let t ∈ T 0 , x ∈ L 0 (t) and n be a large enough integer. We consider

K = sup{n ∈ N, p(n) > 0} < +∞
be the supremum of the support of p. We have

∀ |x| ≤ s ≤ H(t) Z s (t) + Z s-|x| (τ ) ≤ L(t) + K s-|x| < n
We deduce that, L satisfies the property identity, that is:

L(t * x τ ) = n ⇐⇒ L(τ ) = n
According to Theorem 3.1, we deduce the convergence and hence end the proof.

We can generalize these results concerning the width by considering: Let A be a non-empty subset of N. For a tree t and s a non-negative integer, we write The case A = {0}, L (s) A (t) represents the number of leaves in generation s of t and S A (t) the maximum number of leaves in the same generation. We can also have L As, for all t ∈ T 0 , x ∈ L 0 (t) and n be a large enough integer, ∀ |x| ≤ s ≤ H(t) L

L

A

  (t) = {u ∈ t; |u| = s and k u (t) ∈ A} the set of individuals, in generation s, whose number of children belongs to A and L

A

  (t) = Z s (t) and so S A (t) = L(t) the largest generation by taking A = N.Since p(0) + p(1) < 1 and k∈A p(k) > 0, the set {n ∈ N : P(S A (τ ) = n) > 0} is infinite.

  ≤ L(t) + K s-|x| < nThe same work as in the previous proof allows to prove Theorem 4.2. Let τ be a critical GW tree with offspring distribution p satisfying (1) with bounded support and such that k∈A p(k) > 0. Let τ * be a Kesten's tree associated to p. Then as n -→ +∞,dist(τ |S A (τ ) = n) -→ dist(τ * )•where the limit is understood along the infinite sub-sequence {n ∈ N : P(S A (τ ) = n) > 0}. anddist(τ |S A (τ ) ≥ n) -→ dist(τ * )•Note that this type of conditioning is never treated.

Acknowledgements. Sincere thanks to anonymous referees for their comments and suggestions, which improved considerably the presentation of this paper.

• The number of children (k u (τ * ), u ∈ τ * , |u| = h) are independent and distributed according to p if u = v 1 • • • v h and according to p * if u = v 1 . . . v h . • Given {u ∈ τ * ; |u| ≤ h + 1} and (v 1 , . . . , v h ), the integer v h+1 is uniformly distributed on the set of integers {1, . . . , k

Remark 2.2. Notice that by construction, a.s. τ * ∈ T 1 has a unique infinite spine. And following Kesten [START_REF] Kesten | Subdiffusive behavior of random walk on a random cluster[END_REF], the random tree τ * can be viewed as the tree τ conditioned on non extinction. Following [START_REF] Abraham | An introduction to Galton-Watson trees and their locals limits[END_REF], for t ∈ T 0 and x ∈ L 0 (t), we have:

(2)

In the particular case of a critical offspring distribution (µ = 1), we get for all t ∈ T 0 and x ∈ L 0 (t):

Main result

Let A be an integer-valued function defined on T which is finite on T 0 and let n 0 ∈ N ∪ {+∞} be given. We define for all n ∈ N * , the subset of trees

Common values of n 0 that will be considered are 1 and +∞.

The following theorem states a general result concerning the local convergence of critical and sub-critical GW tree τ conditioned on A n toward Kesten's tree and its proof is in fact a straighforward adaptation of the proof of theorem 3.1 in [START_REF] Abraham | Local limits of conditioned Galton-Watson trees : the infinite spine case[END_REF].

We denote by

the conditional law of τ given {τ ∈ A n }.

Theorem 3.1. Assume that assumption (1) hold, µ ≤ 1 and that P(τ ∈ A n ) > 0 for n large enough. Then, if for all t ∈ T 0 and x ∈ L 0 (t),

(3) lim inf n-→+∞

As n -→ +∞, We have

Proof. Since µ ≤ 1, we have a.s τ ∈ T 0 and τ * ∈ T 1 . So we will use Lemma 2.1 to prove the convergence.

Using [START_REF] Abraham | Local limits of conditioned Galton-Watson trees : the infinite spine case[END_REF] or else [START_REF] Bertoin | On largest offspring in a critical branching process with finite variance[END_REF], we have for every t ∈ T 0 , x ∈ L 0 (t) and t ∈ T 0 :

Let t ∈ T 0 and x ∈ L 0 (t). For such n we get:

So, we obtain that 2), we deduce that:

, we deduce that:

(5) lim inf

Furthermore, for all t ∈ T 0 and all n > A(t), we have

Finally, by Lemma 2.1, we have proved this result.

We have several remarks related to Theorem 3.1.

Remark 3.2. However, it seems that the above proof is almost copied from the proof of Theorem 3.1 in [START_REF] Abraham | Local limits of conditioned Galton-Watson trees : the infinite spine case[END_REF], it is easy to see that the sufficient condition in Theorem 3.1 is a particular case of this result.

Remark 3.3. We apply Theorem 3.1 to the Conditioning on the maximal out-degree: Let τ be a critical Galton-Watson tree and let n be a large enough integer, we have for every t ∈ T 0 , x ∈ L 0 (t):