Sami Alkhoury
email: sami.alkhoury@univ-grenoble-alpes.fr

Emilie Devijver
email: emilie.devijver@univ-grenoble-alpes.fr

Marianne Clausel
email: marianne.clausel@univ-lorraine.fr

Myriam Tami
email: myriam.tami@centralesupelec.fr

Eric Gaussier
email: eric.gaussier@imag.fr

Georges Oppenheim
email: georges.oppenheim@gmail.com

Smooth And Consistent Probabilistic Regression Trees

We propose here a generalization of regression trees, referred to as Probabilistic Regression (PR) trees, that adapt to the smoothness of the prediction function relating input and output variables while preserving the interpretability of the prediction and being robust to noise. In PR trees, an observation is associated to all regions of a tree through a probability distribution that reflects how far the observation is to a region. We show that such trees are consistent, meaning that their error tends to 0 when the sample size tends to infinity, a property that has not been established for similar, previous proposals as Soft trees and Smooth Transition Regression trees. We further explain how PR trees can be used in different ensemble methods, namely Random Forests and Gradient Boosted Trees. Lastly, we assess their performance through extensive experiments that illustrate their benefits in terms of performance, interpretability and robustness to noise.

Introduction

Classification and regression trees (CART) [START_REF] Breiman | Classification and Regression Trees[END_REF] and the ensemble methods based on them, as Random Forests [START_REF] Breiman | Random forests[END_REF] and Gradient Boosted Trees [START_REF] Friedman | Greedy function approximation: A gradient boosting machine[END_REF][START_REF] Elith | A working guide to boosted regression trees[END_REF], have been successfully used for regression problems in many applications and machine learning competitions. Indeed, in a 2017 survey conducted by Kaggle 1 , decision trees and random forests are respectively the second and third most used machine learning methods in industries after logistic regression. It is however well known (see, e.g., Irsoy et al. [START_REF] Irsoy | Soft decision trees[END_REF], Linero and Yang [START_REF] Linero | Bayesian regression tree ensembles that adapt to smoothness and sparsity[END_REF]) that standard decision/regression trees, based on piece-wise constant functions with hard assignments of data points to regions, may have difficulties to adapt to the smoothness of the link functions as well as to the noise in the input data.

We address these problems in this study by adapting standard regression trees through smooth predictions based on probability functions that relate each data point to each region of the tree. The trees thus obtained can naturally adapt to noisy input and can be shown to be consistent, a property that was not, to the best of our knowledge, established for previous, similar attempts. In addition, we show how to use such trees in ensemble methods as Random Forests and Gradient Boosted Trees. Our contributions are thus fourfold: [START_REF] Biau | Consistency of random forests and other averaging classifiers[END_REF] we introduce new regression trees, called PR trees for Probabilistic Regression trees, that can adapt to noisy dataset as well as to the smoothness of the prediction function relating input and output variables while preserving the interpretability of the prediction and being robust to noise; [START_REF] Breiman | Random forests[END_REF] we prove the consistency of the PR trees thus obtained, [START_REF] Breiman | Classification and Regression Trees[END_REF] we extend these trees to Random Forests and Gradient Boosted Trees and (4) we show, experimentally, their benefits in terms of performance, interpretability and robustness to noise. There is however no free lunch, and these additional properties come with a computational cost, as described in Appendix A.2.

The remainder of the paper is organized as follows: Section 2 discusses the related work. Section 3 then introduces PR trees and their inference. Section 3 presents the main results concerning their consistency, the complete proof of which can be found in the Supplementary Material. Section 4 presents the experiments conducted on PR trees and discusses their extension to Random Forests and Gradient Boosted Trees. Finally, Section 5 concludes the paper.

Related Work

Several researchers have tried to adapt decision trees so as to explicitly take into account the noise and uncertainty present in the input data. For example, Fuzzy Decision Trees [START_REF] Chang | Fuzzy decision tree algorithms[END_REF][START_REF] Yuan | Induction of fuzzy decision trees[END_REF][START_REF] Janikow | Fuzzy decision trees: issues and methods[END_REF][START_REF] Olaru | A complete fuzzy decision tree technique[END_REF], designed for classification purposes, assume that the values for some features and classes are associated with membership functions that allow one to associate an example to different rules and predict several classes (with various degrees) with a rule. The approach presented in Ma et al. [START_REF] Ma | Online active learning of decision trees with evidential data[END_REF] fits within the same framework and aims at reducing data uncertainty by querying adequate data while learning the tree. The possibilistic trees considered in Elouedi et al. [START_REF] Elouedi | Belief decision trees: theoretical foundations[END_REF], Jenhani et al. [START_REF] Jenhani | Decision trees as possibilistic classifiers[END_REF] also follow the same principles, the uncertainty being this time modeled with belief functions. Another approach to model uncertainty in the input is based on Uncertain Decision Trees [START_REF] Tsang | Decision trees for uncertain data[END_REF][START_REF] Qin | Dtu: A decision tree for uncertain data[END_REF][START_REF] Liang | Decision tree for dynamic and uncertain data streams[END_REF], also designed for classification, that assume that examples can take on several values through explicit probability density functions (pdf) on given intervals (or explicit probability tables for categorical values). An example may belong, with a certain probability, to several leaves depending on its value range. All these approaches rely on assumptions on the information available for each input (membership scores and pdfs or probability tables) and do not directly aim at adapting to the smoothness of the true prediction function.

Closer to our objectives are Soft trees, introduced in Irsoy et al. [START_REF] Irsoy | Soft decision trees[END_REF], Fuzzy Trees, introduced in Suarez and Lutsko [START_REF] Suarez | Globally fuzzy decision trees for classification and regression[END_REF], and Smooth Transition Regression trees (STR trees), introduced in [START_REF] Da Rosa | Tree-structured smooth transition regression models[END_REF]. Soft and Fuzzy trees are very close to each other: they can both be used for classification and regression, and learn a parameter vector at each node, the dimensionality of which equals the one of the input data. This vector is used in a gating function, based on the sigmoid function in both cases, assigning, to each example, a probability to branch to the left and right children of the node. Each example is thus assigned to all leaves with a certain classmembership, and the final prediction is a smooth combination of the prediction at each node. Soft and Fuzzy trees can be seen as a direct extension of hierarchical mixtures of experts (HMEs, Jordan and Jacobs [START_REF] Jordan | Hierarchical mixtures of experts and the em algorithm[END_REF]): indeed, if HMEs either use predefined trees or trees learned by another method (typically a standard decision/regression tree), soft trees are constructed while learning the hierarchy of experts. STR trees follow the same general principle but rely on a single parameter at each node. A sigmoid-based gating function is also used to assign points to different regions of the tree. Adaptative Neural Trees [START_REF] Tanno | Adaptive neural trees[END_REF] and Deep Neural Decision Forests [START_REF] Kontschieder | Deep neural decision forests[END_REF] are both built from decision trees. These models are very close to soft trees; however, the models are enhanced with a neural network representation and suffer from a lack of interpretability (one can even argue that these models are not tree models per se). The main differences between these approaches and PR trees lie in the way the assignment of examples to leaves is done, in the number of parameters used by each method and in their theoretical guarantees. By directly relying on a probability distribution, PR trees rely on less parameters which make them a priori more robust to overfitting. In addition, we show here that PR trees are consistent; we know of no such results for the other models. Frosst and Hinton [START_REF] Frosst | Distilling a neural network into a soft decision tree[END_REF] consider a specific variant of the soft tree model, with knowledge distillation, which is beyond the scope of this work.

Lastly, one of the objectives of this study is to use PR trees in ensemble methods. The most wellknown ensemble methods built on top of regression trees are certainly Random Forests (RF) [START_REF] Breiman | Random forests[END_REF] and Gradient Boosted Trees (GBT) [START_REF] Friedman | Greedy function approximation: A gradient boosting machine[END_REF][START_REF] Elith | A working guide to boosted regression trees[END_REF]. We show in Section 4 how to use PR trees in these ensemble methods and illustrate their behavior compared to standard RF, GBT and BooST [START_REF] Fonseca | Boost: Boosting smooth trees for partial effect estimation in nonlinear regressions[END_REF], a boosting-based ensemble extension of STR trees. More recently, tree-based Bayesian ensemble methods, as BART [START_REF] Chipman | Bart: Bayesian additive regression trees[END_REF] and soft-BART [START_REF] Linero | Bayesian regression tree ensembles that adapt to smoothness and sparsity[END_REF], have been proposed. If these methods are interesting, they fit within a (more complex) Bayesian framework that needs to be adapted to the different trees used (this contrasts with Random Forests and Gradient Boosted Trees which can readily integrate different regression trees). Their consideration is beyond the scope of the current study.

Probabilistic Regression Trees

Let X = (X 1 , • • • , X p ) be a p-dimensional input random vector lying almost surely in a compact subspace X of R p , and let Y be an output random variable linked to X through:

Y = f (X; Θ) + ε Y , ε Y ∼ N (0, σ2 ), (1)
where Θ is the set of parameters on which f relies.

