R Fedele
email: fedele@stru.polimi.it

† B Raka

F Hild

S Roux

Identification of adhesive properties in GLARE assemblies using Digital Image Correlation

Keywords: Cohesive crack model, Digital Image Correlation, Inverse Problems, Parameter identification, Structural adhesives

An experimental-numerical methodology is introduced to identify the parameters of a cohesive law of an adhesive layer within a joined assembly on the basis of kinematic data provided by digital image correlation. Nonconventional experiments on joined samples were designed to generate within the assembly and the adhesive film complex strain and stress states close to those expected in service and up to complete debonding. The modeling is developed with reference to the observed sub-domain in which the experimental boundary conditions are prescribed. The non-linear behavior of the adhesive layer is described as a finite thickness interface endowed with a mixed-mode cohesive law whose parameters are identified so as to match at best the measured displacement field.

Introduction

Joining two materials to obtain a better performance as that of the singles is a general goal in composite design. Various architectures are obtained when dealing with layers or fibers. Among these, GLARE has received some attention in the aerospace industry [START_REF] Vlot | Glare, History of the Development of a New Aircraft Material[END_REF][START_REF] Sinke | Manufacturing of GLARE Parts and Structures[END_REF]. GLARE is a "GLAss-REinforced" fiber metal laminate, composed of several layers of aluminum alloy interspersed with layers of glass-fiber "pre-preg." In the design of stiffened panels for aircraft structures, recourse has also been made to adhesively-bonded assemblies between GLARE laminates with different geometries. To predict the structural response of such assemblies, a good understanding of the adhesive layer properties is required. In this paper, it is proposed to use images of the joined assembly in a non-conventional mechanical test, representative of in-service conditions, to identify the parameters governing a mixed-mode cohesive model for the adhesive layer.

When dealing with model calibration of adhesive joints, point data, e.g., strain and load data, are usually the only experimental information available [START_REF] Derewonko | Strength assessment of adhesive-bonded joints[END_REF]. Pictures at different scales are used in a qualitative way in addition to global data [START_REF] Pardoen | Constraint effects in adhesive joint fracture[END_REF][START_REF] Leffler | Shear behaviour of adhesive layers[END_REF][START_REF] Sørensen | Micromechanical model of cross-over fibre bridging Prediction of mixed mode bridging laws[END_REF][START_REF] Salomonsson | Modeling and parameter calibration of an adhesive layer at the meso level[END_REF], or quantitatively by evaluating deflections [START_REF] Yang | Elastic-plastic mode-II fracture of adhesive joints[END_REF][START_REF] Su | An elasticplastic interface constitutive model: application to adhesive joints[END_REF] and deformed shapes [START_REF] Sargent | Durability studies for aerospace applications using peel and wedge tests[END_REF]. In this study, it is proposed to use full-field measurements provided by Digital Image Correlation (DIC). Its advantage lies in the fact that numerous measurement points are available to analyze an experiment, as opposed to the standard practice of using few data.

During the last years, DIC [START_REF] Sutton | Advances in Two-Dimensional and Three-Dimensional Computer Vision[END_REF] has undergone substantial developments for several reasons. First, its application is generally simple and straightforward. It works effectively under natural light. The resolution is of the same order as many conventional measurement devices. It is employable at different scales, from macro down to micro or even nano scales [START_REF] Sutton | Local crack closure measurements: Development of a measurement system using computer vision and a far-field microscope[END_REF][START_REF] Soppa | Experimental and numerical characterisation of in-plane deformation in two-phase materials[END_REF][START_REF] Chasiotis | A New Microtensile Tester for the Study of MEMS Materials with the Aid of Atomic Force Microscopy[END_REF][START_REF] Forquin | A Method to Determine the Toughness Scatter of Brittle Materials[END_REF][START_REF] Chasiotis | Special issue on nanoscale measurements in Mechanics[END_REF]. In all the previous cases, DIC was based upon local matching of interrogation windows in a series of pictures. An alternative approach, based on a finite-element discretization of the displacement field [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF], will be used herein.

Several attempts have been made to enrich the experimental information gathered by a conventional equipment with full-field measurements [START_REF] Mcneill | Estimation of stress intensity factor by digital image correlation[END_REF][START_REF] Geers | Mixed numerical-experimental identification of non-local characteristics of random-fibre-reinforced composites[END_REF][START_REF] Abanto-Bueno | Investigation of crack growth in functionally graded materials using digital image correlation[END_REF][START_REF] Cho | Measurement of Nanodisplacements and Elastic Properties of MEMS via the Microscopic Hole Method[END_REF][START_REF] Maier | Some innovative industrial prospects centered on inverse analyses[END_REF][START_REF] Hild | Digital image correlation: From measurement to identification of elastic properties -A review[END_REF]. However, very few studies deal with the identification of cohesive models because of experimental challenges and the scales at which data are to be obtained. Abanto-Bueno and Lambros (2005) used a multi-camera system to determine the traction separation law.

The near-field images were used to measure the opening displacements, and the far-field views to evaluate the interface tractions via an elastic solution.

The above procedure, however, is not a reliable route to follow when the adherend behavior is nonlinear and only microscopic measurements are used (as in the present case).

Advanced mechanical models for describing the nonlinear behavior of materials and structures contain numerous parameters, many of which are not easy to identify. Therefore, there is strong interest in industrial environments to develop procedures for estimating these parameters on the basis of innovative experimental tests [START_REF] Yoshida | Inverse approach to identification of material parameters of cyclic elasto-plasticity for component layers of a bimetallic sheet[END_REF]. The latters are designed with the aim of generating stress and strain states close to those expected in-service. This approach leads to difficult challenges both in terms of measurement and identification. On the measurement side, stress and strains are typically very heterogeneous, possibly at the microscopic scale, or difficult to access. For identification, the usual occurrence of different non-linearities requires a full modeling that often makes the procedures quite difficult. The novel developments, which open the way to such delicate task, is the original combination of experimental full field measurement and inverse problem approaches.

The above (inverse) procedure is applied to the characterization of a structural adhesive layer, joining a laminate skin and a Z-shaped reinforcement (stringer), both made of GLARE for aerospace applications. The experiments, specially designed for the present investigation, were monitored by a digital camera, empowered by a long distance microscope (Section 2). A meaningful sample of pictures was analyzed by a finite element based DIC algorithm. The model adopted herein to describe the response of the monitored Region-of-Interest (ROI) is presented in Section 3. The mechanical properties of the adhesive layer, modeled as a cohesive thick interface law, is estimated via an inverse procedure, and the traction profiles along the adhesive layer reconstructed (Section 4).

Notations. Small strain and displacements are assumed. The symbol : denotes the contraction with respect to two indices, namely, A :

B = A ij B ij .
Bold symbols denote vectors, matrices, or tensors, depending on the context.

Debonding experiments

Non-conventional tests are designed to generate complex stress and strain states within a joined assembly. Under these conditions, nothing impedes distinct non-linear mechanical phenomena from occurring during the test, and interact with each other (e.g., fracture of the adhesive layer, simultaneous plastic yielding of the adherends, or local delamination). Therefore, numerical simulations of the overall specimen response should take into account all these phenomena simultaneously. In this context, displacement measurements, say by DIC, turn out to be advantageous when compared with traditional "point" and overall measurements. In particular, DIC allows for the evaluation of the actual kinematic boundary conditions of the modeled domain. The experiments reported herein are applied to the characterization of a structural (60-µm thick) adhesive layer, joining a (1.4-mm thick) laminate skin and a Z-shaped reinforcement (stringer), both made of GLARE (Figure 1). The single epoxy adhesive layer, shown in Figure 2, is made of BMS 5-101 produced by 3M (commercial denomination AF163 degree 10).

Tests were designed and performed on rectangular joined specimens (surface area: 297 × 50 mm 2 ) to lead to complete debonding of the adhesive (Fig-

ure 3). The experiments were monitored by a digital camera equipped with a long distance microscope. A small region on the specimen surface around the adhesive layer (of area approximately equal to 2 × 2.5 mm 2 ) was selected as the ROI and monitored during the tests. The pictures at the microscale (Figure 4) were subsequently processed to determine the displacements ex- The code used herein is based on a Galerkin, finite-element discretization of the displacement field with bilinear shape functions [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF]. It allows one to establish a close connection with subsequent finite element simulations. Figure 5 pictures of Figure 4). The adopted finite element (mechanical) model to be introduced in the sequel is consistent with such measurements.

Mechanical model

The chosen mechanical modeling is matched to the same mechanical problem as in the experiments with the exclusive reference to the considered ROI. Displacement fields measured by means of a DIC procedure are suitable to specify Dirichlet boundary conditions for such a region. The debonding pro- of the debonding process. It represents an important simplification from a numerical standpoint, both in the forward simulations and in the analytical sensitivity analysis.

Aluminum alloy substrate

The aluminum (cold-rolled 2024-T3) alloy substrates are modeled according to Prandtl-Reuss' flow rule [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF])