Due to their success, we are interested in this study in approximations of f based on regression trees, either considered in isolation or aggregated in sum-of-trees models as Random Forests or Gradient Boosted Trees. Regression trees [START_REF] Breiman | Classification and Regression Trees[END_REF] define a partition of R p into K hyper-rectangles, referred to as regions and denoted

R k = p =1 [a k, , b k, ],
∀k, 1 ≤ k ≤ K obtained by dyadic splits2 . A weight γ k is associated to the k-th region R k leading, for observations x ∈ R p , to a predictor of the form:

f (x; Θ) = K k=1 γ k 1 {x∈R k } .
As noted in previous studies [START_REF] Irsoy | Soft decision trees[END_REF][START_REF] Linero | Bayesian regression tree ensembles that adapt to smoothness and sparsity[END_REF], regression trees, as they are based on constant piece-wise functions, may fail to accommodate the smoothness of the link function. To solve this problem, we introduce a prediction function that generalizes the one of standard regression trees by replacing its indicator function with a smooth function Ψ:

T PR (x; Θ) = K k=1 γ k Ψ(x; R k , σ). (2)
The set of parameters Θ = ((R k ) 1≤k≤K , γ, σ) corresponds to the set of regions, associated weights represented by γ ∈ R K and noise in the input variables captured in σ ∈ R p + . In the remainder, we will refer to σ as the noise vector. Note that when

Ψ(x; R k , σ) = 1 {x∈R k } , ∀k, 1 ≤ k ≤ K,
one recovers standard regression trees, and so the consideration of σ extends the standard framework for regression trees, which does not explicitly take into account the fact that the input may be noisy.

The functions Ψ we consider relate, through a sufficiently regular probability density function φ (φ is in L 2 , is of class C 1 and the support of its Fourier transform is R p ), data points to different regions of the tree and smooth the predictions made. They are defined, for all x ∈ X , by:

Ψ(x; R k , σ) = 1 p j=1 σ j R k φ u j -x j σ j 1≤j≤p du. (3)
In practice, experts may have an empirical knowledge on the nature of the errors (for example when measurements are done by calibrated machines) which can help to choose φ. We will make use of different distributions on different datasets in our experiments. As one can expect, the best choice for φ, which can be established via cross-validation, depends on the collection considered. Standard regression trees are obtained by considering a Gaussian distribution for φ with σ j → 0, ∀j, 1 ≤ j ≤ p.

In that case, the distribution of x over regions is concentrated on one region.

Parameter estimation Given a training set D n = (x (i) , y (i) ) 1≤i≤n , with x ∈ R p , y ∈ R, and in accordance with the empirical risk minimization principle with a quadratic loss, the estimation procedure for probabilistic regression trees aims at finding the parameters Θ solutions of:

argmin Θ n i=1 y (i) - K k=1 γ k P ik 2 , with P ik := Ψ(x (i) ; R k , σ). (4)
The n × K matrix P = (P ik ) thus encodes the relations between each training example x (i) and each region R k . It is such that 0 ≤ P ik ≤ 1 and ∀i, 1 ≤ i ≤ n,

K k=1 P ik = 1.
The estimation of the different parameters in Θ alternates in between region and weight estimates, as in standard regression trees, till a stopping criterion is met 3 . During this process, the number of regions is increased and the matrix P and the weights γ are gradually updated.

Estimating γ -When fixing the regions (R k ) 1≤k≤K and the vector σ, minimizing Eq. ( 4) with respect to γ leads to the least square estimator defined by γ = P T P -1 P T y if P T P is not singular. We discuss in the Supplementary Material some conditions on the invertibility of P T P. Estimating (R k ) 1≤k≤K -Let us assume that K regions, referred to as current regions, have already been identified, meaning that the current tree has K leaves. As in standard regression trees, each current region R k , 1 ≤ k ≤ K, can be decomposed into two sub-regions with respect to a coordinate 1 ≤ j ≤ p and a splitting point s j k that minimizes Eq. ( 4). Each split leads to the update of P, that now belongs to M n,K+1 (R), and γ, that now belongs to R K+1 . Substituting γ by its value, the best split for the current region R k is given by:

argmin 1≤j≤p,s∈S j k n i=1 y (i) - K+1 l=1 P T P -1 P T y l P il 2 , (5)
where S j k denotes the set of splitting points for region R k and variable j (more precisely, S j k is the set of middle points of the observations from R k projected on the jth coordinate).

Estimating σ -Lastly, the vector σ can either be based on a priori knowledge or be learned through a grid search on a validation set. We rely on the latter in our experiments.

Illustration As an illustration, we consider a toy example based on Y = cos(X) + ε, with X ∼ U([0, 5]) and ε ∼ N (0, 0.05 2 ). A training set of 100 observations is generated from this model.

For the PR tree, the function Ψ is set to the Gaussian distribution 4 , for 1 ≤ k ≤ K,

Ψ(x; R k , σ) = 1 σ √ 2π b k a k e -(u-x) 2 2σ 2 du. (6)
where the variance is fixed to σ 2 = var((x (i) ) 1≤i≤n )/2 = 0.74 2 . Figure 1 (left) compares the performance of standard tree and Probabilistic Regression tree to the true regression function. Unlike the standard regression tree, which is based on constant piece-wise functions, the PR tree, given in Figure 1 (middle), is able to accurately approximate the true function, even though few observations are available. Lastly, Figure 1 (right) illustrates the weights γ as well as the distribution of Ψ across the regions. As σ 2 is relatively large in this case, there is a large overlap between the functions Ψ across regions, leading to a smooth prediction function.

Consistency We show here that the PR tree learned from a training set of size n, denoted

T (n) PR , is consistent, that is lim n→+∞ E[| T (n) PR (X) -E(Y |X)| 2 ] = 0.
All the proofs, which are partly based on the density results described in Schaback [START_REF] Schaback | Multivariate interpolation and approximation by translates of a basis function[END_REF], Devore and Ron [START_REF] Devore | Approximation using scattered shifts of a multivariate function[END_REF] and the consistency results established in Györfi et al. [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF], Scornet et al. [START_REF] Scornet | Consistency of random forests[END_REF], are fully detailed in the Supplementary Material.

To emphasize the fact that the trees learned depend on a training set of size n, we denote, in this section as well as in the Supplementary Material, the k-th region as R

(n) k with k = 1, • • • , K n .
In the sequel, we assume that the sequence (K n ) is non decreasing. Furthermore, for the sake of simplicity, we consider that the observations lie in [0, 1] p , the extension to any compact subspace X of R p being direct. We also denote the Fourier transform of any function f as Ff and assume that the noise vector σ is fixed and known 5 .

The function Ψ at the basis of PR trees (Eq. ( 3)) relies on a probability density function φ that regulates how each data point is distributed across all regions of the tree. We assume here that φ satisfies the following conditions: the support of its Fourier transform is

R p ; ∃r > 0, sup v∈R p |v| 1+r+p/2 |φ(v)| < ∞ (which ensures that φ ∈ L 2 ) ; and φ ∈ C 1 .
The estimators provided by PR trees belong to the function space V n defined by:

V n = Kn k=1 A k Ψ(•; R (n) k , σ), (A k ) ∈ R Kn .
As usual in consistency theory [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF], we first focus on truncated versions of the different functions. From the remainder of this section, let (β n ) be a sequence of non negative numbers such that lim n→+∞ β n = +∞. We define

T (n) PR = T βn T (n)
PR the truncated estimator as:

T βn u = u if |u| < β n , sign(u)β n if |u| ≥ β n .
We also define T βn V n = {g ∈ V n : g ∞ ≤ β n }, and consider the Sobolev space of functions defined, for s ∈ (1, 2), by

H s ([0, 1] p ) = {f ∈ L 2 ([0, 1] p ), ∃g ∈ H s (R p ) s.t. f = g| [0,1] p }, with H s (R p ) = {f ∈ L 2 (R p ), (1 + || • || 2 2 ) s 2 |Ff (•)| ∈ L 2 (R p )}. Note that H s ([0, 1] p )
, that contains functions from L 2 for which the derivative of order s is also in L 2 , encompasses many usual functions.

As detailed in the Supplementary Material, T βn V n is dense in H s ([0, 1] p ):

Proposition 1. Assume that for some s ∈ (1, 2), E(Y |X = •) ∈ H s ([0, 1] p ).
Let (λ n ) be a non decreasing sequence such that for n sufficiently large ,

β n > λ p/2 n E(Y |X = •) H s ([0,1] p ) inf |ω|≤λn |Fφ(ω)| . (7)
Then, one has for an absolute constant C and for any M > 0:

inf

g∈T βn Vn E(Y |X = •) -g 2 L 2 ([0,1] p ) ≤ C E(Y |X = •) 2 H s ([0,1] p ) a n , (8)
with

a n = 1 (1 + |λ n | 2 ) s/4 + β n M -r + β n max k=1,...,Kn diam(R (n) k ∩ [-M, M ] p ) 2 .
Note that in whole generality, our regions R

k may be unbounded, whereas we want to approximate our unknown link function E(Y |X = •) on the compact set where lie the observations. Proposition 1 states that, up to a controlled error β n M -r , we can replace each region by its restriction to a given compact [-M, M ] p .