= 1 + ν E S + 1 -2ν E pI + 3 2 α q σ 0 n-1 S ( 1
)
where I denotes the second-order identity tensor, the (infinitesimal) strain tensor, S the stress deviator tensor, q = 3 2 S : S its normalized secondorder invariant (von Mises' equivalent stress), p = 1 3 σ : I the hydrostatic stress, E and ν the isotropic elastic constants. Parameters α, σ 0 and n govern the (fully incompressible) nonlinear part of the material response. All parameters describing the response of the aluminum alloy sheets are assumed to be known a priori with sufficient accuracy.

Cohesive layer

To describe the joint response under mixed-mode loading conditions, a modified version of the Xu-Needleman exponential law [START_REF] Xu | Void nucleation by inclusions debonding in a crystal matrix[END_REF], rooted in the atomistic description of binding energies, is chosen. This law, formulated for a zero-thickness interface, specifies a holonomic relationship between the tractions T on the interface, and the corresponding displacement discontinuity vector ∆. Normal and tangential components of traction and displacement vectors are denoted by T n and T t , ∆ n and ∆ t , respectively. Symbols φ n and φ t designate the works-of-separation (for a unit surface) under pure Mode I (peel) and Mode II (shear), respectively.

The (non-dimensional) ratio between the tangential and normal works-ofseparation is denoted by r = φ t /φ n . Originally, the Xu-Needleman interface response assumes that the works-of-separation φ n and φ t are equal (r = 1).

However, in the case of structural adhesives the ratio r generally lies in the range 3 ≤ r ≤ 6 [START_REF] Valoroso | A damage-mechanics-based approach for modelling decohesion in adhesively bonded assemblies[END_REF]. [START_REF] Van Den Bosch | An improved description of the exponential Xu and Needleman cohesive zone law for mixed-mode decohesion[END_REF] proposed an extension of the Xu-Needleman interface law that allows for different normal and tangential works-of-separation

T n = φ n δ n ∆ n δ n exp - ∆ n δ n exp - ∆ 2 t δ 2 t (2a) T t = 2φ t δ t ∆ t δ t 1 + ∆ n δ n exp - ∆ n δ n exp - ∆ 2 t δ 2 t (2b)
The above cohesive relationship is governed by four parameters, namely, φ n , φ t , δ n and δ t that describe the works-of-separation (per unit surface) and the characteristic lengths under mode I and mode II conditions, respectively.

Although a sequence of images will be used for which the mode mixity is not constant, the experienced range of mixity is too narrow to lead to a robust determination of these four parameters as a whole. Complementary results from fracture tests under different loading conditions should be considered (e.g., pure mode I or mode II data). However, if the ratio r is prescribed as an a priori known information, then the problem turns out to be sufficiently well posed for a robust determination of the three remaining free parameters (Section 3.5). Due to the significant scatter usually observed in the mechanical response of adhesive joints, this externally prescribed value is considered as a viable route to follow. In the following, the three unknown parameters governing the adhesive layer behavior are gathered in a vector form, X = {φ n , δ n , δ t } T .

Extension to finite-thickness interfaces

Zero-thickness representation of cohesive interfaces is widely adopted in computational mechanics to simulate quasi-brittle fracture processes at dif-

t n ∆ t ∆ S + S - layer thickness 1 - 1 + 2 - 2 + n T t T n T t T Figure 6:
Interface element with a finite thickness.

ferent scales [START_REF] Bolzon | Parameter identification of a cohesive crack model by Kalman Filter[END_REF][START_REF] Shet | Effect of the Shape of T δ Cohesive Zone Curves on the Fracture Response[END_REF][START_REF] Andena | Cohesive zone modelling of fracture in polybutene[END_REF]. However, this simplification cannot be made in the present application. First, from a geometrical point of view, the thickness of the adhesive layer (approximately 60 µm) is not negligible with respect to the overall size of the modeled region at the microscale (about 600 µm for the shortest size of the ROI). Second, elastic strains are expected to influence the overall response of the considered sub-domain. While the crack propagates, a part of the adhesive joint is compressed (a slightly reduced layer thickness can be noticed on the left side of the ROI in Figure 4).

In this study, the adhesive layer is modeled by means of finite-thickness interface elements. A simplified kinematics is assumed by describing two dominating deformation modes, namely, peel (normal) and shear (tangen-tial), both constant through the thickness (Figure 6). From a mechanical standpoint, corroborated by experimental evidences, an adhesive layer is expected to behave differently from the adhesive as a bulk material [START_REF] Högberg | Mixed mode cohesive law[END_REF]. The constitutive law of the adhesive layer, even with finite thickness, is thus expressed to be a function of normal, ∆ n , and tangential, ∆ t , displacements. The cohesive element adopted herein can be regarded as an interface line element, with an initial gap between opposite nodes, otherwise coincident.

With such an interface element, let us consider a suitable partition for the nodal displacement vector U , namely U = U -, U + T , where superscripts (•) ± denote corresponding nodes belonging to the upper and lower element side, respectively (Figure 6). The relative displacement vector is expressed as

∆ =    ∆ n ∆ t    = [-R T N U -, R T N U + ] = [-B s U -, B s U + ] = B U (3)
where B is the compatibility operator, and R the orthogonal rotation matrix,

int , and the tangent stiffness matrices, K (e) u , are defined as [START_REF] Goyal | Irreversible constitutive law for modeling the delamination process using interfacial surface discontinuities[END_REF])

F (e) int =    -F s F s    F s = +1 -1 B T s T |J| d ξ (4) K (e) u =   K s -K s -K s K s   K s = +1 -1 B T s K u B s |J| d ξ (5)
where ξ is the iso-parametric line coordinate, |J| the jacobian (here constant)

of the iso-parametric (line) mapping, K u = ∂T /∂∆ T the material tangent stiffness matrix. One-dimensional integrals in Equations ( 4) and ( 5) are computed via a Newton-Cotes scheme (with collocation points at the two ends of the master element midline, for ξ = ±1). The formulation (3-5) allows one to obtain an interface response that is independent of the layer thickness, thus avoiding the scaling of work-of-separation, which is instead required when the dissipated energy is smeared over the element area [START_REF] Balzani | An interface element for the simulation of delamination in unidirectional fiber-reinforced composite laminates[END_REF].

The cohesive law (2a-2b) is further modified to take into account compressive normal tractions (T n < 0) generated within the finite-thickness elements. The situation predicted by the cohesive relationship for ∆ n → 0 + , Equations (2a)-(2b), is prolonged when ∆ n < 0 as follows: (i) along the normal direction n, the adhesive response is governed by the tangent stiffness matrix at the origin, as if it were linearly elastic, in a decoupled way with respect to the tangential response; (ii) along the tangential direction t, a zero value of ∆ n is assumed in Equation (2b) and the tangential response T t is the one resulting from a pure mode II. Thus, in the present study, Equations (2a)-(2b) specify the adhesive response for ∆ n ≥ 0, whereas for ∆ n < 0 the following de-coupled equations hold

T n = φ n δ 2 n ∆ n (6a) T t = 2φ t δ t ∆ t δ t exp - ∆ 2 t δ 2 t (6b)
This a priori simplification of the adopted model for ∆ n < 0 assumes no debonding in compression within the adhesive layer, friction or other kinds of coupling between tangential and normal tractions.

Results

n T

[MPa]

10 i = 7 i = 13 i = crack front
joint length [mm] 500 elements 2,000 elements The plane model adopted to simulate the mechanical response of the considered sub-domain is made of 570 finite elements, and the adhesive layer is discretized by means of 39 cohesive elements. Plane stress conditions are assumed since kinematic measurements provided by DIC concern the free surface of the specimen. The computer code for the nonlinear FE analyses and computation of parameter sensitivities was entirely developed in a Matlab environment [START_REF] Bonnet | Analyse des solides déformables par la méthode des éléments finis[END_REF], fully consistent with the DIC software in terms of mesh and shape functions.

It is worth noting that because of measurement uncertainties, the element size cannot be decreased at will since it will be accompanied by higher measurement uncertainties [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF]. In the present case, an element size equal to 16 pixels (or 30.9 µm) is a compromise between uncertainty level and spatial resolution of the DIC procedure. To evaluate the mechanical effect of what may appear as a rather coarse mesh (i.e., 500 elements), a comparison is carried out with a finer mesh (using 2000 elements). layer X and the displacements at measurement instant i inside the ROI, adopted as measurable quantities, namely u comp Ω i , is usually referred to as forward or direct operator (as opposed to inverse), and denoted by

u comp Ω i = H i (X; u exp ∂Ω i ) i = 1, ..., m ≡ 19 (7)
Even for a preliminary validation of the inverse procedure on the basis of pseudo-experimental (i.e., synthetic) data, the prescribed kinematic boundary conditions are exclusively those estimated by DIC, namely u exp ∂Ω i .

Sensitivity analysis

The aim of the sensitivity analysis is to determine which parameters of the chosen model can be identified with the present experimental configuration.

Derivatives of the displacement field u ≡ {u x , u y } T with respect to the model parameters to identify X represent a key ingredient of many identification procedures. Strategies to compute such quantities are commonly referred to as "sensitivity analyses" [START_REF] Kleiber | Parameter Sensitivity in Nonlinear Mechanics. Theory and Finite Element Computations[END_REF]Micharelis et al., 2006). A simple and effective method is adopted based on the direct differentiation of the governing equations (i.e., DDM strategy). The application of DDM to the present problem takes advantage from the holonomy hypothesis assumed for the nonlinear constitutive equations (for the substrates and the adhesive layer), and it is implemented in the same FE code employed for forward simulations. Let us consider the discretized finite element residual equations at the n-th loading step (sensitivities vanish where displacements are prescribed). By differentiating the above residual equation ( 8), the following relationship is satisfied

R n = R n (U n (X), X) = 0 (
∂R n ∂U n T ∂U n ∂X T = - ∂R n ∂X T (9)
Equation ( 9) provides parameter sensitivities for any steady-state mechanical problem (Micharelis et al., 2006). In the present context, residual equations over the ROI Ω particularize as

R n = Ω B T σ n (U n (X), X) dΩ -F ext = 0 (10a) ∂R n ∂U n T = Ω B T ∂σ n ∂U n T dΩ = K u ∂R n ∂X T = Ω B T ∂σ n ∂X T dΩ (10b)
where F ext , B and σ n denote the vector of external forces, the compatibility matrix and the current value of local stress tensor, respectively. An analogy can be proposed between the linearized form of finite element equations, solved iteratively to compute Newton corrections, and the sensitivity analysis problem (9). Both problems are linear with respect to the unknowns and are governed by the assembled tangent matrix K u . For this reason, the RHS of Equations ( 9) and ( 10) is referred to as "pseudo-load" vector. From a numerical standpoint, parameter sensitivities are easily computed at the end of each loading step of the finite element analysis, when convergence by Newton iterations is achieved, by solving only the linear system (9), irrespectively of the parameter vector dimension n X . Since the tangent stiffness matrix is already available at each Newton iteration, the sensitivity analysis requires only the pseudo-load vector as additional datum. Such a vector is computed by assembling contributions from each element, analogously to what is performed for vector F int gathering the internal forces. The above strategy leads to faster identification procedures, thereby avoiding numerous forward analyses required by finite differences schemes.

Figure 9 shows the sensitivity of horizontal displacement fields u x with respect to the parameter δ t , computed by the DDM at different measurement instants i during the test. In spite of the noisy environment, at the beginning of the test the above parameter sensitivity is significant over almost the whole sub-domain, even if it peaks close to the adhesive layer (i = 4).

As the crack propagates, horizontal measurements are less sensitive to such parameter, and non-vanishing sensitivities tend to concentrate around the adhesive joint. The peak values move toward the sound part of the adhesive layer, and indicate approximately the current location of the process zone within the adhesive joint (i = 11, 15, 19). Last, it is worth noting that displacement sensitivities with respect to parameter variations δ n /δ n and δ t /δ t , where overlined quantities denote reference values, are significantly higher (especially the former) than those concerning the normal work-of-separation φ n . This result suggests not to add a further unknown by allowing ratio r to vary.

Reconstruction of interface tractions via parameter identification

Since measured boundary conditions are used, and the pictures at hand (i.e., with a poor texture) will induce high measurement uncertainties, the identification technique needs to be robust enough to allow for low signal / noise ratios. The optimal parameters X are obtained as a minimization of a suitable discrepancy norm [START_REF] Tarantola | Inverse Problems Theory. Methods for Data Fitting and Model Parameter Estimation[END_REF][START_REF] Maier | Some innovative industrial prospects centered on inverse analyses[END_REF]. Such a function quantifies the discrepancy between measured quantities (i.e., displacement of inner nodes of the ROI provided by DIC), and those computed as functions of the unknown model parameters. Let m be the number of available measurement instants, n n the number of nodal displacements selected. The adopted objective function is the least-squares norm

X = arg min X ω(X) = m i=1 ω i (X) (11) ω i (X) = 1 (W i x ) 2 nn j=1 [u exp x -u comp x (X)] 2 j + 1 (W i y ) 2 nn j=1 u exp y -u comp y (X) 2 j
where weights W i x and W i y at each measurement instant i are defined as the maximum displacements (in absolute value, and after the subtraction of the average rigid body motions) measured at that instant, along the horizontal and vertical direction, respectively. The contributions of all m discrepancies ω i are additively included in a unique scalar objective function to minimize. The discrepancy between the experimental and computed displacements is based of n n = 80 nodal values and m = 19 measurement instants, concerning both aluminum alloy strips. The number of used experimental data points, m × n r = 19 × 80 = 1520, is to be compared to the number of unknowns, 3. This large ratio is responsible for the robustness of the procedure.

Artificial data

The displacements inside the discretized domain Ω are generated numerically by using an a priori assumed parameter set X. The parameter estimates X resulting from the identification procedure are thus compared with the original set, and the error is assessed. The reference parameters are φ n = 1.0 N mm -1 , δ n = 5 µm and δ t = 50 µm. This parameter set specifies an adhesive layer capable of developing maximum normal and tangential tractions equal to T n max = 74 MPa and T t max = 103 MPa, under pure mode I and II, respectively. The mechanical parameters for the aluminum alloy sheet are as follows: Young's modulus E = 72 GPa, Poisson's ratio ν = 0.33, yield stress σ 0 = 180 MPa, yield offset α = 0.06 and hardening exponent n = 7.

The minimization in the parameter space of the objective function ω(X) is performed through gradient-based (first-order) algorithms. The most recent Trust Region Method [START_REF] Coleman | An interior trust region approach for nonlinear minimization subject to bounds[END_REF]Matlab, 2003) allows for box-constraints by an interior-point strategy. This provision is particularly important when an unknown model parameter entering constitutive models and input of finite element codes belong to a feasible interval, and the initialization set is not necessarily close to the minimum point.

To validate the above procedure, first artificial data are considered with no noise. In this situation, the three parameters are estimated correctly in a few iterations, up to machine precision. No dependence of the parameter estimates X on the algorithm initialization, which might have suggested the presence of secondary minima, was observed. Thereafter, to investigate the robustness of the inverse analysis procedure, displacements inside the subdomain Ω are corrupted by means of an additive, uncorrelated (over time and space) Gaussian noise e

y exp i = H i (X; u exp ∂Ω i ) + e i (12)
with zero average and standard deviation s. For the sake of simplicity, prescribed boundary displacements provided by DIC, namely u exp ∂Ω i , which are already affected by measurement uncertainties, are left unchanged. If the parameter estimates X resulting from the inverse procedure are far from the correct values (a conventional error threshold of 15 % is given), the inverse problem is expected to be ill-posed, and non-identifiability should be invoked with regard to one or more parameters. Values of standard devia- s (pixels) 0.1 0.6 1.0 2.0 φ n (N mm -1 ) 1.0 (2 %) 0.9 (10 %) 1.0 (3 %) 0.9 (6 %)

δ n (µm) 5.0 (1 %) 5.1 (1 %) 5.6 (12 %) 7.4 (48 %)

δ t (µm) 50 (1 %) 48 (5%) 52 (3 %) 49 (2 %) ω( X) 0.4 7.5 14.3 30.6
tions considered for white noise are as follows: s = 0.1, 0.6, 1.0, 2.0 pixels.

Table 1 shows the parameter estimates resulting from data with increasing level of noise, together with their percentage error (in parentheses) and values of objective function at the minimum point X. Values of the objective function at the minimum point, ω( X), increase significantly with the noise level (Figure 10). In the presence of noise levels with a standard deviation s ≤ 1.0 pixel, the resulting parameter estimates X provided by the above inverse procedure are rather acceptable (i.e., the relevant percentage error does not exceed the chosen threshold of 15 %). For s ≥ 2 pixels instead, the parameter δ n is no longer identifiable by the above numerical procedure. It is worth emphasizing that a standard deviation s = 2 pixels implies that the single Gaussian perturbation term belongs to the interval e ∈ [-6, +6] pixels, with a probability of 99 %. With reference to the DIC procedure, such a noise intensity on kinematic data is very high [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF], and therefore the inverse procedure can be regarded as sufficiently robust in the present case.

Experimental data

Experimental data of Section 2 are considered. Through the bilinear shape functions defined over each finite element, kinematic data provided by the DIC algorithm are interpolated at the discretization nodes belonging to the mechanical model. Parameter estimates obtained by the inverse procedure are φ n = 2.5 N mm -1 , δ n = 8.5 µm, δ t = 85 µm, and hence φ t = rφ n = 9.8 N mm -1 . The value of the objective function at the minimum point is equal to ω( X) = 7.5. It is worth noting that the latter has a level close to the noisy case resting on independent Gaussian perturbations, with standard deviation s = 0.6 pixel (Table 1) for which the identification was deemed acceptable.

Although the cohesive model adopted herein describes a nonlinear behavior of the adhesive layer, useful indications are obtained concerning the elastic properties of the adhesive layer. From the estimated parameters, the (de-coupled) tangent stiffness at the origin (∆ n , ∆ t → 0) reads