For the estimation error, one can first bound the difference between the empirical and real errors for functions in T βn V n through: Proposition 2. Let (X (1) , Y (1) ), • • • , (X (n) , Y (n) ) be i.i.d. copies of the vector (X, Y ). For any n ∈ N, any g ∈ T βn V n and for any η > 0,

E sup g∈T βn (Vn) 1 n n i=1 |g(X (i) ) -T βn Y (i) | 2 -E[|g(X) -T βn Y | 2 ] ≤ η + 64β 2 n exp [-c n ] , with c n = n β 4 n η 2 2048 - β 4 n (K n -1) log(pn) n - 2β 4 n K n n log C β 2 n η ,
and where C > 0 is a universal constant.

Let us now set

β n = E(Y |X = •) ∞ + σ√ 2(log n) 2 .
Consider λ n such that Eq (7) is satisfied. Moreover, assume that: K n -----→ n→+∞ +∞, and K n (log n) 9 /n -----→ n→+∞ 0, which means that the number of regions needed to explain an infinite number of data points is also infinite, however still largely (by a factor (log n) 9 ) dominated by the number of points. Then the above results lead to the consistency of the truncated estimator T (n) PR . The extension to the consistency of the untruncated estimator

T (n) PR follows from lim n→+∞ E[| T (n) PR (X) - T (n) PR (X)| 2 ]
= 0, derived from arguments similar to the ones developed in Scornet et al. [START_REF] Scornet | Consistency of random forests[END_REF]. This finally leads to: Theorem 1. Set M > 0. With K n as above, assume that for some s

∈ (1, 2), E(Y |X = •) ∈ H s ([0, 1] p ) and that max k=1,••• ,Kn diam(R (n) k ∩ [-M, M ] p ) -----→ n→+∞ 0. Then lim n→+∞ E[| T (n) PR (X) -E(Y |X)| 2 ] = 0.
The assumption on the diameter of the regions is reasonable for data points lying in a compact subspace: as the number of regions grows (to infinity) with the number of data points, their diameter will decrease, data points being more and more concentrated in each region. An important difference with respect to Biau et al. [START_REF] Biau | Consistency of random forests and other averaging classifiers[END_REF], Scornet et al. [START_REF] Scornet | Consistency of random forests[END_REF] is that we are not averaging over independent classifiers as regions are dependent on each other. Lastly, note that the only assumption made on the link function is that it belongs to a particular Sobolev space; no modeling assumption, as the additive assumption made in Scornet et al. [START_REF] Scornet | Consistency of random forests[END_REF], is required here.

Experimental Validation

We illustrate in this section the behavior of PR trees on several collections, and compare them with standard regression trees, Soft regression trees [START_REF] Irsoy | Soft decision trees[END_REF] and STR trees [START_REF] Da Rosa | Tree-structured smooth transition regression models[END_REF].

Methods For Soft trees, we use the implementation available in github 6 with the default parameters.

For STR trees and BooST, its extension to GBT, we use the implementation available in github 7 . For standard regression trees and their ensemble extensions, we use the implementation from Scikit-Learn [START_REF] Varoquaux | Scikit-learn: Machine learning without learning the machinery[END_REF]. PR trees are built on top of this implementation, and the code for this implementation is available in https://gitlab.com/sami.kh/pr-tree. To compute the weights γ, we rely on the Moore-Penrose pseudo-inverse. In addition, unless otherwise specified, we use the Normal distribution for Ψ (Eq. 3). Lastly, for both PR and standard regression trees, the stopping criterion is the same in all experiments: all leaves should contain at least 10% of the training data. The stopping criterion for STR trees is based on the number of regions, which is chosen to be equal to the ones obtained for PR and standard trees. For Soft trees, the growth is stopped when the performance on a validation set decreases, which gives them a slight advantage in terms of performance. All the results are evaluated using the Root Mean Squared Error (RMSE).

Complexity As the complexity for training PR trees is O(pSK 2 n), where S is a constant upper bounding the number of considered splitting points (see the Supplementary Material), we consider a subset of splitting variables consisting of the top V variables according to the splitting criterion of standard regression trees (this last step is negligible when the tree is of depth 2 or more). This finally leads to an overall complexity of O(V SK 2 n). In the remainder, V is set to 3. Note that the results obtained with this and higher values are close to the ones obtained when considering all variables (see the Supplementary Material).

Data sets We make use here of 13 data sets of various size, namely (ordered by increasing sample size) Riboflavin (RI), Ozone (OZ), Diabetes (DI), Abalone (AB), Boston (BO), Bike-Day (BD), E2006, Skill (SK), Ailerons (AL), Bike-Hour (BH), Super Conductor (SC), Facebook Comments (FC) and Video Transcoding (VT), all commonly used in regression tasks. In the experiments reported here, we use the original data sets, without any modification, to illustrate the fact that one can gain by treating real data sets as noisy. Full details on these datasets (sample size and number of variables, location) are given in the Supplementary Material.

Experiments We use stratified 10-fold cross-validation to evaluate the performance of each method. Each fold is divided into 80% for train and 20% for test, except for Soft trees and PR trees and their gradient boosted extension (see below) for which each fold is divided into 65% for train, 15% for validation and 20% for test. For Soft trees, the validation set is used for the stopping criterion. For PR trees and their gradient boosted extension, it is used to estimate the noise vector σ through a grid search taking values, for each variable j, 1 ≤ j ≤ p, in the interval [0, 2σ j ] with a step of σj /4, where σj denotes the empirical standard deviation of variable j. All experiments were conducted on a 256 GB RAM server with 32 CPUs at 2.60GHz. Lastly, a two-sided t-test at the level 0.05 is used to assess whether differences are significant or not. Results for tree-based methods Figure 2 displays the results obtained with the different regression trees (full performances stand in Supplementary Material). As one can note, PR trees significantly outperform standard, Soft and STR trees, on almost all data sets. Standard regression trees never significantly outperform PR trees whereas Soft trees are only significantly better than PR trees on one collection. In contrast, PR trees are significantly better than Soft trees on eight collections and significantly better than standard and STR trees on eleven collections. When the number of examples is sufficiently large, e.g., on datasets containing more than 5,000 examples (as AL, BH, SC, FC, VT), the difference between PR and Soft and STR trees is more important. This indicates that the models learned by Soft and STR trees, which are more complex as they rely on more parameters, have a tendency to overfit in this case. This indicates that the models learned by Soft and STR trees, which are more complex as they rely on more parameters, have a tendency to overfit in this case. Lastly, when the number of observations is small compared to the number of features, as on RI and E2006 Interpretability, adaptability and robustness As mentioned before, PR trees can be easily interpreted. In addition, they are adaptable in the sense that different distribution functions can be used to smooth the prediction, and they are relatively robust to noise. We illustrate these points below and refer the reader to the Supplementary Material for more details.

Interpretability To illustrate how observations are linked to the different regions, one can see in the figure on the right, for OZ, DI and BD, the boxplots for all observations of the probability to belong to the three most probable regions (K * 1 denotes the most probable region for any observation, K * 2 the second most probable region and K * 3 the third most probable region). For the three data sets and for most observations, roughly 75% of their distribution is concentrated on these three regions, which are by far the most important ones and can be used to provide a first explanation for the values predicted.

Adpatability to the underlying probability distribution Another advantage of PR trees is that different distribution functions can be used to smooth the prediction (function φ in Eq. 3). The choice of φ can be made according to some a priori knowledge on the nature of the errors or by testing different distributions and selecting the best one on a validation set. We provide in the Supplementary Material the results of an experiment conducted on seven datasets with six different distributions (two variants of the Gamma distribution, the Laplace distribution, the Lognormal distribution, the Normal distribution and two variants of the Student distribution). As one can expect, the choice of the best distribution depends on the collection, the Student distribution with 3 degrees of freedom being, for example, particularly adapted to the AL, BO and DI collections. This said, the Normal distribution behaves well on all collections, hence its choice in our experiments.

Robustness Because they allow data points to be probabilistically assigned to different regions, PR trees are more robust to noise than their standard counterpart. To illustrate this point, we added, for each variable j, 1 ≤ j ≤ p, some random Gaussian noise N (0, σ 2 τ,j ) with σ τ,j = (τ σj ) 1≤j≤p , τ ∈ {0.1, 0.25, 0.5, 1, 1.25}, corresponding to 5 noise levels on the empirical standard deviation σj . We then evaluated how the RMSE evolves according to the noise level on three collections, OZ, DI and BD. The results are summarized in Figure 3. As one can note, PR trees are more robust to noise than standard trees, and to a lesser extent than Soft trees. STR trees are relatively stable as well (as Soft trees on two datasets) but with worse performance, limiting their interest.

Ensemble extensions

The extension of PR trees to Random Forests, denoted by PR-RF is straightforward. Note that the noise vector is fixed to the one obtained for a single tree as each tree aims at predicting the output variable Y . On all the collections described above, we ran both PR and standard Random Forests with 100 trees. The default parameters (consisting of all variables with a sampling strategy consisting of bootstrap with replacement) are used in both cases. The results obtained (presented in the Supplementary Material) show that PR-RF significantly outperform their standard counterpart on seven collections (out of 13), are on a par on five collections and are significantly worse on only one collection. Random Forests aim at reducing the variance (at the expense of a small increase in the bias), whereas our adaptation to smoothness aims at reducing the bias. Combining both, as in PR-RF, reduces both bias and variance and leads to a method which significantly outperforms RF.