• along the normal direction, ∂T n /∂∆ n = φ n /δ 2 n = 34.1 N mm -3 ,
• and along the tangential direction, ∂T t /∂∆ t = 2rφ n /δ mm -3 [START_REF] Högberg | Mixed mode cohesive law[END_REF], along the normal and tangential direction, respectively. These results show that the loading conditions generated at the microscale during the test proposed herein are rather severe in order to calibrate a mixedmode cohesive model. Mode mixity varies continuously during the debonding experiment along the adhesive layer. Interface tractions computed by the identified model (Figure 11) are reasonable from an engineering standpoint, when compared to experimental data on structural adhesives available in the literature [START_REF] Högberg | Constitutive behaviour of mixed mode loaded adhesive layer[END_REF][START_REF] Leffler | Shear behaviour of adhesive layers[END_REF][START_REF] Hua | Continuum damage modelling of environmental degradation in joints bonded with EA9321 epoxy adhesive[END_REF] when taking into account the relevant layer thickness. Under pure mode I and mode II, the calibrated model predicts maximum tractions T n max = 107 MPa and T t max = 99 MPa. To assess the role played by the dimensionless T n max , predicted by the model with the above estimates (Table 2), are also varying with r. Further experimental tests under different loading conditions are needed to tune properly the parameter r.

10 i = 7 i = 13 i = 0 n ∆ < 0 n ∆ ≥ crack front joint length [mm] t T [MPa] 10 i = 7 i = 13 i = crack front joint length [mm]

Closing remarks and outlook

In this study the mechanical properties of an adhesive layer in a GLARE assembly, manufactured with an industrial process, have been estimated The analysis performed herein, and resting only on measured displacements, is quite general in terms of applicability. For instance, delamination in composite materials, friction and sliding in metal/metal or hybrid assemblies may be analyzed with the same type of procedure. Last, the measured displacement field may be enriched as in extended finite elements [START_REF] Black | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] to better capture displacement discontinuities [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF][START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF].

Figure 1 :

 1 Figure 1: Joined sample (surface area: 297 × 50 mm 2 ) cut from the original laminate.The epoxy adhesive layer between the plate and the Z-shaped stringer is the focus of the present study.

  perienced within the ROI by DIC.

Figure 2 :

 2 Figure 2: Detail of the adhesive joint (in red), showing a tapered profile to allow for the insertion of the Z-shaped reinforcement on the skin laminate. Both laminates are made of three aluminum alloy (light gray) and two epoxy pre-preg (dark gray) layers. The small sub-domain monitored by DIC is located at the right end of the picture, where the joint thickness stabilizes to 60 µm (constant along the remaining length). The magnification used herein is less than that adopted in the pictures analyzed with DIC.

Figure 3 :Figure 4 :

 34 Figure 3: Non-conventional test on a joined sample, up to complete debonding.

Figure 5 :

 5 Figure 5: Measurements at the microscale provided by DIC through a finite element discretization. Deformed meshes of joined substrates are shown at given measurement instants i during the debond test.

  transforming the global coordinate system to the local one (in the present case R = I). Shape functions, gathered in matrix N , vary linearly along the element length, and are constant through the layer thickness. The traction T at the joint interface is expressed as a function of the relative displacement vector ∆ [Equations (2a)-(2b)]. At each interface element e in the adopted discretization, the vector of internal nodal forces, F

Figure 7 :

 7 Figure 7: Profiles of normal tractions along the modeled interface, at different measurement instants i, computed with coarse (solid line) and fine (open circles) meshes.

Figure 7 Figure 8 :

 78 Figure 7 shows the normal traction profiles along the modeled joint interface computed with the coarse and fine discretizations. The two profiles are very close. This result shows that the constraints imposed by DIC do not degrade significantly the results in terms of traction profiles since the measurements are performed close enough to the cohesive zone. Twenty digital images are considered to analyze the debonding process within the ROI. Full field measurements are available at m = 19 measurement instants, since the first image (i = 0) is the reference. Each forward analysis requires approximately 15 minutes on a standard PC (2.2 GHz Quad Core, and 3 Gigabyte RAM). The finite element simulation of the debonding process within the monitored sub-domain is shown in Figure 8. For those results, a fixed ratio of r = 6 was chosen. Displacements provided by DIC along the boundary of the monitored zone, u exp ∂Ω i , are prescribed at different instants i (i = 1, ..., 19) during the test. Figure 8 shows deformed meshes for four different steps of the debonding process.The nonlinear relationship between the model parameters for the adhesive

Figure 9 :

 9 Figure 9: Sensitivity of the horizontal displacement field u x with respect to parameter δ t , ∂u x /∂δ t , at different instants i

Figure 10 :

 10 Figure 10: Validation of the identification procedure on the basis of artificial data. Minimization process of the objective function by Trust Region algorithm, in the presence of a Gaussian (non-correlated) noise on measured displacements, with increasing standard deviation s, namely s = 0, 0.1, 0.6, 1 and 2 pixels.

Figure 11 :

 11 Figure 11: Change of normal and tangential traction profiles at the adhesive interface, T n and T t , resp., re-constructed by the proposed identification procedure at different measurement instants i, with the a priori assumed ratio r = 4.

Figure 11

 11 Figure 11 shows the profiles of tangential and normal tractions along the adhesive layer computed by the identified model at different measurement instants i (8 ≤ i ≤ 19) during the test. The profiles of normal tractions T n exhibit a classical bell shape with a softening branch for tensile values when ∆ n ≥ 0, while compressive values are observed in the remaining part of the adhesive layer (T n < 0 when ∆ n < 0). During the debonding process, profiles of normal tractions T n translate progressively towards the left part of the monitored sub-domain as the crack propagates. Simultaneously, their peak values change, since the mode mixity varies continuously during the experiment along the adhesive layer. The traction-free part of the interface, corresponding to a macroscopic crack observable in digital images (see Figure 4), becomes larger, starting from the right end of the ROI, and the compressed part (T n < 0) at the left side progressively reduces until it completely vanishes. As a consequence of measurement uncertainties, small values of tensile tractions (≤ 10 MPa) are observed on the left end of the ROI. The profiles of tangential tractions T t appear to be smoother with respect to the normal ones, since the horizontal displacements detected by DIC and prescribed along the model boundaries are larger and yield better signal to noise ratios. Moreover, horizontal boundary displacements prescribed at each aluminum alloy strip exhibit a rather limited variation along the layer. Therefore, bell-shape profiles for tangential tractions are less obvious. The pre-peak response is clearly distinguished only at the first measurement instants.

Figure 12 :

 12 Figure 12: Interface tractions T n at measurement instants i = 8, 11, 14, 17, identified by the proposed inverse procedure for different values of the parameter 3 ≤ r ≤ 6.

  through an inverse method based exclusively on full field displacements provided by DIC. Non-conventional tests were designed to generate mixed-mode fracture conditions within the adhesive layer. Displacements at the microscale estimated by DIC are adopted to specify Dirichlet boundary conditions for a small sub-domain, including part of the adhesive layer subjected to debonding and part of the aluminum alloy substrates.Local kinematic measurements, provided by DIC in a format that is fully consistent with conventional finite element discretizations, enabled us to focus the inverse analysis on the sole ROI of the DIC analysis. This procedure allows one to model only a small part of the experiment since the effect of the remote boundary conditions is accounted for by the measured displacements on the boundary of the ROI. The inner displacements are then used to tune the parameters of the cohesive law. Kinematic data are thus adopted simultaneously to drive the mechanical simulation of the test (Dirichlet conditions prescribed along the external boundary), and as comparison term (inside the domain) included in the objective function. The parameters were identified and a good match was found between measured and computed (inner) displacements at the optimum point.

5

  Measurements at the microscale provided by DIC through a finite element discretization. Deformed meshes of joined substrates are shown at given measurement instants i during the debond test. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 6 Interface element with a finite thickness. . . . . . . . . . . . . 13 7 Profiles of normal tractions along the modeled interface, at different measurement instants i, computed with coarse (solid line) and fine (open circles) meshes. . . . . . . . . . . . . . . . 16 8 Mechanical model (in plane stress) to simulate debonding within the monitored domain, in which boundary displacements estimated by DIC are prescribed. Deformed mesh at given measurement instants i during the tests. . . . . . . . . . . . . . . 18 9 Sensitivity of the horizontal displacement field u x with respect to parameter δ t , ∂u x /∂δ t , at different instants i during the test simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 10 Validation of the identification procedure on the basis of artificial data. Minimization process of the objective function by Trust Region algorithm, in the presence of a Gaussian (noncorrelated) noise on measured displacements, with increasing standard deviation s, namely s = 0, 0.1, 0.6, 1 and 2 pixels. . . 24 11 Change of normal and tangential traction profiles at the adhesive interface, T n and T t , resp., re-constructed by the proposed identification procedure at different measurement instants i, with the a priori assumed ratio r = 4. . . . . . . . . . . . . . . 27

Table 1 :

 1 Identified parameters as functions of the standard uncertainty level s of the displacement data.

Table 2 :

 2 Values of the different identified parameters as a function of the r ratio.

	r	3.0 4.0	5.0	6.0
	φ n (N mm -1 )	3.0 2.5	2.2	1.9
	δ n (µm)	9.4 8.5	7.9	7.5
	δ t (µm)	80	85	90	94
	φ t = rφ n (Nmm -1 ) 9.1 9.8 10.6 11.4
	T n max (MPa)	118 107	99	94
	T t max (MPa)	97	99	101 103
	∂ Tn ∂ ∆n | 0 (GPa mm -1 ) 34	34	34	34
	∂ Tt ∂ ∆t | 0 (GPa mm -1 ) 2.8 2.7	2.6	2.6
	ω( X)	7.5 7.5	7.5	7.5
	ratio r = φ t /φ n , the same identification process was repeated by assuming
	different values for r. Table 2 shows for each value of r various parameter