For using PR trees in Gradient Boosted Trees, an extension we refer to as PR-GBT, assuming that (T -1) trees have been built so far, we simply replace y (i) in Eq. 4 by its residual (y (i) -T -1 t=1 Kt k=1 γ k P ik ). One can directly use the procedure given in Section 3 to estimate the new tree in the ensemble. However, this time, the noise vector may change from one tree to another as a different output is predicted. We thus tune the noise vector on each tree independently, using the same grid search as the one detailed in the previous section. Figure 4 compares PR-GBT with standard Gradient Boosted Trees, denoted by GBT, and BoosT, the extension of STR trees to Gradient Boosted Trees. As one can note, PR-GBT is the best performing method overall, even though all methods are on a par on five collections and BoosT obtains results close to PR-GBT (the latter being significantly better than the former on 4 collections, and significantly worse on 3 collections). Figure 4: Normalized results obtained with 10 stratified cross-validation on PR GBT, standard GBT and BooST. '+' means the best result, and '*' corresponds to results are not significantly worst than the best one, according to a t-test at 5%. Values above 1.1 and below 0.9 have been truncated.

Conclusion

We have studied here approximations, based on regression trees, of smooth link functions. Standard regression trees, while being a basic building block for two of the most popular and efficient ensemble methods, Random Forests and Gradient Boosted Trees, are based on constant piece-wise functions and may fail to accommodate the smoothness of the link function. To solve this problem, we have introduced functions that relate, through sufficiently regular probability density functions, data points to different regions of the tree and smooth the predictions made. We have then shown that the probabilistic regression trees thus obtained are consistent, i.e., lim n→+∞ E[ T (n) s (X) -E(Y |X)] 2 = 0, with T (n) s the probabilistic regression tree learned from a training set of size n, and outperform, on a variety of collections, previously proposed trees (as standard, Soft and Smooth Transition Regression trees) in terms of performance and robustness to noise. Lastly, we have proposed versions of Random Forests and Gradient Boosted Trees based on Probabilistic Regression trees and shown that these versions outperform the state-of-the-art.

In the future, we want to extend the consistency results of probabilistic regression trees to their ensemble extensions, and derive Bayesian versions of these extensions, following the work conducted in Chipman et al. [START_REF] Chipman | Bart: Bayesian additive regression trees[END_REF] and Linero and Yang [START_REF] Linero | Bayesian regression tree ensembles that adapt to smoothness and sparsity[END_REF]. We also plan to investigate knowledge distillation, as introduced in Frosst and Hinton [START_REF] Frosst | Distilling a neural network into a soft decision tree[END_REF].

Broader Impact

This work is mainly theoretical, with standard regression applications.

Figure 1 :

 1 Figure 1: Left: plot of the true regression function, the probabilistic regression tree prediction and the standard regression tree prediction. Middle: PR tree learned from 100 observations where the stopping criteria consists in having at least 20% of the observations in each leaf. Right: description of the probabilistic regression tree parameters: regions (corresponding to intervals as p = 1), γ in each region and the functions Ψ associated to each region.

Figure 3 :

 3 Figure 3: RMSE for the different trees across 5 noise levels on OZ (left), DI (middle) and BD (right).

  Normalized results obtained with 10 stratified cross-validation on PR, standard, Soft and STR trees. '+' means the best result, and '*' corresponds to results that are not significantly worst than the best one according to a two-sided t-test at 5%.

		2
																	PR Tree	Std Tree	Soft Tree	STR
	Relative RMSE	0.5 1 1.5	*	+	+	*	*	+	*	**	+	+	*	+	+	+	+	+	+	+	*	+
		0
				RI		OZ		DI		AB			BO	BD	E2006	SK	AL	BH	SC		FC	VT
	Figure 2:

34th Conference on Neural Information Processing Systems (NeurIPS 2020), Vancouver, Canada.

The State of Data Science & Machine Learning; https://www.kaggle.com/surveys/2017

For 1 ≤ k ≤ K, (a k,j , b k,j ) ∈ (-∞, +∞)[START_REF] Breiman | Random forests[END_REF] and the segments are extended in the obvious way when a k,j = ±∞ or b k,j = ±∞, 1 ≤ j ≤ p.

Any standard stopping criterion, as tree depth or number of examples in a leaf, can be used here.

Note that in this case, defining the distance between an observation x and a regionR k by d(x, R k ) = inf z∈R k ||x -z|| 2 2, the closer x is to region R k , the most important the contribution of R k is in the prediction of Y given x.

Note that this assumption is in line with the procedure we rely on to estimate σ.

https://github.com/oir/soft-tree

https://github.com/gabrielrvsc/BooST

Acknowledgments and Disclosure of Funding

This work has been partially supported by MIAI@Grenoble Alpes (ANR-19-P3IA-0003) and by the french PIA project Lorraine Université d'Excellence (ANR-15-IDEX-04-LUE).


Konstantinos Nikolakopoulos

Jean-Philippe Crete
email: jean-philippe.crete@supmeca.fr

Patrice Longere

Volume averaging based integration method in the context of XFEM-cohesive zone model coupling

Keywords:

a b s t r a c t

The main issue of the extended finite element method (XFEM) is the numerical integration of the system of equilibrium equations. Indeed, in order to have a correct displacement jump vector, the integration needs to be achieved on both sides of the discontinuity and thus requires the existence of integration points on both sides of the discontinuity. A volume averaging based integration method is developed in the present work alleviating this constraint and applied to XFEM coupled with cohesive zone model in a three-dimensional formulation. Moreover, unlike other widely used integration methods, the proposed method does not require the a priori knowledge of the position of the discontinuity inside the finite element nor the projection of the state variables.

Introduction

The main issue of the extended finite element method (XFEM) is the numerical integration of the system of equilibrium equations, in particular when XFEM is coupled with a cohesive zone model (CZM). Indeed, CZM is generally governed by a relationship between the displacement jump vector across the discontinuity and the cohesive traction forces that prevent the crack from opening. If there is no integration point on one of the two sides of the crack, classic Gauss quadrature does not provide an accurate enough integration leading to wrong calculations of the stiffness matrix, the internal forces and then the displacement jump vector. In literature there exist several approaches aiming at palliating this deficiency.

The most commonly used technique is the subdivision of elements, with examples in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF][START_REF] Mariani | Extended finite element method for quasi-brittle fracture[END_REF] , wherein the element sub-volumes V + and V -, separated by the crack, are triangulated into subdomains, viz. triangles (if R 2 ) or tetrahedra (if R 3 ), each containing one or more Gauss points. The standard Gauss quadrature is then applied for the integration of the weak form in every subdomain. This method should not be confused with the mesh modification techniques, see [START_REF] Mediavilla | A robust and consistent remeshing-transfer operator for ductile fracture simulations[END_REF] , as the said partitioning is realized only for the integration and thus does not introduce any supplementary degrees of freedom. The element subdivision based technique is well adapted for elastic materials or for structures with initial defects. It implies the projection of state variables on the new integration points. This procedure can be both computationally cumbersome and leads to errors, especially in the case of elasto-plastic materials whose response is path-dependent.

To address the aforementioned issue, Elguedj et al. [START_REF] Elguedj | Appropriate extended functions for x-fem simulation of plastic fracture mechanics[END_REF] instead proposed to raise significantly the number of Gauss points (64 in 2D), as a means to ensure that there is always at least one integration point on both sides of the discontinuity (if the discontinuity crosses the element near a node, the discontinuity is slightly shifted without substantial modification of the discontinuity path). With this technique, there is no projection requirement but the extension in 3D requires 512 integration points, increasing drastically the cost and the time of calculation.

Martin et al. [START_REF] Martin | About the use of standard integration schemes for x-fem in solid mechanics plasticity[END_REF] proposed the substitution of the discontinuous Heaviside function, in the elements cut by the discontinuity, by a new continuous enrichment function. Using this technique, there is no issue of projection of the state variables and no modification of the crack path, for a number of Gauss point of 27 for a linear hexahedron, which also leads to an increase in the computation time. Moreover, the authors do not discuss the use of the B-Bar approach (see [START_REF] Hughes | Generalization of selective integration procedures to anisotropic and nonlinear media[END_REF] ) to prevent volumetric locking in the case of a strongly nonlinear behavior material.

The volume averaging based integration method developed in the present work is inspired from the approach developed in Belytschko [START_REF] Belytschko | A finite element with embedded localization zones[END_REF] . The method allows for alleviating the need for the existence of integration points on both sides of the discontinuity and is applied to XFEM coupled with cohesive zone model in a https://doi.org/10.1016/j.mechrescom.2020.103485 three-dimensional formulation. Moreover, unlike other widely used integration methods, the proposed method does not require the a priori knowledge of the position of the discontinuity inside the finite element nor the projection of the state variables.