Interface tractions T n at measurement instants i = 8, 11, 14, 17, identified by the proposed inverse procedure for different values of the parameter 3 ≤ r ≤ 6. . . . . . . . . . . . . . . . . . 29

Acknowledgements

The Authors thank ALENIA for providing composite laminates investigated in a research project between LMT-Cachan and Politecnico di Milano, supported by the EU project "Knowledge-based Multicomponent Materials for Durable and Safe Performance" (KMM-NoE, NMP3-CT-2004-502243), and by the "young researcher fund 2006" of Politecnico di Milano.


Gilles Audemard

Steve Bellart

Louènas Bounia
email: bounia@cril.fr

Frédéric Koriche
email: koriche@cril.fr

Jean-Marie Lagniez
email: lagniez@cril.fr

Pierre Marquis
email: marquis@cril.fr

Trading Complexity for Sparsity in Random Forest Explanations

Random forests have long been considered as powerful model ensembles in machine learning. By training multiple decision trees, whose diversity is fostered through data and feature subsampling, the resulting random forest can lead to more stable and reliable predictions than a single decision tree. This however comes at the cost of decreased interpretability: while decision trees are often easily interpretable, the predictions made by random forests are much more difficult to understand, as they involve a majority vote over multiple decision trees. In this paper, we examine different types of reasons that explain "why" an input instance is classified as positive or negative by a Boolean random forest. Notably, as an alternative to prime-implicant explanations taking the form of subset-minimal implicants of the random forest, we introduce majoritary reasons which are subset-minimal implicants of a strict majority of decision trees. For these abductive explanations, the tractability of the generation problem (finding one reason) and the optimization problem (finding one minimumsized reason) are investigated. Unlike prime-implicant explanations, majoritary reasons may contain redundant features. However, in practice, prime-implicant explanationsfor which the identification problem is DP-complete -are slightly larger than majoritary reasons that can be generated using a simple linear-time greedy algorithm. They are also significantly larger than minimum-sized majoritary reasons which can be approached using an anytime PARTIAL MAXSAT algorithm.

Introduction

Over the past two decades, rapid progress in statistical machine learning has led to the deployment of models endowed with remarkable predictive capabilities. Yet, as the spectrum of applications using statistical learning models becomes increasingly large, explanations for why a model is making certain predictions are ever more critical. For example, in medical diagnosis, if some model predicts that an image is malignant, then the physician may need to know which features in the image have led to this classification. Similarly, in the banking sector, if some model predicts that a customer commits fraud, then the banker might want to know why. Therefore, deriving explanations for why certain predictions are made is essential for securing user confidence in machine learning technologies [START_REF] Miller | Explanation in artificial intelligence: Insights from the social sciences[END_REF][START_REF] Molnar | Interpretable Machine Learning -A Guide for Making Black Box Models Explainable[END_REF].

This paper focuses on classifications made by random forests, a popular ensemble learning method that constructs multiple randomized decision trees during the training phase, and predicts by taking a majority vote over the base classifiers [START_REF] Breiman | Random Forests[END_REF]. Since decision tree randomization is achieved by essentially coupling data subsampling (or bagging) and feature subsampling, random forests are easy to implement, with few tuning parameters. Furthermore, they often make accurate and robust predictions in practice, even for small data samples and high-dimensional feature spaces [START_REF] Biau | Analysis of a Random Forests Model[END_REF]. For these reasons, random forests have been used in various applications including, among others, computer vision [START_REF] Criminisi | Decision Forests for Computer Vision and Medical Image Analysis[END_REF], crime prediction [START_REF] Bogomolov | Once Upon a Crime: Towards Crime Prediction from Demographics and Mobile Data[END_REF], ecology [START_REF] Cutler | Random Forests for Classification in Ecology[END_REF], genomics [START_REF] Chen | Random forests for genomic data analysis[END_REF], and medical diagnosis [START_REF] Azar | A random forest classifier for lymph diseases[END_REF].