In Section 2 the formulation of the X-FEM and the coupling with a cohesive zone model are recalled. The numerical integration of the equilibrium equations, including the B-Bar approach, is detailed in Section 3 . Two applications, viz. with and without initial discontinuity, are presented in Section 4 , in view of comparing the proposed technique with classical approaches. Concluding remarks are given in Section 5 .

X-FEM -CZM coupling

The eXtended Finite Element Method [1 , 9] has become one of the most widely used methods for the simulation of crack propagation in engineering structures. It is based on the idea of embedding the crack within the finite element. The kinematics of the crack is accounted for by enriching the regular displacement field of the finite element with additional degrees of freedom. One of its most important advantages is that the crack can propagate independently of the meshing, without the need for the a priori knowledge of the crack path. XFEM has demonstrated its performances for brittle and quasi-brittle materials, see [3 , 10] . Now, in the case of strongly non-linear elastoplastic materials subject to ductile damage, this method may lead to unrealistic results, see [START_REF] Crété | Numerical modelling of crack propagation in ductile materials combining the GTN model and X-FEM[END_REF] . Indeed, the activation of the enriched dofs in an element leads to a sudden (traction free) crack opening. Yet, in these materials, the failure should be gradual. To tackle this problem, the concept of combining XFEM with a cohesive segment [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF] is adopted in the present work.

X-FEM formulation

In this work, the "shifted basis" formulation of the X-FEM [START_REF] Zi | New crack-tip elements for xfem and applications to cohesive cracks[END_REF] without singular functions is used. The discretized total displacement field u(x) accordingly reads

u ( x ) = n N i a i + m H ( x ) -H j N j b j (1.1)
where Ni is the i th standard FE shape function, n the number of nodes, a i the i th standard displacement degree of freedom, m the number of enriched nodes and b j the j th additional degree of freedom associated to the j th node. H is the Heaviside function, which equals to + 0.5 if a point is located at the "positive" side of the discontinuity and -0.5 else and H j the value of the Heaviside function at each node.

The coupling between XFEM and CZM [START_REF] Wolf | Numerical modeling of strain localization in engineering ductile materials combining cohesive models and X-FEM[END_REF] consists in inserting a cohesive segment in the strong discontinuity in order to form a "cohesive strong discontinuity". The XFEM formulation is adopted and the system to be solved reduces to the following linear form:

[ K ] da db = - ⎧ ⎨ ⎩ ∫ V e B T σdV ∫ V e B * T σd V + ∫ D N T T loc d ⎫ ⎬ ⎭ , ( 1.2)
where

[ K ] = ⎡ ⎣ ∫ V e B T D t B dV ∫ V e B T D t B * dV ∫ V e B * T D t B dV ∫ V e B * T D t B * d V + 4 ∫ D N T C loc Nd ⎤ ⎦ (1.3)
Where σ is the stress vector, D t the elastic-plastic tangent operator, T loc the cohesive tractions vector, C loc the cohesive stiffness matrix, B the spatial derivatives matrix of the shape functions, B * the spatial derivatives matrix of the ( H ( x ) -H j ) N j functions and D the discontinuity. The subscript loc refers to the local frame of the localization band.

Adopted CZM approaches

As can be seen in Eqs. (1.2) and ( 1.3 ), the coupling of CZM with the XFEM implies the incorporation of the cohesive tractions vector T loc and the cohesive stiffness matrix C loc , in the finite element formulation. The method of calculation of the above physical quantities, via the traction separation law (TSL), strongly depends on the material and engineering application under consideration [START_REF] Li | Analysis of crack growth and crack-tip plasticity in ductile materials using cohesive zone models[END_REF] . In the present work, two different applications are considered ( Sections 4.1 and 4.2 ), leading to the use of two different cohesive laws ( Sections 2.2.1 and 2.2.2 ).

Both laws are activated as soon as the deterioration criterion is met. Their common basis is inspired from the approach in Camanho [START_REF] Camanho | Mixed-mode decohesion finite elements for the simulation of delamination in composite materials[END_REF] where the evolution of the softening branch of the cohesive law is governed by a de-cohesion variable D ∈ [0,1], where 0 denotes a healthy element and 1 a fully damaged element. The argument of the function D = f ( ) is an equivalent displacement that takes into account the displacement jump vector δ, with δ 1 the normal component and δ 2,3 the shear components (expressed in the local frame of the localization band):

= δ 1 2 + δ 2 2 + δ 2 3
, where δ 1 = max ( 0 , δ 1 )

(1.4)
For the sake of simplicity, contact is not treated here. The evolution law of the de-cohesion variable is assumed herein of the form

D = 0 , i f ≤ 0 -0 c , i f 0 < < 0 + c 1 , i f else (1.5)
0 is the equivalent displacement at the degradation occurrence and 0 + c the critical value of the equivalent displacement at the complete fracture of the element.

Pre-existing discontinuity in the structure

Adhesively joined assemblies and laminate composites contain pre-existing discontinuities. The numerical analysis of such bonded structures at the meso-scale is generally conducted by means of cohesive interfaces/elements whose behavior obeys a TSL generally characterized by a linear ascending branch, possibly a plateau, and then a linear or nonlinear softening branch. Depending on the authors, the criteria for the loss of (initial) linearity (deterioration onset) and ultimate failure (complete de-cohesion) are expressed mostly as a combination of two material parameters, e.g. in terms of critical displacement and critical traction force [START_REF] Lorentz | A mixed interface finite element for cohesive zone models[END_REF] or critical traction force and energy release rate [START_REF] Alfano | On the influence of the shape of the interface law on the application of cohesive-zone models[END_REF] . For simplicity, we are here considering a bi-linear cohesive law as depicted in Fig. 1 . This cohesive law will be used later in Section 4.1 .

The cohesive tractions vector corresponding to Fig. 1 is expressed as: with

T loc = [ C loc ] × δ (1.6)
C loc ( i, j ) = ( 1 -D ) P k ( i, j ) (1.7)
where P k ( i,i ) = E coh , the value of the cohesive law slope. The other components of the matrix are zero, leading to a non-coupling of the normal and shear effects.

The deterioration and complete failure onsets are here defined by critical displacements. The cohesive law accordingly requires three independent constant quantities, viz. E coh , 0 , c .

No pre-existing discontinuity in the structure

In continuous bodies initially exempt of defects, the crack formation is generally preceded by a strain/damage localization possibly induced by microvoid/microcrack coalescence. The progressive deterioration of the mechanical properties of the localization band containing material may be described via a cohesive law. In this case the elastic response is the one of the bulk material [START_REF] Elices | The cohesive zone model: advantages, limitations and challenges[END_REF] , and the cohesive law possesses only a plateau and a softening branch. Once again, criteria are needed to trigger the cohesive law and ultimate fracture. A particular attention must be paid to ensure the stress continuity at the cohesive law onset. For example, the initial (maximum) traction can be a fixed value or a value that is determined in the element at the onset as a function of the current stress state.

In Fig. 2 is depicted an example of a plateau/linear softening cohesive law; as used in this work in the case of no pre-existing discontinuity in Section 4.2 .

The initial traction at localization occurrence in the element is calculated through the use of a stress tensor ˜ σ averaged over all Gauss points of the element as follows

T 0 = ˜ σ n loc (1.8)
T 0 is thus not a constant quantity and is different for each element since it results from its current stress state at the moment of localization. It may also be defined arbitrarily or result from a bifurcation analysis.

The cohesive tractions vector is expressed as:

T loc = ( 1 -D ) T 0 (1.9)
The tangent matrix of the cohesive law in the softening regime is obtained by

C loc ( i, j ) = ∂ T loc ( i ) ∂ δ j = - ∂D ∂ δ j T 0 ( i ) = - 1 c δ j T 0 ( i ) (1.10) with C loc ( i ,1) = 0 if δ 1 < 0
The cohesive law accordingly requires two independent constant quantities, viz. 0 , c .

As soon as the critical value of the de-cohesion variable D is reached the cohesive tractions are automatically set to zero, leading to a traction free crack. The value D c at failure can thus be lower than 1 depending on the material behavior. In the present work it is taken equal to 0.5.

Numerical integration

In the following are outlined the three integration techniques considered in the present work.

Standard Gauss integration

It is reminded that the discrete standard Gauss integration consists in computing [START_REF] Crété | Numerical modelling of crack propagation in ductile materials combining the GTN model and X-FEM[END_REF] where I can be the stiffness matrix or the internal forces vector, nint the number of Gauss points, ξ i the vector containing the local coordinates of the Gauss point, w i the weight of the Gauss point, J the determinant of the Jacobian and f a function of the position of the integration point with respect to the discontinuity.