However, random forests are often considered far less interpretable than decision trees. Thus, while many XAI queries are tractable for decision trees, they are not tractable for random forests [START_REF] Audemard | On the Computational Intelligibility of Boolean Classifiers[END_REF]. Especially, the prediction made on a given data instance can be easily interpreted by reading the direct reason furnished by the classifier. For a decision tree, it is the unique root-to-leaf path that covers the instance (also known as its path-restricted explanation [START_REF] Izza | On Explaining Decision Trees[END_REF]. By contrast, no such direct reason has been defined so far for a random forest. More generally, a key issue in random forests is to infer abductive explanations, that is, to capture in concise terms why a data instance is classified as positive or negative by the model ensemble.

Related Work. Explaining random forest predictions has received increasing attention in recent years [START_REF] Bénard | Interpretable Random Forests via Rule Extraction[END_REF][START_REF] Choi | On Symbolically Encoding the Behavior of Random Forests[END_REF][START_REF] Izza | On Explaining Random Forests with SAT[END_REF]. In the classification setting, [START_REF] Choi | On Symbolically Encoding the Behavior of Random Forests[END_REF] and [START_REF] Izza | On Explaining Random Forests with SAT[END_REF] have focused on prime-implicant explanations, also known as sufficient reasons [START_REF] Darwiche | On the Reasons Behind Decisions[END_REF]. Informally, if we view any random forest classifier as a Boolean function f , then a prime-implicant explanation for classifying a data instance x as positive by f is a subsetminimal implicant t of f covering x. By construction, t is an abductive explanation involving only relevant features, since removing any feature from t would question the fact that t explains the way x is classified by f . Note that if f is described by a single decision tree, then generating a primeimplicant explanation for any data instance x can be done in linear time. Yet, in the general case where f is represented by an arbitrary number of decision trees, the problem of determining whether a given term is a prime-implicant explanation for an instance given a random forest has recently been shown DP-complete [START_REF] Izza | On Explaining Random Forests with SAT[END_REF]. Despite this intractability statement, the empirical results reported by the authors show that a MUS-based algorithm for computing such explanations can prove efficient in practice.

In addition to model-based explanations described above, model-agnostic explanations can be applied to random forests. Notably, the LIME method [START_REF] Ribeiro | Explaining the Predictions of Any Classifier[END_REF] extrapolates a linear threshold function g from the behavior of the random forest f around an input instance x. For the ANCHOR method [START_REF] Ribeiro | Anchors: High-Precision Model-Agnostic Explanations[END_REF], the extrapolated model g takes the form of a decision rule. Yet, even if in both cases a prime implicant of g can be easily computed, the resulting explanation is not guaranteed to be an abductive explanation for x given f since g is an approximation of f . Contributions. In this paper, we introduce new notions of abductive explanations: direct reasons, which extend to the case of random forests the corresponding notion defined for decision trees, and majoritary reasons, which are abductive explanations taking into account the averaging rule of random forests. Informally, a majoritary reason for classifying an instance x as positive by some random forest F is a term t that covers x and is a subset-minimal implicant of a strict majority of decision trees in F .

What makes direct and majoritary reasons valuable is the possibility of inferring them in a tractable way, whilst there is no similar tractability result when dealing with sufficient reasons, unless P = NP. In the following, we examine the tractability of both the generation problem (finding one explanation) and the optimization problem (finding one minimum-sized explanation) for direct and majoritary reasons. Direct reasons (that coincide with minimum-sized direct reasons) and majoritary reasons can be derived in polynomial time. By contrast, identifying minimum-sized majoritary reasons is NP-complete, and identifying minimumsized prime-implicant explanations is Σ p 2 -complete. Based on these results, we provide algorithms for deriving random forest explanations, enabling an empirical comparison. Our experiments made on standard benchmarks assess both the runtime complexity of finding abductive explanations and the sparsity of such explanations (i.e., how much parsimonious they are). In a nutshell, majoritary reasons and minimum-sized majoritary reasons offer interesting compromises in comparison to, respectively, prime-implicant explanations and minimum-sized prime-implicant explanations. Indeed, even if in theory for every majoritary reason, there exists a prime-implicant explanation that is not larger, in practice the size of majoritary reasons are generally smaller than those of prime-implicant explanations. Furthermore, the computational effort to be spent for deriving ma-joritary reasons is smaller than the one required by primeimplicant explanations. Similarly, in practice, minimumsized majoritary reasons outperform minimum-sized primeimplicant explanations, in the sense that deriving minimumsized prime-implicant explanations is often too computationally demanding. Using an anytime PARTIAL MAXSAT solver for minimizing the size of majoritary reasons, we show how to derive abductive explanations which are typically shorter than all other forms of abductive explanations considered in the paper. A full-proof version of the paper is available at www.cril.univ-artois.fr/expekctation/papers. html.

Preliminaries

For an integer n, let [n] = {1, • • • , n}. By F n we denote the class of all Boolean functions from {0, 1} n to {0, 1}, and we use

X n = {x 1 , • • • , x n } to denote the set of input Boolean variables. Any Boolean vector x ∈ {0, 1} n is called an in- stance. For any function f ∈ F n , an instance x ∈ {0, 1} n is called a positive instance of f if f (x) = 1, and a negative instance if f (x) = 0.
We refer to f as a propositional formula when it is described using the Boolean connectives ∧ (conjunction), ∨ (disjunction) and ¬ (negation), together with the constants 1 (true) and 0 (false). f is satisfiable if it has a positive instance, and it is unsatisfiable otherwise. A literal l i is a variable x i or its negation ¬x i , also denoted x i . A term t is a conjunction of literals, and a clause c is a disjunction of literals. A DNF formula is a disjunction of terms and a CNF formula is a conjunction of clauses. The set of variables occurring in a formula f is denoted Var (f ). In the following, we shall often treat instances as terms, and terms as sets of literals. For an assignment z ∈ {0, 1} n , the corresponding term is

t z = n i=1
x zi i where x 0 i = x i and x 1 i = x i A term t covers an assignment z if t ⊆ t z . An implicant of a Boolean function f is a term that implies f , that is, a term t such that f (z) = 1 for every assignment z covered by t.

A prime implicant of f is an implicant t of f such that no proper subset of t is an implicant of f . A (Boolean) decision tree on X n is a binary tree T , each of whose internal nodes is labeled with one of n input variables, and whose leaves are labeled 0 or 1. Without loss of generality, every variable is supposed to occur at most once on any root-to-leaf path. The value T (x) of T on an input instance x is given by the label of the leaf reached from the root as follows: at each node go to the left (resp. right) child if the input value of the corresponding variable is 0 (resp. 1). A (Boolean) random forest on X n is an ensemble

F = {T 1 , • • • , T m }, where each T i (i ∈ [m]
) is a decision tree on X n , and such that the value F (x) is given by

F (x) = 1 if 1 m m i=1 T i (x) > 1 2 0 otherwise. The size of F is given by |F | = m i=1 |T i |,
where |T i | is the number of nodes occurring in T i . The class of decision trees

x 4 T 1 0 x 2 1 x 3 0 x 1 0 1 x 2 T 2 x 1 1 0 x 4 0 1 x 3 T 3 x 2 x 2 x 1 x 4 0 x 1 0 1 0 x 4 0 1 0 1
Figure 1: A random forest for recognizing Cattleya orchids.

The left (resp. right) child of any decision node labelled by x i corresponds to the assignment of x i to 0 (resp. 1).

on X n is denoted DT n , and the class of random forests with at most m decision trees (with m ≥ 1) 1 is composed of three decision trees. It separates Cattleya orchids from other orchids using the following features: x 1 : "has fragrant flowers", x 2 : "has one or two leaves", x 3 : "has large flowers", and x 4 : "is sympodial".

over DT n is denoted RF n,m . Finally, RF n = m≥1 RF n,m and RF = n≥1 RF n . Example 1. The random forest F = {T 1 , T 2 , T 3 } in Figure
We conclude this section with two important properties of random forests, that are useful to prove forthcoming results. The first property is related to the fact that any decision tree T can be transformed into its negation ¬T ∈ DT n , by simply inverting the labels of its leaves. Negating a random forest can also be achieved in polynomial time: Proposition 1. There exists a linear-time algorithm that computes a random forest ¬F ∈ RF n,m equivalent to the negation of a given random forest F ∈ RF n,m .

Furthermore, it is well-known that any decision tree T can be encoded in linear time into an equivalent disjunction of terms DNF(T ), where each term coincides with a 1-path (i.e., a path from the root to a leaf labeled with 1), but also as a conjunction of clauses CNF(T ), where each clause is the negation of a term describing a 0-path. Yet, when switching to random forests, the picture is quite different: Proposition 2. Any CNF or DNF formula can be converted in linear time into an equivalent random forest, but there is no polynomial-space translation from RF to CNF or to DNF.

Random Forest Explanations

The key focus of this study is to explain why a random forest classifies some data instance as positive or negative. This calls for a notion of abductive explanation.1 Specifically, an abductive explanation for an instance x ∈ {0, 1} n given a Boolean function f ∈ F n is a term t that covers x and is an implicant of f (resp. ¬f ) if f (x) = 1 (resp. f (x) = 0). Such an abductive explanation always exists, since t = t x is such a (trivial) explanation. So, in the rest of this section, we shall mainly concentrate on sparse forms of abductive explanations.

Direct Reasons. For a decision tree T ∈ DT n and a data instance x ∈ {0, 1} n , the direct reason for x given T is the term t T

x corresponding to the unique root-to-leaf path of T that covers x. This simple form of abductive explanation can be extended to random forests as follows: Definition 1. Let F = {T 1 , . . . , T m } be a random forest in RF n,m , and x ∈ {0, 1} n be an instance. The direct reason for x given F is the term t F

x defined by

t F x = Ti∈F :Ti(x)=1 t Ti x if F (x) = 1 Ti∈F :Ti(x)=0 t Ti x if F (x) = 0
By construction, t F

x is an abductive explanation for x given F that can be computed in

O(n|F |) time. Example 2. Based on Example 1, consider the instance x = (1, 1, 1, 1). Since F (x) = 1, x is recognized as a Cattleya orchid. The direct reason for x given F is t F x = x 1 ∧ x 2 ∧ x 3 ∧ x 4 . Consider now the instance x ′ = (0, 1, 0, 0), which is not recognized as a Cattleya orchid, since F (x) = 0. The direct reason for x ′ given F is t F x ′ = x 2 ∧ x 3 ∧ x 4 . Prime-Implicant Explanations.
Another valuable notion of abductive explanation is the one of prime-implicant explanation [START_REF] Shih | A Symbolic Approach to Explaining Bayesian Network Classifiers[END_REF]) also referred to as sufficient reason (Darwiche and Hirth 2020). In the setting of random forests, such explanations can be defined as follows: Definition 2. Let F ∈ RF n be a random forest and x ∈ {0, 1} n be an instance. A prime-implicant explanation for x given F is a prime implicant t of F (resp. ¬F ) if F (x) = 1 (resp. F (x) = 0) such that t covers x. Example 3. For our running example, x 2 ∧ x 3 ∧ x 4 and x 1 ∧ x 4 are the prime-implicant explanations for x given F . x 4 and x 1 ∧ x 3 are the prime-implicant explanations for x ′ given F .

Importantly, all features occurring in a prime-implicant explanation t are relevant. Indeed, removing any literal from t would question the fact that t implies F . Note that the direct reason t F

x for x given F may contain arbitrarily many more features than a prime-implicant explanation for x given F , since this is already known in the case where F consists of a single decision tree [START_REF] Izza | On Explaining Decision Trees[END_REF].

The problem of identifying a prime-implicant explanation t for an input instance x ∈ {0, 1} n given a random forest F ∈ RF n , has recently been shown DP-complete (Izza and Marques-Silva 2021). In fact, even the apparently simple task of checking whether t is an implicant of F is already hard: Proposition 3. Let F be a random forest in RF n and t be a term over X n . Deciding whether t is an implicant of F is coNP-complete.

The above result is in stark contrast with the computational complexity of checking whether a term t is an implicant of a decision tree T . This task can be solved in polynomial time, using the fact that T can be converted (in linear time) into a clausal form CNF(T ), together with the fact that testing whether t implies CNF(T ) can be done in O(n|T |) time. That mentioned, in the case of random forests, the implicant test can be achieved via a call to a SAT oracle: Proposition 4. Let F = {T 1 , . . . , T m } be a random forest of RF n,m , and t be a (satisfiable) term over X n . Let H be the CNF formula

{(y i ∨ c) : i ∈ [m], c ∈ CNF(¬T i )} ∪ CNF m i=1 y i > m 2
where {y 1 , . . . , y m } are fresh variables and CNF

m i=1 y i > m 2 is a CNF encoding of the cardi- nality constraint m i=1 y i > m 2 . Then, t is an implicant of F if and only if H ∧ t is unsatisfiable.
Based on such an encoding, the prime-implicant explanations for an instance x given a random forest F can be characterized in terms of MUS (minimal unsatisfiable subsets), as suggested in [START_REF] Izza | On Explaining Random Forests with SAT[END_REF]. This characterization is useful because many SAT-based algorithms for computing a MUS (or even all MUSes) of a CNF formula have been pointed out for the past decade [START_REF] Audemard | Improving Glucose for Incremental SAT Solving with Assumptions: Application to MUS Extraction[END_REF][START_REF] Liffiton | Fast, flexible MUS enumeration[END_REF][START_REF] Marques-Silva | Minimal sets on propositional formulae. Problems and reductions[END_REF]. They can be leveraged for computing prime-implicant explanations.

Going one step further, a natural way to improve the intelligibility of prime-implicant explanations is to focus on parsimonious ones. 2 Specifically, a minimum-sized primeimplicant (minPI) explanation for x ∈ {0, 1} n given F ∈ RF n is a prime-implicant explanation for x given F , that is of minimal size. Example 4. For our running example, x 1 ∧ x 4 is the unique minPI explanation for x given F , and x 4 is the unique minPI explanation for x ′ given F .