I = nint i =1 f H ξ i , ξ i , . . . w i J ξ i (1.

Sub-division technique

As was stated in Section 1 the subdivision of hexahedral elements into tetrahedra is the most widely used technique for the numerical integration of elements crossed by discontinuities. After the subdivision, the field variables are projected from the current integration points onto new ones belonging to the tetrahedra, and then a standard Gauss integration rule is applied on each subelement.

Volume averaging integration (VAI)

The volume of the element crossed by the discontinuity is split into two sub-volumes V + and V -separated by the discontinuity, see Fig. 3 left. The proposed volume averaging integration method makes use of a modified quadrature rule that performs the integration twice over the standard Gauss points of the element by means of an averaging of the contributions of the two sub-volumes V + and V -:

I = V - V e nint i =1 f -0 . 5 , ξ i , . . . w i J ξ i + V + V e nint i =1 f +0 . 5 , ξ i , . . . w i J ξ i (1.12)
According to ( 1.12 ), the contribution of a given sub-volume to I consists in a sum over all integration points of the element while assigning the value + 0.5 ( V + ) or -0.5 ( V -) to H in ( 1.11 ) and weighting the result by the ratio of the sub-volume over the total volume of the element. VAI principle is depicted in the visual representation of a cut element with 8 integration points (red marks) in Fig. 3 .

It is noteworthy that this method has been originally proposed by Belytschko [START_REF] Belytschko | A finite element with embedded localization zones[END_REF] and has been very recently applied by Jin et al. [START_REF] Jin | A comparative study of shear band tracking strategies in three-dimensional finite elements with embedded weak discontinuities[END_REF] within an embedded-band based approach. The authors do not, however, address the case where the band crosses the element in such a way that one sub-volume does not contain any integration points, a case that will be precisely addressed in the present work. The VAI technique is applied here in a different context, i.e. to XFEM coupled with a cohesive strong discontinuity and compared to other integration methods.

The proposed scheme must be completed by the use of the B-Bar approach in order to deal with incompressibility. Indeed, in order to alleviate volumetric locking in nearly-incompressible media the B-bar approach is employed (see Hughes [START_REF] Hughes | Generalization of selective integration procedures to anisotropic and nonlinear media[END_REF] ). The matrix B containing the spatial derivatives of the shape functions can be decomposed into a deviatoric and a dilatational component, i.e. B dev and B dil respectively:

B = B dil --- B de v (1.13)
The dilatational part is responsible for an over-stiff response. To palliate this deficiency, the component B dil in ( 1.13 ) is replaced by a new B dil component in a way that volume growth is each integration point is replaced by a value averaged over all points:

B = B dil B de v ⇒ (1.14) B i = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ N i,x + ˜ N i,x ˜ N i,y ˜ N i,z ˜ N i,x N i,y + ˜ N i,y ˜ N i,z ˜ N i,x ˜ N i,y N i,z + ˜ N i,z N i,y N i,x 0 N i,z 0 N i,x 0 N i,z N i,y ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (1.15)
where ˜ N i,x = 1 3

( N i,x -N i,x ) (1.

16)

N i,x is the averaged value of N i,x over all integration points. For the enriched nodes, the B-bar matrix takes the form

B * i = ( H( x ) -H i ) B i .
For more details about the implementation of the Bbar approach in the framework of ABAQUS the reader can refer to Shi et al. [START_REF] Shi | Abaqus implementation of extended finite element method using a level set representation for three-dimensional fatigue crack growth and life predictions[END_REF] .

Application

The methodology coupling XFEM and CZM outlined in Section 2.1 has been implemented as a user element subroutine (UEL) in the commercial finite element computation code ABAQUS. In order to verify the ability of VAI to perform an accurate integration in the case of embedded cohesive strong discontinuities, two different applications are taken into consideration: (i) the single lap adhesive joint in Fig. 4 and (ii) the tensile specimen in Fig. 7 .

In both numerical tests the infinitesimal strain theory is assumed and the tetrahedral elements created during the subdivision stage have one integration point. As shown in 2D in Seabra et al. [21] , by choosing a single point integration rule, more incompressible modes are reproducible, which can alleviate volumetric locking.

In the sequel, the developed 'volume averaging integration' or VAI technique described in Section 3.2 is notably assessed with respect to the subdivision technique above.

Single lap joint test -Pre-existing discontinuity

The single lap specimen, see Fig. 4 , used in this work was taken from Tsai [START_REF] Tsai | An experimental investigation of nonlinear deformations in single-lap joints[END_REF] . Two different finite element models were created, see Fig. 5 . The first model, designated SL.1, contained an explicitly given surface of zero-thickness ABAQUS cohesive elements between the adherents, see Fig. 5 left.

For the second model, designated SL.2, the discontinuity was introduced as the level-set of a cohesive pre-crack crossing the elements in a way that one of the created sub-volumes does not contain any integration points, see Fig. 5 right. This choice is made in order to verify the accuracy of the proposed scheme in the most extreme conditions.

Both models were discretized using a structured threedimensional meshing of full integration hexahedral elements (C3D8 ABAQUS elements) of 1 mm length, therefore for "SL.1 the length of each cohesive element is also of 1 mm.

The material considered for the adherents is a rate and temperature independent mild steel with E = 210 GPa and v = 0.33. The plastic flow condition is defined via the Von Mises yield criterion:

f σ , σ Y = σ eq -σ Y = 0 (1.17)
where σ is the Cauchy stress tensor and σ eq the equivalent stress defined by where s is the deviatoric part of the stress tensor and I the identity tensor. The yield stress σ Y in ( 1.17 ) accounts for isotropic hardening and takes the form of Voce law:

σ eq = 3 2 s : s ; s = σ - 1 3 T r σ I (1.18)
σ Y = R 0 + R ∞ [ 1 -exp ( -k ε p ) ] β (1.19)
where R 0 , R ∞ , k , β are material constants describing the isotropic hardening behavior (see Table 1 ) and ε p the cumulated plastic strain.

The linear cohesive law used in this part of the study is described in Section 2.2.1 and is similar to the cohesive law governing the behavior of the ABAQUS cohesive elements. The properties of the cohesive law are indicated in Table 1 . Since the cohesive crack is inserted in the structure from the very beginning of the analysis, the law is activated within the elements containing the meso-crack from the first time increment.

In the following are conducted four numerical simulations (see Fig. 5 for SL1 and SL2):

• zero-thickness ABAQUS cohesive elements SL.1 • 8-Gauss point standard integration technique SL.2 • sub-division integration technique SL.2 • volume averaging integration (VAI) technique SL.2 SL.1 is used here as a reference solution since ABAQUS cohesive elements are a tool of known accuracy for crack propagation in interfaces.

The numerical results for the four numerical simulations are superimposed in Fig. 6 in terms of reaction force vs. displacement.

According to Fig. 6 , the standard Gauss quadrature (8 GPs) is unable to provide a correct prediction of the reaction force (leading to 0 N all along the loading), as expected. On the other hand, the VAI and subdivision schemes give the same response as the cohesive zone model.

Tensile test -No pre-existing discontinuity

For the tensile test the specimen of Fig. 7 is used. It is discretized using a structured three-dimensional meshing of complete integration hexahedral elements (C3D8 ABAQUS elements) of 0.5 × 0.5 × 0.5 mm 3 (in the area of interest).

The material taken into consideration here is the same as the one of the adherents in Section 4.1 , see ( 1.17 -1. [START_REF] Jin | A comparative study of shear band tracking strategies in three-dimensional finite elements with embedded weak discontinuities[END_REF] ) and Table 1 .

The critical value to be attained by an element for localization (activation of XFEM enriched degrees of freedom) is tentatively an arbitrary critical cumulated strain ε p c , see Table 1 , giving the value of T 0 in ( 1.8 ). The linear/plateau cohesive law presented in Section 2.2.2 is then used to describe the meso-crack deterioration effect.

In the following are conducted three numerical simulations:

• elasto-plastic response (FEM only)

• sub-division integration technique • volume averaging integration (VAI) technique

For the last configurations involving the formation and development of a cohesive strong discontinuity, three through-thickness orientations (with respect to the loading direction) of the discontinuity are considered, i.e. 0 °(brittle-like fracture), 54.7 °(ductile-like fracture [START_REF] Pineau | 2.06-failure of metals[END_REF] ) and 45 °(worst numerical case). For the 0 °case, the crack is forced to appear extremely close to the localized elements' side, ensuring again that one of the created sub-volumes does not contain any integration points. For the two other cases, the crack position within the elements is not imposed since the number of cut elements ensures the existence of multiple sub-volumes that will satisfy the above mentioned extreme condition. The angle of 45 °is the most critical numerical case in the sense that, for regular meshing, it guarantees that one out of two elements will not have Gauss points on one side of the crack. The 54.7 °case does not necessarily imply the aforementioned problem but it is addressed for physical reasons. The fractured specimens are shown in Fig. 8 .

The two integration procedures are detailed below.

• Volume averaging integration (VAI) technique

The localization band propagates throughout the whole width of the specimen (with the chosen angle) when one element attains the critical cumulated strain in one increment. This increment denoted kinc LOC is saved for later use.

• Sub-division integration technique

A band is defined at the very beginning of the analysis with the same path as the one created considering the previous integration technique, resulting in a subdivision in tetrahedra (whose number depends on the way the element is cut into the two sub-volumes) and a distribution of Gauss points from the first increment of the analysis, although the enriched dofs of the localized elements are activated only when the analysis reaches kinc LOC . Indeed, in order to make a successful and pertinent comparison between VAI and the sub-division schemes it is important to control the localization in terms of (i) time of occurrence, (ii) path and (iii) tetrahedra created. It is noted here that in order to calculate the sub-volume values in VAI a triangulation technique is used, dividing the volume into tetrahedra whose volumes can be easily calculated and then summed (these tetrahedra have no other participation in the integration procedure). For both cases, by forcing the band to appear throughout the whole width of the specimen in one increment (kinc LOC ) and on the same path, it is possible to control (i), (ii) and (iii) in a way that they are exactly the same for the two schemes. Having isolated these three factors we can compare the schemes solely on their ability to integrate the system of equilibrium equations.

The numerical results for the three numerical simulations are superimposed in Figs. 9 , 10 and 11 for the three discontinuity orientations in terms of reaction force vs. displacement. According to Figs. 9 , 10 and 11 , VAI and subdivision techniques give similar responses. A very slight divergence can be seen in the cases of 45 °and 54.7 °, but this can be attributed to a difference in the last time step of the analysis.

It should be noted that the standard 8 Gauss point scheme is not presented in any of the two cases because when it was used no analysis was capable to converge after the onset of the localization, so it is omitted.

Discussion

It is reminded that the subdivision of elements consists in triangulating the sub-volumes V + and V -, separated by the discontinuity in several tetrahedra and projecting the state variables known at the current integration points onto the newly formed integration points of the tetrahedra. These operations are costly in terms of computation time and potential sources of numerical errors.

The volume averaging integration (VAI) technique does need simply the values of the sub-volumes V + and V -. The cost in this case is very low and even negligible.

The numerical examples treated in Sections 4.1 and 4.2 have evidenced that both techniques lead to the same results, demonstrating the interest of using the volume averaging integration technique.

It is noteworthy that for tension loading the softening part of the CZM law gives the shape of the overall response in the softening regime. The apparent similarity in the brittle and ductile-like responses in Section 4.2 is due to the fact that the same law was utilized for the sake of simplicity and must not be considered as a generality.

Concluding remarks

A volume averaging integration technique or VAI is proposed for the integration of the equilibrium equations when using the XFEM coupled with cohesive zone model. Two different numerical tests, i.e. a single lap joint test and a tensile test, have been realized in order to verify the accuracy of the proposed scheme. It has been shown that the proposed VAI scheme leads to the same level of accuracy as the well-known subdivision technique, while being much simpler to implement and less costly to use.

Fig. 1 .

 1 Fig. 1. Bi-linear cohesive law.

Fig. 2 .

 2 Fig. 2. Linear/plateau cohesive law.

Fig. 3 .

 3 Fig. 3. Visual representation of the "volume averaging" integration scheme.

Fig. 4 .

 4 Fig. 4. Single lap specimen dimensions [22] .

Fig. 5 .

 5 Fig. 5. Visualization of the two different single lap models.

Fig. 6 .

 6 Fig. 6. Single lap test reaction force graph.

Fig. 7 .

 7 Fig. 7. Traction specimen dimensions.

Fig. 8 .

 8 Fig. 8. Tensile specimen cut by 0 °, 45 °and 54.7 °cracks (cumulated plastic strain distribution).

Fig. 9 .

 9 Fig. 9. Tensile specimen reaction force graph (0 °case).

Fig. 10 .

 10 Fig. 10. Tensile specimen reaction force graph (45 °case).

Fig. 11 .

 11 Fig. 11. Tensile specimen reaction force graph (54.7 °case).

Table 1

 1 Hardening law & cohesive laws properties.

Acknowledgement

The authors would like to acknowledge the French Direction Générale de l'Armement (DGA) for its support.

Declaration of Competing Interests

None.


Joseph Marae Djouda
email: joseph.marae_djouda@epf.fr

Donato Gallittelli
email: donato.gallitelli@gmail.com

Marouene Zouaoui
email: marouene.zouaoui@utt.fr

Ali Makke

Julien Gardan
email: julien.gardan@epf.fr

Naman Recho
email: naman.recho@epf.fr

Jérôme Crépin
email: jerome.crepin@mines-paristech.fr

Local scale fracture characterization of an advanced structured material manufactured by fused deposition modeling in 3D printing

Keywords: Additive manufacturing, 3D printing, Digital Image Correlation, Micro speckle, Fracture mechanics

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

dditive manufacturing (AM) is a promising principle able to shape complex structures through a layer by layer manufacturing. In recent years, AM has been the subject of several developments both by academic and/or industrial researchers. AM has great potential for the manufacture of several structures that can be integrated into the industrial sectors such as aeronautics, automotive and others [START_REF] Gardan | Additive manufacturing technologies: State of the art and trends[END_REF]. To benefit from its advantages, it is important to characterize accurately the obtained structures to allow a better integration in industrial applications. In the specific case of fused deposition modeling (FDM) technologies, the mechanical behavior of the obtained structures suffers from anisotropy, which is introduced by the weld lines and the air gap in the printed part [START_REF] Ahn | Anisotropic material properties of fused deposition modeling[END_REF]. Many studies aiming to improve the mechanical properties of these structures have been conducted. Some of these studies focus on the modification of the shape of the printed patterns [START_REF] Prall | chiral honeycombe. Poisson's ratio. 1. INTRODUCTION Cellular solids are used widely in a variety of engineering applications. In particular[END_REF] and others are based on improving 3D printing algorithms to enhance the local fracture toughness of printed thermoplastics [START_REF] Gardan | Improving the fracture toughness of 3D printed thermoplastic polymers by fused deposition modeling[END_REF]. Besides of all these processes of pre-structured material depositions, accurate characterization of the mechanical behaviour of the obtained structure is determinant. The knowledge of their mechanical properties is important for use in the technological application. The mechanical characterization of a material goes through a knowledge of its properties both at macroscopic and local scales. Many technics have been developed to analyze strain evolutions at the local scale. These technics are based on the used periodic or random gratings and have been applied in most of the case in metal specimens [START_REF] Allais | Experimental characterization of the local strain field in a heterogeneous elastoplastic material[END_REF]- [START_REF] Clair | Strain mapping near a triple junction in strained Ni-based alloy using EBSD and biaxial nanogauges[END_REF]. The periodic gratings consist in arrays of cross lines or of nanoparticles and are deposited by lithographic technics (electron beam lithography, holographic or photolithography). The periodic gratings can be coupled with the technics like Moiré Interferometry [START_REF] Du | Large-deformation analysis in microscopic area using micro-moiré methods with a focused ion beam I milling grating[END_REF] or with a specific approach like the one developed in our previous studies [START_REF] Djouda | Investigation of Nanoscale Strains at the Austenitic stainless steel 316L surface : Coupling between Nanogauges Gratings and EBSD Technique during in situ Tensile Test[END_REF], [START_REF] Maurer | The beginnings of plasmomechanics: towards plasmonic strain sensors[END_REF]. The technics very often associated with the random gratings are Moiré interferometry and digital image correlation (DIC) [START_REF] Tanaka | Development of a pattern to measure multiscale deformation and strain distribution via in situ FE-SEM observations[END_REF]. In these technics, the roughness of the samples can be used in certain cases. In other cases, paint speckle pattern, deposition and annealing of thin metal layer in order to obtain metal nanoparticles, and lithography technics can be used to deposit gratings [START_REF] Joo | Method for measuring nanoscale local strain in a dual phase steel using digital image correlation with nanodot patterns[END_REF]. In a very recent study, DIC has been used to measure strains at the surface of 3D printed structures with a large dimension [START_REF] Schnittker | Integrating digital image correlation in mechanical testing for the materials characterization of big area additive manufacturing feedstock[END_REF]. The DIC analyses the grey levels between two images, one of which is considered as "reference" and the other is "deformed". In this technic, the images used can be obtained by various imaging systems: optical or digital microscopes, scanning electron microscope, atomic force microscope, and others. In addition to the surface kinematic fields, DIC could allow to access the out-of-plane displacements. In the basic version of DIC, the roughness of the surface can be used like a speckle pattern but sometime the contrast of images is not enough. Simple surface preparation may be appropriate depending on material (polymers or alloys) and imaging technique used [START_REF] Allais | Experimental characterization of the local strain field in a heterogeneous elastoplastic material[END_REF], [START_REF] Joo | Method for measuring nanoscale local strain in a dual phase steel using digital image correlation with nanodot patterns[END_REF]. The spatial resolution of strains measured by technics using gratings is strongly related to the grating parameters: the size of particles and pitch. With the progress in nanotechnology and nanofabrication, it is possible to deposit gratings with dozens of nanometre of pitch and diameter of nanoparticles. However, it is difficult to use lithographic techniques to deposit gratings on the surface of a polymer structure, because of the solutions used during the nanofabrication process. Moreover, it could be straightforward to image polymer materials by SEM. The key challenge for better understanding mechanisms leading to structure damages or failures remains the development of tools giving direct access to the local scale material behaviour. For a good characterization of polymer structures obtained by AM, it is necessary to develop new techniques to measure local deformations.

A

In this study, an original process is developed to analyze strain evolutions at the surface of polymers structures obtained by AM. It consists in deposing speckle pattern with micron dimension. An in situ tensile test under a digital microscope in which images of the specimen surface are acquired is conducted until the specimen failure. The images obtained are analyzed with DIC tools. The specimen consists in a mini SENT. The deformation evolutions are analyzed for different loading levels and local mechanisms of deformation are highlighted. The local strain evolutions are confronted with the macroscopic behavior of the specimen.

METHODS AND RESULTS

n this section, the detail of specimens' preparation and the implementation of tests are presented. Moreover, the main results will be discussed.

a) b) c) d)

Specimen preparation

The specimen is printed by adding melted layers of Acrylonitrile butadiene styrene (ABS) using a makerbot replicator 2X. The layers were deposed at a temperature of 235 °C with a +45°/-45° oriented filaments as described in Fig. 1a). A layer thickness of 0.25 mm was used except for the first one (0.3 mm) in order to ensure good adhesion with the manufacturing tray. A mini SENT specimen is 3D printed with the dimension showed in Fig. 1a). The notch was also made by AM. However, because of the temperature during the manufacturing, the layers have closed at the notch and the latter is not clearly visible. The yellow arrow in Fig. 1-b) illustrates its position on the specimen surface. The notch was made to meet the condition 𝑎 𝑤 ⁄ 0.45 as recommend in ASTM E1820 standard; 𝑎 represents the initial length of the notch and 𝑤 the width at the useful part of specimen. It is important to notice that the specimen like design herein is not standardized. The initial thickness of the specimen 3D printed was of 6 mm and by mechanical polishing each side, the final thickness considered was 3 mm in order to reduce the influence of surface defects on the mechanical feature. The final surface of the specimen and also to obtain a soft surface in order to deposit the speckle pattern to allow the use of the DIC. Some defects (hole, void ...) are visible on the specimen surface.

Speckel characterization

DIC consists in a segmentation of an image into subsets. The changes in the different subsets on the successive images are followed. The dimensions of the speckles and their number in each subset are determinant to obtain a good resolution. In this study, an original speckle pattern with micron dimensions is deposited. In contrary of the conventional process of speckle pattern deposition with the white base coat, in the present study the surface of the sample (red) was used as the base. White speckles of micrometric dimensions were deposited using an airbrush. This results in a speckle pattern with the particles of 20 microns of average diameter (see Fig. 1-c)).