As a by-product of the characterization of a primeimplicant explanation in terms of MUS (Izza and Marques-Silva 2021), a minPI explanation for x given f may be viewed as a minimum-sized MUS. Thus, we can exploit algorithms for computing minimum-sized MUSes (see e.g., [START_REF] Ignatiev | Smallest MUS Extraction with Minimal Hitting Set Dualization[END_REF]) in order to infer minPI explanations. However, identifying a minPI explanation is computationally harder than identifying a prime-implicant explanation: Proposition 5. Let F ∈ RF n , x ∈ {0, 1} n , and k ∈ N. Deciding whether there exists a minPI explanation t for x given F such that t contains at most k features is Σ p 2 -complete. Majoritary Reasons. Based on the above considerations, a natural question arises: does there exist a middle ground between direct reasons, which may include many irrelevant features but are easy to calculate, and prime-implicant explanations, which only contain relevant features but are potentially much harder to infer? Inspired by the way prime 2 Everything else being equal, shortest explanations can be viewed as better than longer explanations because they are easier to understand. Of course, sparsity is only one aspect of the intelligibility of an explanation. The quality of an explanation typically depends on the explainee (i.e., the person who asked for an explanation) (Doshi-Velez and Kim 2017).

implicants can be computed when dealing with decision trees, we can reply in the affirmative using the notion of majoritary reasons, defined as follows.

Definition 3. Let F = {T 1 , . . . , T m } be a random forest in RF n,m and x ∈ {0, 1} n be an instance. A majoritary reason for x given F is a term t covering x, such that t is an implicant of at least ⌊ m 2 ⌋ + 1 decision trees T i (resp. ¬T i ) if F (x) = 1 (resp. F (x) = 0), and for every l ∈ t, t \ {l} does not satisfy this last condition.

Example 5. Based on our running example, the majoritary reasons for x given F are x 1 ∧ x 2 ∧ x 4 , x 1 ∧ x 3 ∧ x 4 , and x 2 ∧x 3 ∧x 4 . Each of these explanations is more concise than the direct reason t F

x . For x ′ , the majoritary reasons given F are x 1 ∧x 4 , x 2 ∧x 4 , and x 1 ∧x 2 ∧x 3 . Note that the each of the majoritary reasons x 1 ∧x 2 ∧x 4 and x 1 ∧x 3 ∧x 4 for x given F includes an irrelevant literal for the task of classifying x using F since x 1 ∧ x 4 is a prime-implicant explanation for x given F . Similarly, every majoritary reason for x ′ given F contains an irrelevant literal for the task of classifying x ′ using F .

In general, the notions of majoritary reason and of primeimplicant explanation do not coincide. Indeed, a primeimplicant explanation is a prime implicant (covering x) of the forest F , while a majoritary reason is an implicant t (covering x) of a majority of decision trees in the forest F , satisfying the additional condition that t is a prime implicant of at least one of these decision trees.

A key observation justifying this difference is that even if every implicant of a Boolean function f is an implicant of the function f ∨ g, it is not always the case that every prime implicant of f is a prime implicant of f ∨ g. To this point, consider our running example and take the term t = x 1 ∧x 3 ∧x 4 . Here, t is a majoritary reason for x = (1, 1, 1, 1) given F , since t covers x, t is a prime implicant of T 1 , and t is an implicant of T 2 . Thus, t is an implicant of f = T 1 ∧ T 2 (it is a prime one), and hence an implicant of F , using the fact that F is logically equivalent to (T 1 ∧ T 2 ) ∨ (T 1 ∧ T 3 ) ∨ (T 2 ∧ T 3 ). However, t is not a prime implicant of F . Indeed, the sub-term x 1 ∧ x 4 is a prime-implicant explanation for x given F , since it is a prime implicant of F that covers x.

Viewing majoritary reasons as "weak" forms of primeimplicant explanations, they can include irrelevant features: Proposition 6. Let F = {T 1 , . . . , T m } be a random forest of RF n,m and x ∈ {0, 1} n . Unless m < 3, it can be the case that every majoritary reason for x given F contains arbitrarily many more features than any prime-implicant explanation for x given F .

What makes majoritary reasons valuable is that they are abductive and can be generated in linear time. The evidence that any majoritary reason t for x given F is an abductive explanation comes directly from the fact that if t implies a majority of decision trees in F , then it is an implicant of F (note that the converse implication does not hold in general).

The tractability of generating majoritary reasons lies in the fact that they can be found using a simple greedy algorithm. For the case when F (x) = 1, start with t = t x , and iterate over the literals l of t by checking whether t deprived of l is an implicant of at least ⌊ m 2 ⌋ + 1 decision trees of F . If so, remove l from t and proceed to the next literal. Once all literals in t x have been examined, the final term t is by construction an implicant of a majority of decision trees in F , such that removing any literal from it would lead to a term that is no longer an implicant of this majority. So, t is by construction a majoritary reason. The case when F (x) = 0 is similar, by simply replacing each T i with its negation in F . This greedy algorithm runs in O(n|F |) time, using the fact that, on each iteration, checking whether t is an implicant of T i (for each i ∈ [m]) can be done in O(n|T i |) time.

By analogy with minimum-sized prime-implicant explanations, a natural way of improving the quality of majoritary reasons is to seek for the most parsimonious ones. Formally, a minimum-sized majoritary (minMAJ) reason for an instance x ∈ {0, 1} n given a random forest F ∈ RF n is a majoritary reason for x given F , that is of minimal size.

Example 6. For our running example, the three majoritary reasons for x given F are minMAJ reasons. Contrastingly, among the majoritary reasons for x ′ given F , only x 1 ∧ x 4 and x 2 ∧ x 4 are minMAJ reasons.

Unsurprisingly, the optimization task for majoritary reasons is more demanding than the generation task. Still, min-MAJ reasons are easier to identify than minPI reasons:

Proposition 7. Let F ∈ RF n , x ∈ {0, 1} n , and k ∈ N.
Deciding whether there exists a minMAJ reason t for x given F such that t contains at most k features is NP-complete.

A common approach for handling NP-optimization problems is to rely on modern constraint solvers. From this perspective, recall that a PARTIAL MAXSAT problem consists of a pair (C soft , C hard ) where C soft and C hard are (finite) sets of clauses. The goal is to find a Boolean assignment that maximizes the number of clauses c in C soft that are satisfied, while satisfying all clauses in C hard .

Proposition 8. Let F ∈ RF n,m and x ∈ {0, 1} n . Let (C soft , C hard ) be an instance of the PARTIAL MAXSAT problem such that:

C soft = {x i : x i ∈ t x } ∪ {x i : x i ∈ t x } C hard = {(y i ∨ c |x ) : i ∈ [m], c ∈ CNF(T ± i )} ∪ CNF m i=1 y i > m 2
where c |x = c ∩ t x is the restriction of c to the literals in t x , {y 1 , . . . , y m } are fresh variables,

T ± i = T i (i ∈ [m]) if F (x) = 1, T ± i = ¬T i if F (x) = 0, and CNF( m i=1 y i > m 2 ) is a CNF encoding of the constraint m i=1 y i > m 2
. Let z * be an optimal solution of (C soft , C hard ). Then, t x ∩ t z * is a minMAJ reason for x given F .

Thanks to this characterization result, one can leverage the numerous algorithms that have been developed so far for PARTIAL MAXSAT (see e.g. [START_REF] Ansótegui | SAT-based MaxSAT algorithms[END_REF][START_REF] Morgado | MSCG: Robust Core-Guided MaxSAT Solving[END_REF][START_REF] Narodytska | Maximum Satisfiability Using Core-Guided MaxSAT Resolution[END_REF][START_REF] Saikko | LMHS: A SAT-IP Hybrid MaxSAT Solver[END_REF]) in order to compute minMAJ reasons.

Experiments

Experimental Setup. The empirical protocol we considered was as follows. We have focused on 15 datasets for binary classification, which are standard benchmarks from the repositories Kaggle (www.kaggle.com), OpenML (www.openml.org), or UCI (archive.ics.uci.edu/ml/). These datasets are compas, placement, recidivism, adult, ad data, mnist38, mnist49, gisette, dexter, dorothea, farm-ads, higgs boson, christine, gina, and bank. mnist38 and mnist49 are subsets of the mnist dataset, restricted to the instances of 3 and 8 (resp. 4 and 9) digits. Additional information about the datasets, (especially the number and types of features, the number of instances), and about the random forests that have been trained (especially, the number of Boolean features used, the number of trees, the depth of the trees, the mean accuracy) are available at www.cril.univ-artois.fr/ expekctation/.

Categorical features have been treated as arbitrary numbers (the scale is nominal). As to numeric features, no data preprocessing has taken place: these features have been binarized on-the-fly by the random forest learning algorithm. For this learner, we have used the version 0.23.2 of the Scikit-Learn library [START_REF] Pedregosa | Scikitlearn: Machine Learning in Python[END_REF]. The maximal depth of any decision tree in a forest has been bounded at 8. All other hyper-parameters of the learning algorithm have been set to their default value, except the number of trees. We made some preliminary tests for tuning this parameter in order to ensure that the accuracy was good enough.

For every benchmark b, a 10-fold cross validation process has been achieved: a set of 10 random forests have been computed and evaluated from the labelled instances of b, partitioned into 10 parts. One part was used as the test set and the remaining 9 parts as the training set for generating a forest. The classification performance on b was measured using the mean accuracy obtained over the 10 random forests. For each benchmark b, each random forest F , and a pool of 25 instances x drawn at random from the test set (leading to 250 instances per dataset), we have run the algorithms described in the previous section for deriving the direct reason for x given F , a prime-implicant explanation for x given F , a majoritary reason for x given F , a minPI reason for x given F , and a minMAJ reason for x given F .

For computing prime-implicant explanations and min-MAJ reasons, we took advantage of the Pysat library (Ignatiev, Morgado, and Marques-Silva 2018) (version 0.1.6.dev15) that provides the implementation of the RC2 PARTIAL MAXSAT solver and an interface to MUSER [START_REF] Belov | MUSer2: An Efficient MUS Extractor[END_REF]. For majoritary reasons, we picked up uniformly at random 50 permutations of the literals describing the instance and tried to eliminate those literals (within the greedy algorithm) following the ordering corresponding to the permutation. We kept a smallest explanation among those derived (of course, the corresponding runtime that has been measured is the cumulated time over the 50 tries). Prime-implicant explanations have been computed using MUSes, as explained before.

We also derived a "LIME explanation" for each instance. Such an explanation has been inferred as follows. Given an input instance x, we first used LIME (Ribeiro, Singh, and Guestrin 2016) to generate a linear zero-threshold function w x ∈ R n . The value w x (z) of w x on any instance z is given by w x (z) = 1 if w ⊤ x z > 0, and w x (z) = 0 otherwise. Now, when x is classified positively by w x , in order to derive an explanation, it is enough to sum in a decreasing way the positive weights occurring in w x until this sum exceeds (the opposite of) the sum of all the negative weights occurring in w x . The term t composed of the variables x i associated with the positive weights which have been selected is, by construction, a minPI reason for x given w x since for every x ′ covered by t, the inequality w ⊤

x x ′ > 0 holds. Indeed, the inequality w ⊤

x x ′ > 0 holds in the worst situation when all the variables associated with a positive weight in w x and not belonging to t are set to 0, whilst all the variables associated with a negative weight in w x are set to 1 (see also [START_REF] Marques-Silva | Explaining Naive Bayes and Other Linear Classifiers with Polynomial Time and Delay[END_REF]). Instances that are classified negatively can be handled in a similar way.

Experiments have been conducted on a computer equipped with Intel(R) XEON E5-2637 CPU @ 3.5 GHz and 128 Gib of memory. A time-out (TO) of 600s has been considered for each instance and each type of explanation, except LIME ones (no time bound has been used for them).

Experimental Results. A first conclusion that can be drawn from our experiments is the intractability of computing minPI reasons in practice; this coheres with the complexity result given by Proposition 5. Thus, we have been able to compute within the time limit of 600s a minPI reason for only 10 instances and a single dataset (compas).

Due to space limitations, we report hereafter empirical results about three datasets only, namely placement, gisette and dorothea. The results obtained on the other datasets are similar and available on line. placement is a small dataset about the placement of 215 students in a campus; students are described using 13 features, related to their curricula, the type and work experience, and the salary. An instance is labelled positive when the student gets a job. The random forests consist of 25 trees, and their mean accuracy is 97.6%. gisette is much larger, including 5000 features and 7000 examples. Each feature is a pixel or a combination of pixels, and the task is to separate the digits 4 and 9. The random forests consist of 85 trees, and their mean accuracy is 96%. Finally, dorothea is a high-dimensional dataset, with 100,000 features and 1950 examples. Each instance is an organic molecule, and the goal is to discriminate binding compounds from non-binding ones. Here, the random forest consists of 71 trees, with a mean accuracy of 93%.

Figure 2 provides the results obtained for placement, using four plots. Each dot represents an instance. The first plot shows the time needed to compute a reason on the x-axis, and the size of this reason on the y-axis. On this plot, there are no dots for minPI reasons, because their computation did not terminate before the time-out. The plot also highlights that all other reasons have been computed within the time limit, and in general using a small amount of time. In particular, it shows that the direct reason can be quite large, that the computation of LIME explanations is usually more expensive than the ones of the other explanations, and that LIME explanations can be very short. 3 A box plot about the sizes of all the explanations is reported (the LIME ones and the direct reasons are not presented for the sake of readability). The figure also provides two scatter plots, aiming to compare the sizes of majoritary reasons with the sizes of prime-implicant explanations, and the sizes of the minMAJ reasons with the sizes of prime-implicant explanations.

These plots clearly show the benefits with respect to size reduction that can be offered by considering majoritary reasons and minMAJ reasons instead of prime-implicant explanations. At first sight, these empirical results may look surprising since, by construction, for any majoritary reason t for x given F (including the shortest ones) there exists at least one prime-implicant explanation for x given F that is implied by t (hence that cannot be larger). As to majoritary reasons, one must keep in mind that the result that is reported is a shortest reason out of a set of 50 majoritary reasons that are computed for each x (so to say, we leverage the tractability of computing such reasons to tackle the size issue). For minMAJ reasons, the PARTIAL MAXSAT algorithm used to compute them aims to minimize the size of the reason that is derived, while MUS algorithms for computing prime-implicant explanations do not focus on the size. Indeed, computing minimum-sized MUSes is much harder, as explained previously (see Proposition 5).

Figures 3 and4 synthesize the results obtained for gisette and dorothea, respectively. Conclusions similar to those drawn for placement can be derived for gisette and dorothea, with some exceptions. First of all, no dots have been drawn for minMAJ reasons because the computation of such reasons did not terminate within the time limit. Furthermore, LIME explanations are much longer. This can be partly explained by the fact that the computation achieved by LIME relies on a binary representation of the instance that is quite different (and possibly much larger) than the one considered in the representation of the random forest. Indeed, each decision tree in the forest focuses only on a subset of most important features (in the sense of Gini criterion) found during the learning phase. In our experiments, the size of LIME explanations was typically large for high-dimensional datasets.

When minMAJ reasons are difficult to calculate (as it is the case for gisette and dorothea), a natural approach is to look for approximations. From this perspective, we took advantage of the incremental PARTIAL MAXSAT algorithm called LMHS [START_REF] Saikko | LMHS: A SAT-IP Hybrid MaxSAT Solver[END_REF] to do the job. Specifically, the result given in Proposition 8 provides a way to derive abductive explanations for an instance x given a random forest F in an anytime fashion. Basically, using LMHS, a Boolean assignment z satisfying all the hard constraints of C hard and a given number, say k, of soft constraints from C soft is looked for (k is set to 0 at start). If such an assignment is found, then one looks for an assignment satisfying k + 1 soft constraint, and so on, until an optimal solution is found or a preset time bound is reached. By construction, every z that is generated that way is such that t x ∩ t z is an implicant of F that covers x (and hence, an abductive explanation). In practice, the approximation z of a minMAJ reason for x given F , which is obtained when the time limit is met, can be significantly shorter than the prime-implicant explanation for x given F that has been derived.

In our experiments, we used three time limits: 10s, 60s, and 600s. The results are reported in the box plots and the scatter plots in Figures 2, 3, and4. As illustrated by the box plots, the sizes of the approximations z which are derived gently decrease with time. The scatter plots indicate that significant size savings can be achieved even for the smallest time bound of 10s that has been considered.

Conclusion

We have introduced, studied, and evaluated new notions of abductive explanations suited to random forests, namely direct reasons, majoritary reasons and minimum-sized majoritary reasons. Unlike prime-implicant explanations, majoritary reasons and their minimum-sized counterparts may contain irrelevant features. Nevertheless, in practice, majoritary reasons and minMAJ reasons appear as valuable alternative to prime-implicant explanations. Indeed, majoritary reasons can be computed in polynomial time while prime-implicant explanations cannot (unless P = NP). In addition, in most of our experiments, majoritary reasons slightly smaller than prime-implicant explanations have been computed thanks to a simple greedy algorithm with random permutations of literals. minMAJ reasons can be looked for when majoritary reasons are too large, but this is at the cost of an extra computation time that can be important, and even prohibitive in some cases. However, minMAJ reasons can be approximated using an anytime PARTIAL MAXSAT algorithm. Empirically, approximations can be derived within a small amount of time and their sizes are significantly smaller than the ones of prime-implicant explanations.

Figure 2 :

 2 Figure 2: Empirical results for the placement dataset.

Figure 3 :

 3 Figure 3: Empirical results for the gisette dataset.

Figure 4 :

 4 Figure 4: Empirical results for the dorothea dataset.

Unlike(Ignatiev, Narodytska, and Marques-Silva

2019), we do not require those explanations to be minimal with respect to set inclusion, in order to keep the concept distinct (and actually more general) than the one of prime-implicant explanations.

Recall that LIME explanations are not guaranteed to be abductive. See also[START_REF] Narodytska | Assessing Heuristic Machine Learning Explanations with Model Counting[END_REF]) that reports some experiments about ANCHOR explanations.

Acknowledgements

Many thanks to the anonymous reviewers for their comments and insights. This work has benefited from the support of the AI Chair EXPEKCTATION (ANR-19-CHIA-0005-01) of the French National Research Agency (ANR). It was also partially supported by TAILOR, a project funded by EU Horizon 2020 research and innovation programme under GA No 952215.