In situ tensile test

Once the speckle grating deposited at the surface of the sample, an in situ tensile test under a digital microscope was conducted. The experimental set up consists in a numerical microscope Keyence VHX-1000 made by Keyence Corporation TSE for the surface observation, a tensile micro-machine and a trigger system (see Fig. 1-d)). The trigger system allows recording of images at a specific rate when the specimen is continuously solicited during a tensile test and related each image to the corresponding applied load in the macroscopic behavior of the specimen. One image was recorded per second during the test until the failure of the specimen.

DISCUSSION OF THE RESULTS

he process developed herein allows obtaining both macroscopic and local evolutions of the deformations of the material. Fig. 2 shows the strain -load curve. At the elastic-plastic transition (around 3.6% of macroscopic deformation), an important decrease of the applied load is observed, and this materializes the sudden opening of the notch. This last is initially enclosed due to the deposition temperature of neighboring filaments during the deposition process. After this load variation, the applied load increases again, reach to one threshold, and decrease until the failure of the specimen. Some images are shown in Fig. 2 in order to illustrate the main local features at the specimen surface. The applied load leads to a gradual deformation at the tip of the crack. Once the crack reaches the design length 𝑎 , the concentric curved filaments at the crack tip begin to experience deformation. The deformation at the crack tip corresponds to the strain lobe visible in Fig. 2 -Step 725). Another's strain concentration zones (zones (2) and ( 3)) are also visible. The zone (1) corresponds to the crack tip, zones (2) and ( 3) are equidistant from the perpendicular to the direction of traction passing through the initial crack direction. The zones (2) and (3) could correspond to the boundaries between the two concentric curved filaments materializing the crack tip and layers oriented +45°/-45° of the specimen (see complementary material). However, it is important to notice that the positions of strain initiation at the surface of the specimen do not correspond to visible defects (hole, void). They could correspond to defects in the volume of the material, to filling defects between curved filaments and oriented layers +45°/-45°. Once the crack has spread and is beyond the two curved filaments, layers oriented +45°/-45° begins to deform. A V-shape notch tip is observed at the crack tip and the deformation of the layers oriented +45°/-45° highlighted tears marks oriented +45°. The strain marks between the crack tip and zone (2) increase gradually until they overlap (See Fig. 2 -Step 780). With the applied load, zone (1) and zone (2) form now a unique zone of strain concentration, the strain at the zone (3) have also greatly increased and the cavities developed in the volume of the material are visible at the surface. Ncorr software was used for the strain quantitative analysis. Ncorr is an open-source 2D-DIC MATLAB Software developed at Georgia Institute of Technology. Furthermore, this software integrated the modern algorithm of DIC [START_REF] Blaber | Ncorr: Open-Source 2D Digital Image Correlation Matlab Software[END_REF]. Fig. 3 presents the strain maps corresponding to images from the previous section. The strain 𝜀 was computed using the Green-Lagrangian formalism in which the four-displacement gradient is used in order to reduce noise from differentiation. The strain concentration and magnification are in good agreement with the deformation evolutions observed. At step 551, no strain concentration is visible at the surface of the specimen. This step is macroscopically located in the elastic domain of the material. The plastic domain is accompanied by the sudden opening of the notch. The strain map of step 582 shows clearly strain concentrations at the crack tip. These strains could correspond to the beginning of deformation of the curved filaments as observed in the images. Another strain concentration zone (zone (2)) is visible. The strains are quite uniform in the rest of the specimen. By increasing the applied load, strain at the crack tip increases in dimension and in magnification (Fig. 3 -Step 725). Locally strain is around 25% and macroscopic strain value is less than 5%. Then local strain can reach very significant values compared to macroscopic one. Similar observations were made in the case of austenitic stainless steel 316L by applying the nanogauges technic [START_REF] Djouda | Nanogauges gratings for strain determination at nanoscale[END_REF]. It is important to notice that strain at the crack tip grows in the direction of the zone (2). Another T observation is the strain concentration materializing the zone (3) like indicated in the previous section. The beginning of this zone below to the crack could also correspond to the frontier between the curved filaments and filaments oriented 45°. At step 780, strain at the crack has increased so that they start to overlap with those of the zone [START_REF] Ahn | Anisotropic material properties of fused deposition modeling[END_REF]. The strain at the zone (3) have also increased and at the different strain concentration area, the strain magnification is up than 40%. Step 780 is macroscopically at the beginning of one plateau. This last could express the competition between cavities growth and necking. At step 885, the crack tip and zone (2) overlap and form a unique strain concentration zone. This unique zone contains the crack notch and with the zone (3), they now form the two mains strain concentration zones. They are all oriented almost perpendicular to the tensile direction and the crack tip seems to interact with the zone (3). The local strain intensity reaches the values up to 100%. Macroscopically, this step is in the necking domain, a few moments before the failure of the specimen.

CONCLUSIONS

n this study, an original approach is developed for local strain characterization at the surface of polymer materials produced by Fused Depositing Modeling (FDM) in 3D printing. This approach allows a comprehensive qualitative analysis of the local deformation from the optical image recorded during the in-situ tensile test. The micro speckle pattern deposited at the surface of the specimen ensures the strain quantification from a DIC software. Ncorr 2D-DIC MATLAB Software was used. The strain evolutions in the strain maps are in good agreement with the strain evolutions observed in the optical images. The surface strain evolutions are quantified and the interactions between strain concentration zones are clearly evidenced. These strain concentration zones could correspond to filling defects between curved filaments and oriented layers +45°/-45°. The local strain values are very important compared to the macroscopic ones and the gap increases with the applied load. The approach developed herein presents great potential for local characterization of materials in general and for polymer materials specifically because of the limitations they present in regard of the classical approach. It also brings precious details which could help for the optimization of smart materials or pre-structured materials by additive manufacturing.

Figure 1 :

 1 Figure 1: a) Plan of the sample used in the study. It consists to SENT specimen with a notch set as 𝑎 𝑤 ⁄ 0.45. b) The sample is 3D printed layer by layer and oriented +45°/-45°. After surface preparation of the sample, the micro-speckle pattern is deposited. c) Histogram shows the diameter distribution of speckles at the sample surface. More than 70% of speckles have their diameter less than 20 µm. d) The experimental set up consist in a digital microscope for the surface observation, a tensile micromachine and trigger system.

IFigure 2 :

 2 Figure 2: Strain evolutions from the optical images. The image at step 551 is recorded at the magnification x50, the notch is not visible at this loading level. The applied load leads to the sudden opening of the notch. The images recorded after the notch opening is at the magnification x100 to allow a better visualisation of the local strain evolutions. The scale bar is 1 mm in all the images.

Figure 3 :

 3 Figure 3: Quantitative analysis of strain evolutions. Step 551 shows the region where the strains are analyzed in front of the design notch materialized by the red arrow. The scale bar at each image is 1 mm.

COMPLEMENTARY MATERIAL