J Réthoré

Lamcos

Automatic crack tip detection and stress intensity factors estimation of curved cracks from digital images

Keywords: curved crack, stress intensity factors, crack tip detection, finite element, Williams' series, digital image correlation

With the development of full-field measurement techniques, it has been possible to analyze crack propagation experimentally with an increasing level of robustness. However, the analysis of curved cracks is made difficult and almost unexplored because the possible analysis domain size decreases with crack curvature, leading to an increasing uncertainty level. This paper proposes a digital image correlation technique, augmented by an elastic regularization, combining finite element kinematics on an adapted mesh and a truncated Williams' expansion. Through the analysis of two examples, the proposed technique is shown to be able to address the experimental problems of crack tip detection and stress intensity factors estimation along a curved crack path.

INTRODUCTION

With the development of full-field measurement techniques, it has been possible to analyze experimentally crack propagation with an increasing level of robustness. The improvement of the robustness is directly related to the increasing amount of data available for the analysis. For crack tip fields analysis, this robustness basically depends on the size of the domain used for experimental stress intensity factors extraction.

The analysis of curved cracks is made difficult and almost unexplored because the possible extraction domain size, and thus the maximum robustness, decreases as the crack curvature is increasing. The experimental estimation of the crack tip position is also a difficult question in this context.

The most common practice to estimate stress intensity factors (SIF) at the tip of a crack relies on a post-processing of the displacement or strain/stress fields. Analyzing the displacement field leads to a least-squares projection onto an appropriate basis. This technique has been used either in numerical simulations e.g. [START_REF] Duarte | A generalized finite element method for the simulation of three-dimensional dynamic crack propagation[END_REF] or in the post-processing of experimental 2D or 3D full displacement field measurements e.g. [START_REF] Mcneill | Estimation of stress intensity factor by digital image correlation[END_REF][START_REF] Abanto-Bueno | Investigation of crack growth in functionally graded materials using digital image correlation[END_REF][START_REF] Yoneyama | Automatic evaluation of mixed-mode stress intensity factors utilizing digital image correlation[END_REF][START_REF] Roux | Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches[END_REF][START_REF] Lopez-Crespo | Study of a crack at a fastener hole by digital image correlation[END_REF][START_REF] Réthoré | Digital volume correlation analyses of synchrotron tomographic images[END_REF].

Extracting SIF from strain/stress fields is usually based on Linear Elastic Fracture Mechanics theory and path independent integrals. This kind of technique is extremely popular for numerical simulations and it has been used in a number of different cases in 2D or in 3D: quasi-static crack propagation, cohesive fracture, confined plasticity, fretting with contact and friction between the crack faces, dynamic crack propagation.

The interaction integral has also been used to extract SIF in 2D [START_REF] Réthoré | Estimation of mixed-mode stress intensity factors using digital image correlation and an interaction integral[END_REF] and in 3D [START_REF] Rannou | Three dimensional experimental and numerical analysis of a fatigue crack[END_REF] from experimental displacement fields measured by digital image or volume correlation.

However, due to the derivation of the displacement field to obtain the strain field, the noise sensitivity of this class of strain/stress based techniques is higher than for displacement based methods. Even if the noise sensitivity of the interaction integral can be optimized [START_REF] Réthoré | Noise-robust stress intensity factor determination from kinematic field measurements[END_REF], it has been shown in [START_REF] Réthoré | Optimal and noise-robust extraction of fracture mechanics parameters from kinematic measurements[END_REF] that displacement based techniques have lower noise sensitivity. The conclusion of this paper is that the best accuracy is obtained when SIF are directly estimated, i.e. when they are among the degrees of freedom of the problem to solve as in [START_REF] Roux | Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches[END_REF]. This remark applies not only for experimental displacement fields but also for numerical simulations for which different strategies have been proposed for direct estimation of SIF. These strategies have been developed in the context of the eXtended Finite Element Method (X-FEM) [START_REF] Liu | Xfem for direct evaluation of mixed mode sifs in homogeneous and bi-materials[END_REF][START_REF] Chahine | Crack-tip enrichment in the xfem method using a cut-off function[END_REF] and they recover optimal convergence in terms of SIF if the underlying partition of unity finite element functions are inactivated over the SIF extraction zone [START_REF] Xiao | Implementation of hybrid crack element on a general finite element mesh and in combination with xfem[END_REF][START_REF] Réthoré | Hybrid analytical and extended finite element method (HAX-FEM): A new enrichment procedure for cracked solids[END_REF][START_REF] Passieux | Direct estimation of generalized stress intensity factors using a three-scale concurrent multigrid x-fem[END_REF]. The latter two references make use of the coupling between an X-FEM model and an analytical basis based on a truncation of the Williams' series [START_REF] Williams | On the stress distribution at the base of a stationary crack[END_REF]. The coupling of the two models is achieved either by a partition of the energy over an overlapping domain [START_REF] Réthoré | Hybrid analytical and extended finite element method (HAX-FEM): A new enrichment procedure for cracked solids[END_REF] or with a Mortar method [START_REF] Passieux | Direct estimation of generalized stress intensity factors using a three-scale concurrent multigrid x-fem[END_REF]. These two methods were shown to be able to simulate the propagation of a crack under loading conditions leading to a curved crack path [START_REF] Réthoré | Mixed-mode crack propagation using a hybrid analytical and extended finite element method[END_REF][START_REF] Passieux | Local/global non-intrusive crack propagation simulation[END_REF].

The use of the Williams' series (obtained for a straight crack) to track curved cracks requires that the size of the extraction domain over which they are used as displacement basis functions is small compared to the local crack curvature. For numerical simulations, the performance of the direct SIF extraction is related to the ratio between the extraction domain size and the mesh size in the vicinity of the crack tip.

The mesh is thus refined to obtain an objective performance level given the extraction domain size appropriate to the crack curvature. From the experimental point of view, the existing technique for direct SIF estimation from digital images was proposed in [START_REF] Roux | Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches[END_REF].

In the latter, the optical flow equation is solved over a circular domain centered at the presupposed crack tip and Williams' series are used as displacement basis functions.

The performance of the SIF extraction in terms of accuracy and noise sensitivity is again linked to the extraction domain size which has to be as large as possible. This constraint is not compatible with the analysis of curved cracks for which the extraction domain has to be compared to the crack curvature. For this reason, the method proposed in [START_REF] Roux | Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches[END_REF] is not suitable for the experimental analysis of curved cracks. This is also the case for the methods cited above that can be used for extracting mixed-mode SIF, but not along curved crack.

Concerning the determination of the crack tip position from experimental displacement fields, the crack tip coordinates can be included in the least-squares minimization process between the experimental displacement field and the analytical basis [START_REF] Yoneyama | Automatic evaluation of mixed-mode stress intensity factors utilizing digital image correlation[END_REF]. However, a more robust approach relies on the concept of an elastic crack tip and its detection using the first super-singular term of the Williams' expansion [START_REF] Hammam | Stress intensity factor gauging by digital image correlation: Application in cyclic fatigue[END_REF]. This approach reveals robust with respect to a non-linear process zone due to plasticity [START_REF] Henninger | Enriched kinematic fields of cracked structures[END_REF],

diffuse damage [START_REF] Leplay | Identification of damage and cracking behaviours based on energy dissipation mode analysis in a quasi-brittle material using digital image correlation[END_REF] or cohesive fracture [START_REF] Réthoré | Identification of a cohesive zone model from digital images at the micron-scale[END_REF]. It has been shown to be robust enough for real time detection of crack propagation for load shedding in fatigue crack propagation experiments [START_REF] Durif | Controlling stress intensity factors during a fatigue crack propagation using digital image correlation and a load shedding procedure[END_REF]. This paper is aimed at proposing a methodology that allows for the analysis of curved cracks using digital images. It includes the detection of the crack path, the detection of the crack tip and the estimation of SIF and higher order terms like the T -stress.

The proposed method uses a truncated Williams' expansion as the displacement basis function in the crack tip vicinity. The direct estimation of the Williams' expansion terms is obtained through a condensation of the degrees of freedom of a finite element model using mesh adaptation at the crack tip. The issue of decreasing noise robustness due decreasing extraction domain size, required by the crack curvature, is solved by adding a supplementary elastic regularization term to the optical flow equation. The paper is organized as follows: first the methodology for direct SIF estimation is described, then the formulation of the optical flow equation with elastic regularization is detailed. Then, the implementation of the methodology is described. Last, two examples are analyzed.

The first example consists of the analysis of artificially generated images and the second example considers an experiment on a sample geometry and loading conditions leading to a curved crack propagation.

DIRECT SIF ESTIMATION

For the direct estimation of SIF, a truncation of the Williams' expansion is used in combination with a finite element model. This section details how these two descriptions are combined to obtain a direct SIF estimation. These descriptions require that two coordinate systems are defined (see Figure 1):

• the global frame (X,Y ) associated with the structure • the local frame at the crack tip (x, y) aligned with the crack

In the following, we adopt a finite element discretization for the displacement field. The discretized displacement field reads:

u u u FE (X) = ∑ k=1...2N u k N N N k (X) = [N] {U}, (1)
where N is number of nodes, leading to 2N degress of freedom (DOF), and N N N k the finite element shape functions. Odd, respectively even, k index corresponds to a DOF along X, respectively Y . A compact notation is adopted in the following: to collect the values of the displacement at the evaluation points, a matrix form of this equation is also used:

[N] {U}, where [N] is a matrix that collects the value of the 2N shape functions at the evaluation points and {U} the vector that collects the values of the displacement DOF.

Note that [N] has 2N columns and as many lines as points where the shape functions are evaluated. For clarity, the positions of these evaluation points are omitted by the notation.

Williams' series

For a 2D elastic body with a semi-infinite straight crack subjected to a mechanical load, an analytical solution for the displacement field was proposed by [START_REF] Williams | On the stress distribution at the base of a stationary crack[END_REF]. The displacement is decomposed over a double series of φ φ φ n j solution fields that satisfy stress-free conditions along the crack faces. Two fracture modes (indexed by j) are usually distinguished: mode I (opening mode) and mode II (shear mode). These elementary solutions are indexed by n, and are written as follows:

φ φ φ n I (r, θ) = r n/2 κe inθ/2 - n 2 e i(4-n)θ/2 + ( n 2 + (-1) n )e -inθ/2 (2)
and

φ φ φ n II (r, θ) = ir n/2 κe inθ/2 + n 2 e i(4-n)θ/2 -( n 2 -(-1) n )e -inθ/2 . (3)
In these equations, i is the pure imaginary number, r the distance to the crack tip and θ the angle with respect to the coordinate system (x, y) aligned with the crack. Hence, adopting complex notations, the position of a material point is defined by z which can be expressed as x + iy or by r e iθ . κ is the Kolossov's constant, namely, κ = (3ν)/(1 + ν)

for plane stress or κ = (3 -4ν) for plane strain conditions, where ν is the Poisson's ratio.

This family of solutions is the appropriate basis to describe the displacement field around the tip of a crack. Usually, only n ≥ 0 terms are considered because n < 0 terms exhibit an infinite strain energy density at the crack tip. An actual displacement field u u u w , determined over the previously introduced basis is:

u u u w = ∑ j=I,II ∑ n=n m ...n M a n j φ φ φ n j (z) . (4)
In this equation, n m and n M define the truncation of the Williams' expansion that is considered and φ φ φ n j the projection of φ φ φ n j onto the X,Y coordinate system.

Condensation

On the sample surface, the displacement field is described using a finite element model.

Within the extraction domain Ω K surrounding the crack tip, a condensation technique is used in order to reduce the finite element kinematics to the truncated Williams'

expansion. The condensation is written as a linear operator between the finite element DOF of the entire mesh U and the reduced basis U r :

{U} = [L] {U r }. (5)
For the nodes outside the extraction domain, the condensation operator L reduces to a diagonal matrix. For the nodes within the extraction domain, the displacement is given by the φ φ φ n j functions evaluated at the corresponding node position. The expression of the condensation operator L is thus given by

L =                    1 . . . 1 φ φ φ n m I (r i o , θ i o ) . . . φ φ φ n M I (r i o , θ i o ) φ φ φ n m II (r i o , θ i o ) . . . φ φ φ n M II (r i o , θ i o ) . . . . . . . . . . . . . . . . . . φ φ φ n m I (r N , θ N ) . . . φ φ φ n M I (r N , θ N ) φ φ φ n m II (r N , θ N ) . . . φ φ φ n M II (r N , θ N )                    . (6)
In this equation, i o is the index of the first node in the extraction domain which contains all the nodes with an index between i o and N. Using this description, the 2(i o -1) first degrees of freedom in U r are standard finite element degrees of freedom but the 2N w last degrees of freedom correspond to the amplitude of the N W functions for mode I and mode II:

U r =                                u 1 . . . u 2(i o -1)
a n m I . . .

a n M I a n m II . . . a n M II                                . ( 7
)

Mesh adaptation

The Williams' expansion is written in a polar coordinate system. It is thus convenient to use this description over a circular domain. The finite element mesh is thus adapted to have an explicit description of a circle of radius r tip centered at the crack tip (see to be a good compromise between accuracy and computational cost.

DIGITAL IMAGE CORRELATION

To follow crack propagation, a series of digital images is captured and then analyzed.

However, at each step (image), the analysis reduces to the study of a pair of images: f , the reference image associated with the reference state of the structure, and g namely the deformed image associated with the deformed state. A matrix-vector format is also adopted for the images data. For an image f , F is a column vector that contains the grey level of each pixel in the area covered by the mesh. Within this formalism, the pixel integration (discrete summation over the pixels) simply writes as a scalar product between two column vectors.

Optical flow

The optical flow equation states that the grey level variation between the reference state f and the deformed state g is due to only the passive (no source term) advection of image texture ( the grey level is an intrinsic property of a material point):

f (X p ) = g(X p + u(X p )), (8)
X p giving the position of a pixel in the global coordinate system. A non-linear leastsquares resolution scheme of this equation adopted in the following [START_REF] Besnard | finite-element' displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF]. After the finite element discretization of the displacement has been introduced and the equation is linearized with respect to the solution increment dU = U -U p , where U p is the solution at the previous iteration, a linear system is obtained:

[M] {dU} = {b}. (9)
In this equation, the following expression is adopted for M:

[M] = [N] T [∇F] [∇F] [N] , (10)
and b is given by

{b} = [N] T [∇F] {F -G} (11)
In these equation ∇F is the vector containing the values of the image gradient , computed by a central finite difference method, at each pixel.

Elastic regularization

To compensate for the use of a small extraction domain, an elastic regularization of the optical flow is adopted. Let K be the elastic finite element stiffness matrix calculated for the adapted mesh, then Equation ( 9) is augmented by adding a weighted equilibrium gap energy [START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis of fractured samples[END_REF] [

M] + ω( c ) [K] T [P] T [P] [K] {U} = {b} + [M] {U p }. (12)
In this equation, ω( c ) is a scalar parameter weighting the regularization. It has been shown that the regularization introduces a mechanical low-pass filter which cut-off wave length c can be controlled [START_REF] Leclerc | Digital volume correlation: what are the limits to the spatial resolution[END_REF]. P is a selection matrix that allows to select the DOF to be controlled by the elastic regularization. It formally reads as a diagonal matrix with 1 if the corresponding degree of freedom is regularized, and 0 if not. In the present case the degrees of freedom hold by the nodes inside the extraction domain Ω K are not regularized as they are constrained by the condensation operator L that intrinsically prescribes divergence free elastic stress. The nodes on the boundary of the structure that are submitted to non stress free conditions are not regularized as well (see Figure 1). This system is solved very efficiently with a preconditioned conjugate gradient solver using U p as an initial solution. The convergence is typically reached within 3 to 5 iterations.

Using the condensation, the resolution algorithm is written in term of U r by introducing Equation (5) into Equation ( 12):

[L] T [M] + ω( c ) [K] T [P] T [P] [K] [L] {U r } = [L] T {b} + [L] T [M] [L] {U p r }. ( 13
)
After convergence is reached, the resolution directly provides estimates of the selected Williams' amplitudes a n j . They are then scaled to obtain the meaningful quantities as SIF or T -stress:

K I = 2µ √ 2π a 1 I , K II = 2µ √ 2π a 1 II , T = 2µ √ 2π a 2 I , ( 14
)
where µ is the first Lamé constant given by µ = E 2(1+ν) , E being the Young's modulus.

IMPLEMENTATION

In this section, all the steps required to obtain the crack path, crack tip position and SIF estimates are detailed. A graphical overview of the methodology is presented in Figure 2.

Crack path detection

As proposed in [START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF], the correlation residual map η

= | f (X p ) -g(X p + u(X p ))| allows
for a robust positioning of the crack surface. This was used in 2D but also in 3D e.g. [START_REF] Rannou | Three dimensional experimental and numerical analysis of a fatigue crack[END_REF][START_REF] Réthoré | Digital volume correlation analyses of synchrotron tomographic images[END_REF]. The crack surface is defined from the residual map either manually or automatically by tracking the surface of maximum error [START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF][START_REF] Rannou | Three dimensional experimental and numerical analysis of a fatigue crack[END_REF] or using a virtual image correlation technique [START_REF] Réthoré | Curve and boundaries measurement using b-splines and virtual images[END_REF]30]. To get a residual map η suitable for using one of these detection algorithms, a displacement measurement analysis is run for all the images of the experiment using a non-conforming regular mesh. Since the displacement discontinuity at the crack faces are not described, large error levels are obtained on the crack faces.

After the crack surface detection is performed, either for a set of points or a level set function, a mesh adapted to the sample boundaries and the crack faces is generated. It also fulfills the requirements for direct SIF estimation using the condensation method proposed above (explicit description of a circle of radius r tip to define the extraction domain and a circle of radius r int for element removal).

Crack tip detection

Once the crack surface is defined, the crack tip position is tracked along the crack support. The image sequence is analyzed starting from the last state using the adapted mesh and direct SIF estimation. The Williams' expansion is truncated between n m = -3

and n M = 7. Super-singular terms are used following [START_REF] Hammam | Stress intensity factor gauging by digital image correlation: Application in cyclic fatigue[END_REF]. The n = -1 term can be used to detect the equivalent elastic crack tip position. For this purpose, the following property of Williams' functions is used:

∂φ φ φ n j ∂x = n 2 φ φ φ n-2 j . (15)
Then, suppose that the estimated crack tip position is shifted by a small distance d. If the displacement u u u w is expressed as a function of the Williams' series in the shifted

coordinate system (z = z + d) u u u w = ∑ n≥0, j a n j φ φ φ n j (z) = ∑ n, j ãn j φ φ φ n j (z) , (16)
using a first order Taylor expansion of φ φ φ n j (z) and Equation ( 15) allows consideration of not only n ≥ 0 terms but also n < 0 terms. Then, one can write the following by considering that the expression of u u u w from Equation 4 has a zero term for n = -1 and mode I:

- d 2 ã1 I + ã-1 I = 0 . ( 17
)
The crack tip position along the crack faces and coordinate system is progressively corrected until sufficient precision (for the crack tip position) is reached. Because the first order correction only is considered in order to obtain Equation 17, the relation between the actual shift, d, and its estimation from -2 ã-1 I / ã1 I is linear. In practice, this linear dependence is observed to occur for quite a large range of crack tip shifts and convergence is usually obtained in 2 or 3 iterations.

To avoid expensive computational cost due to remeshing, a mesh morphing technique is used to update the mesh between two iterations of the crack tip detection. An elastic finite element simulation, with 0 Poisson's ratio, is run under the following boundary conditions:

• all the nodes inside the extraction domain have a prescribed displacement of d along the local tangent to the crack tip x,

• the displacement of the nodes on the crack faces (except those inside the extraction domain) is along the local tangent to the crack surface only,

• the non-regularized nodes on the mesh boundary are fixed,

• the displacement of all the other nodes on the mesh boundary is along the local tangent to the mesh boundary only.

Once convergence is reached on the crack tip position, i.e. the correction d being lower than 1 pixel in amplitude, the values of the meaningful Williams' amplitudes are stored as well as the position of the crack tip. Then, a new mesh is generated and the previous image is analyzed. This avoids large mesh distortion and preserves computational cost.

The performance of the proposed methodology will now be evaluated for two examples. Especially, it will be assessed that small extraction domains can be used with the proposed strategy because the elastic regularization preserves from increasing uncertainty on the SIF and estimated crack tip position.

EXAMPLES

The method proposed in the present paper is now tested with two examples. The first example is based on the analysis of artificially deformed images. As the crack tip position and SIF amplitudes are prescribed, this example allows for a quantitative evaluation of the performance of the method. Then the ability of the method to track curved crack paths is illustrated through the analysis of an experiment.

Synthetic images analysis

Starting from a real 128 × 128 pixels image duplicated four times in each direction (see • the RMS relative SIF error that is calculated by taking the standard deviation (RMS) over the propagation steps of the relative gap in SIF (the scatter between the estimated value and the prescribed one normalized by the prescribed value)

• the RMS crack tip position error that is calculated by taking the standard deviation over the propagation steps of the gap between the estimated crack tip position and the prescribed one.

These indicators allow for quantifying the uncertainty on these two output variables (SIF and crack tip position) of the analysis. The evolution of these indicators with respect to c and r tip is detailed below.

Influence of r tip For these analyses, c is set to 400 pixels but r tip is varying between 30 and 100 pixels. Analyses are also conducted with the integrated approach in [START_REF] Roux | Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches[END_REF]. They are labelled as W-DIC. These analyses were carried out to show how the uncertainty increases when the extraction radius is decreasing. Indeed, in Figure 5, it is observed that the uncertainty of W-DIC analyses exponentially increases when r tip is decreasing.

Conversely, thanks to the elastic regularization, the extraction radius does not influence the results obtained by the proposed method. The averaged RMS relative SIF error is about 1.5 × 10 -3 and the average RMS crack tip position error is about 0.5 pixel. Note that the convergence criterion for crack tip detection is set to 1 pixel. The error obtained for these analyzes on synthetic images should have been lower with a lower convergence criterion. However, a value lower than 1 pixel may be unrealistic for the real experimental case.

This first series of analyses shows that a very small extraction domain size can be used without deteriorating the performances of the method. Indeed, the elastic regularization allows the outer domain (outside the extraction domain) to share far-field information with the analytical model whatever the extraction domain size, thus keeping the uncertainty levels very low. The key parameter thus seems to be the cut-off wave length c .

Influence of c

The extraction domain size is now kept at 40 pixels and c is varying from 0 to 400 pixels. The case c = 0 pixel corresponds to the integrated method W-DIC [START_REF] Roux | Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches[END_REF]. The evolution of the error indicators is shown in Figure 6. Concerning the crack tip positioning, it is obtained that the elastic regularization is very efficient, meaning that even for low values of c (from 10 pixels), the uncertainty is lower than or about 1 pixel.

For the SIF, the evolution of uncertainty with the cut-off wave length is also sharp and an improvement of the sensitivity is obtained when increasing c . A fit, on a power law, in a logarithmic scale gives a slope of -1 meaning that the SIF uncertainty behaves as

1/ c .
From these two series of analyzes, it is concluded that the proposed method with direct SIF estimation by condensation and elastic regularization allows for estimating SIF and crack tip position even for small extraction domains. The extraction domain can be as small as required by the crack curvature provided there is a layer of about 10 elements along its radius and the cut-off wave length of the elastic regularization is about the extraction domain size.

Mixed mode crack propagation experiment

An experiment is conducted using a conventional tensile electromechanical device. The geometry of the sample, inspired from [START_REF] Grégoire | Dynamic crack propagation under mixed-mode loading -comparison between experiments and x-fem simulations[END_REF], is a rectangular plate of 140 mm×70 mm with a hole of 30 mm in diameter, see Figure 7. The hole center is aligned along the transverse (Y ) direction and it is shifted by 25 mm from the specimen center along X direction. The specimen thickness is 10 mm, leading to an approximation of the 3D state by a 2D plane strain hypothesis. The material is a commercial PolyMethyl Methacrylate (PMMA) (Perspex). Its elastic properties, within the strain rate range of the quasi-static experiment conducted herein, are taken to E = 5GPa for the Young's modulus and ν = 0.32 for the Poisson's ratio.

A pre-notch is machined manually before the specimen is submitted to a compressive load along X at a constant speed of 0.1 mm/min. The pre-notch of about 35 mm length from the hole center, with an angle of 30 o with respect to the X axis, is chosen for generating mixed-mode conditions at initiation and growth. It is expected that the crack runs along X and progressively reaches the specimen axis, leading to a pure mode I situation.

During the experiment, grey scale images of 4872 × 3248 pixels, with 8-bit digitization, are acquired. The camera is mounted with a 200 mm Nikon lens, resulting in a physical pixel size of 32.5 µm. The sample surface is prepared using white and black spray paints, see Figure 7. 113 images are captured before complete failure. Crack initiation occurs at images number 27. In the results presented below, the time scale is shifted so that crack initiation occurs at t=0. The time is in fact a virtual time standing for the image number.

First, an analysis with a regular mesh of quadrangular elements of 64-pixel side is performed. The two components of the displacement field for the last image are displayed in Figure 8. The discontinuity of the displacement field across the crack faces is clearly evidenced for the two displacement components, unlike it is smoothed by the regular finite element mesh. Crack openings of about 20 pixels are obtained at this loading step. From this analysis, an correlation error map η is computed and then used to determine the crack path as shown in initiation and it changes sign during propagation. However, the values for K II are very low in amplitude and this is difficult to interpret. The values of K Ieq just after initiation are slightly lower than the value obtained at initiation. This is due to the finite thickness of the pre-notch that leads to a lower singularity and thus an estimated value for K Ieq at initiation higher than the one produced by a sharp crack. Then the value of K Ieq is slightly decreasing as the crack propagates. This should be linked to the variation in the crack tip speed as obtained in Figure 16 but also to the variation of the T -stress induced by structural effects as observed in Figure 15.

CONCLUSION

In this paper, a methodology is proposed to experimentally estimate crack tip position and SIF along curved crack paths using digital images. The method is based on a digital image correlation technique using finite elements with mesh adaptation and a truncated i in the proposed strategy, there is no need for boundary conditions other than stress free conditions. The biases that may be introduced by considering idealized, Dirichlet or non-zero Neumann, conditions are thus avoided. This problem was already emphasized in the literature for 3D stationary cracks [START_REF] Roux | Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2d and 3d cracks[END_REF] or for the identification of material parameters [START_REF] Réthoré | A fully integrated noise robust strategy for the identification of constitutive laws from digital images[END_REF].

ii performing a finite element simulation of crack propagation requires a criterion for the crack angle and a crack growth law. Concerning the crack angle, the max hoop stress criterion has shown its ability to capture the crack bifurcation in a brittle material such PMMA. Common pratices for crack growth are to suppose that the crack propagates at a constant value of K Ieq , or to adopt a Paris' like regularization of the Griffith's law. It is clearly shown by the results of this paper that, for the analyzed case, the relation between K Ieq and the crack tip velocity is much more complex. In Figure 15, it is observed that K Ieq starts decreasing from step 50 whereas it is observed on Figure 16 that the crack is accelerating from this step. This point was not emphasized in this paper for which the focus is kept on the technical aspects. However, the results presented (which are probably among the first ones for a curved crack) migh open a discussion e.g. on the effect of the crack curvature or of the T -stress on the crack velocity.

  It contains: n = 0 rigid body translations n = 1 usual asymptotic terms that are weighted by stress intensity factors n = 2 rigid body in-plane rotation and T -stress n > 2 sub-singular terms, i.e. the corresponding elastic stress field is not singular at the crack tip.

Figure 1 )

 1 Figure 1). The nodes along the crack path are split to allow displacement discontinuity across the crack faces. For using the first super-singular Williams' component to detect the crack tip position along the crack faces (see Section 4.2), it is required that a small disc around the crack tip is discarded from the analysis. The finite element mesh thus also explicitly describes a circle of radius r int . The elements inside this circle are removed from the initial finite element mesh. The element size is refined in the extraction domain Ω K to control the interpolation error between the finite element description and the analytical Williams' functions in the condensation operator. A factor of 10 between the element size and the radius of the extraction domain r tip has shown, in the present work,

Figure 3 (

 3 Figure 3(a)), a series of 7 1024 × 1024 pixels images is generated. For generating these images, the same sub-pixel grey level interpolation as for digital image correlation is used, i.e. cubic splineinterpolation. The displacement field used to deforme the reference image consists of a mode I singular field, n = 1. The amplitude a 1 I of the prefactor of

Figure 9 .Figure 12 , 5 . 2 . 2 .

 912522 Figure 12, negative values of K I are reported. These negative values correspond to a globally compressive loading of the specimen and they are allowed by the non-zero finite thickness of the machined pre-notch. Indeed, crack faces are allowed to close due to the pre-notch thickness resulting in negative K I values. The equivalent mode I SIF K Ieq , in the sense of the maximum hoop stress, is also plotted in Figure12. A values of K Ieq at initiation is 2.23 MPa.m 1/2 . After this initiation phase, the crack starts propagating with a kink due to the high mode mixity of the loading. The propagation phase is analyzed below.

  Williams' expansion. The two descriptions are coupled through a condensation of the finite element degrees of freedom onto the Williams' coefficients. The analysis of a curved crack path requires a small extraction domain compared to the crack curvature.The decrease of the robustness of the extraction of the crack tip position and the SIF is compensated herein by using an elastic regularization over the finite element adapted mesh. It is shown by the first example in the paper, that the robustness of the proposed methodology does not depend on the extraction domain size, which can thus be chosen as small as required by the crack curvature, but on the cut-off wave length of the elastic regularization. However, a cut-off wave length of the order of the extraction domain size ensures robust extraction of crack tip position and SIF. The ability of the proposed methodology to track curved cracks is then illustrated through the analysis of an experiment. Initiation and crack propagation are analyzed and robust estimates of SIF, T -stress and crack length are obtained. A fast and complete implementation of the method is also detailed in the paper. Compared to direct finite element simulations of the experiment, the approach proposed in the paper offers the following capabilities:

Figure 1 .Figure 2 . 10 Figure 3 .

 12103 Figure 1. Shematic representation of a sample with a crack. Definition of: the global coordinate system (X,Y ), the local coordinate system (x, y), the extraction domain Ω K and its radius r tip and the radius r int of the zone of removed elements at the crack tip.

Figure 4 .

 4 Figure 4. SIF (a) and crack tip position (b) obtained for r tip = 40, h = 4 and c = 400, compared to the prescribed solution.

Figure 5 .

 5 Figure 5. Relative RMS SIF error and RMS crack tip position error as a function of the extractiondomain size for the proposed method and the integrated method (W-DIC) in[START_REF] Roux | Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches[END_REF].

Figure 6 .

 6 Figure 6. Relative RMS SIF error and RMS crack tip position error as a function of the cut-off wave length c of the elastic regularization.

Figure 8 . 10 Figure 9 .Figure 10 .Figure 12 . 10 Figure 13 .Figure 14 .Figure 15 .Figure 16 .

 810910121013141516 Figure 8. Horizontal and vertical displacement map (in pixel) for the last image and a regular finite element mesh with 64-pixel elements.


P Fogel
email: paulfogel@hotmail.com

C Geissler
email: cgeissler@advestis.com

P Cotte
email: pcotte@advestis.com

G Luta

Applying separative non-negative matrix factorization to extra-financial data

Keywords: Interpretability, Features Engineering, ESG data, Clustering, Dimension reduction, Machine Learning

We present here an original application of the non-negative matrix factorization (NMF) method, for the case of extra-financial data. These data are subject to high correlations between co-variables, as well as between observations. NMF provides a much more relevant clustering of co-variables and observations than a simple principal component analysis (PCA). In addition, we show that an initial data separation step before applying NMF further improves the quality of the clustering.

Introduction

A question frequently encountered in economics is that of business classification. The financial industry (in particular index providers) has developed several standards for classifying companies based on their dominant industrial activity, such as the Global Industry Classification Sector (GICS R ). These classifications are widely used by investment fund managers as a tool for diversifying economic risk, or alternatively to create thematic portfolios focused on one or a small number of industry sectors.

The purpose of this paper is to outline a method for classifying companies based on their extra-financial profiles. To this end, we will use scores from extra-financial rating agencies. These scores quantify the behaviour of issuers in the three main pillars of Environment (E), Social (S) and Governance (G).

For the data-scientist, these data constitute a rich and complex field of study. The richness is linked to the variety of themes noted, ranging from the treatment of employees and suppliers, to the carbon footprint and the independence of the Board of Directors. These ratings provide information that cannot be derived from the industrial sector alone, and provide a new viewpoint to discriminate between companies.

The complexity stems from the fact that those data are not provided uniformly for all issuers, as some specific scores only make sense for very specific sectors. On the other hand, behavioural scores are both numerous (typically between 50 and 100) and sometimes highly correlated.

We describe here an unsupervised clustering experiment based on a matrix of extra-financial data. The data constituting the matrix are normalized between 0 and 100, and therefore lend themselves well to the use of non-negative matrix factorization (NMF) methods. This factorization provides a grouping of both co-variables and companies. The results are compared with those of a Principal Component Analysis (PCA). We also show, following a path initiated by [START_REF] Fogel | Applications of a novel clustering approach using non-negative matrix factorization to environmental research in public health[END_REF], that a tensor variant of NMF can be used to isolate groups of issuers characterized by either high or low scores on some co-variables, while ordinary NMF would mainly focus on groups of high scores.

Presentation of the data

We use ESG scores from the research of Vigeo-Eiris, whom we thank here for making these data available for this study.

The data consist of 77 scores spread over the three pillars E, S and G. They concern 699 European actions. The scores are positive numbers between 0 and 100, with 0 representing the worst possible score and 100 the best. These quantities do not have absolute significance but are essentially used to establish a ranking between issuers.

Nomenclature of variables

We will denote by X ESG the matrix formed by the 77 variables or 'features' as columns. This matrix is of type (699, 77). The generic term X ESG i,j is the value of the feature j for stock i.

Vigeo's behavioural scores are grouped into 6 areas, with taxonomy being reflected in the denomination of variables:

• Environment: variables names starting with env .

• Corporate Governance: " " cg • Social: Business Behaviour: " " c s

• Social: Human rights: " " hrts • Social: Human Resources: " " hr • Social: Community involvement: " " cin More information about the scoring methodology can be found at: https://vigeo-eiris.com/about-us-2/ methodology-quality-assurance/.

The interesting feature here is that we have a prior fundamental knowledge about variables that will be used as a benchmark to evaluate the relevance of any NMF-based clustering. Similarly, the fundamental information about the industrial sector the companies belong to will also be used as a prior knowledge.

Missing data

Variables are scores that rate the company's activity according to different criteria. Some of these criteria have no relevance outside a particular industry sector, and are therefore not reported for all issuers. The next two charts give a visual summary of the average fill rates grouped by sector (Figure 2) and by features (Figure 3). The sectors are presented along lines, by decreasing rate, while the variables are presented in columns, also by decreasing rate. As shown on chart 2, the consumer staples sector has the least amount of missing data. In contrast, several environment-related scores in the finance or telecommunications sector have little or no information at all, i.e. a fill rate close to zero. twenty two features out of 77 have a fill rate less than 75 %, the average being 81.1 %. Before performing any statistical analyses such as linear regressions, it is necessary to apply a missing values imputation strategy. A default strategy is to fill in the missing data for a variable with the mean or median, taken over the entire sample, of the non-missing values of that variable.

A first analysis using Principal Components Analysis (PCA)

Principal component analysis of a X data matrix focuses on the uncorrelated directions that explain the highest proportion of variance in the data. The directions are linear combinations of the original variables, which are mutually orthogonal. These principal directions are none other than the eigenvectors of the M = X t X matrix.

In the case of the X ESG matrix, we need to fill in the missing data with non-empty numerical values, because PCA algorithms are not meant to handle missing values. To do this, we apply the most neutral imputation strategy with respect to the objective of explaining variance: for each variable, a missing observation is replaced by the average of the variable, calculated on the set of non-missing data.

A first look at the results of a standard PCA allows the following observations to be made:

• The observation of the cumulated variance explained according to the rank of the component shows that it increases rather slowly. At least twenty two components are needed to explain 90 % of the total variance (Figure 4).

• The first component has coefficients that are all positive except one. All other components, which are mutually orthogonal, have positive and negative coefficients, at least half of which are greater in absolute value than 10 % of the largest coefficient. (Figure 5).

Except for the first component (thick blue line) for which all coefficients but one are positive, all other components have coefficients of both signs, which makes their interpretability difficult. 4 Factoring the data through Non negative Matrix Factorization (NMF)

Principles

Let us assume the available data are represented by an X matrix of type n × f , i.e. n rows and f columns.

We also assume that these elements are positive or null and bounded. With that in mind, we can re-write X (here shown as a 3 × 4 example matrix, with rows and columns indices starting at 0) as

  x 0 0 x 1 0 x 2 0 x 3 0 x 0 1 x 1 1 x 2 1 x 3 1 x 0 2 x 1 2 x 2 2 x 3 2   X =   w 0 0 w 1 0 w 0 1 w 1 1 w 0 2 w 1 2   W × h 0 0 h 0 1 h 0 2 h 0 3 h 1 0 h 1 1 h 1 2 h 1 3 t H +   0 0 1 0 2 0 3 0 0 1 1 1 2 1 3 1 0 2 1 2 2 2 3 2   =X-W × t H (1)
A non-negative factorization of X (Lee D. [1999]) then solves the following problem: Find two matrices W of type n × c and H of type f × c minimizing the quadratic norm of the last term of the right hand of Equation ( 1)

|| || 2 2 := X -W × t H 2 2
(2)

H and W being subject to the constraints of only having non-negative entries.

A few remarks:

• The quadratic norm used here is restricted to numeric entries (i.e. represents the square root of the sum of squared non-missing values). This slight extension of the definition allows to deal with input matrices having some missing values.

• The objective function in 2 may also include additional by L1 and L2 penalization terms applied to the W and H matrices. The aim of this penalization is to eliminate noise coming from small and generally unstable coefficients, and to concentrate the signal on the most significant contributors. This objective is called sparsity.

• This decomposition is an approximation, not an equality. Any matrices W and H that are solutions to (2) minimize the quadratic error of the approximation.

• It is common in the literature to refer to rows of W matrix as scores, and to those of H as loadings.

• The decomposition is not necessarily unique, even with the positivity constraint. One can indeed use any D matrix having non-negative coefficients, as well as its inverse D -1 , and then obtain identical decompositions of the form:

X = W DD -1t H + .
Different examples of such a transformation would be:

-A change of unit in the variables, with a diagonal D matrix representing the equivalence coefficient between two measuring systems.

-A permutation in the components, D being then represented by a permutation matrix.

• c ≤ min(n, f ) is an integer representing the number of selected components, also called the rank of the decomposition. Like in most clustering methods, this parameter represents a degree of freedom; determining the most adequate number of components requires an additional optimality criterion than the one for error, as the error generally decreases with the rank.

• The n rows of X will be referred to as observations, the f columns of X as variables or features, and the c columns of W as meta-variables or meta-features.

It is customary to consider the rows of X as n observations of objects from the same category (e.g. customers, patients, companies, etc). Each object has f attributes, and the data contained in the matrix correspond to the values taken by the attributes for the corresponding objects.

The main effect of this decomposition is to decrease the information necessary to describe an observation. The original observations of the X matrix can be recovered, at the cost of an approximation, as linear combinations with positive coefficients of the rows of the H matrix. A reduction in size has thus been achieved: part of the information contained in the original matrix, comprising f features per individual, can be approximately summarized by a smaller set of c features per individual.

A first NMF run

We start with a first NMF decomposition on the 500 × 77 scores matrix X. Two prior steps are handling missing values in X, and choosing a rank c.

We use the Python package nmtf developed by one of the authors. This package is public and available at https://pypi.org/project/nmtf/. It can be directly installed with the pip command: pip install nmtf.

The optimisation algorithm used in this package takes advantage of the extended definition of the quadratic norm and can handle missing values in input matrices, as opposed to the NMF version 0.17 of Scikit-Learn, which raises an error in case of missing entries in the matrix.

This ability to handle missing values avoids having to deal with the delicate issue of imputation of these missing values. In the field of extra-financial scores, a missing value is not simply data for which a measurement is missing because it did not take place. It can be data whose measurement simply does not make sense. To give a concrete example from the environmental domain, scores for waste reprocessing or contribution to deforestation are not relevant for some industrial sectors.

Another difference between nmtf and Scikit-Learn lies in the way sparsity is specified. In Scikit-Learn, the objective function is modified by addition of penalty coefficients that have an indirect effect on the sparsity of the W and H matrices. In nmtf, on the contrary, the sparsity objective on H or on W is provided as a parameter to the model, as a positive coefficient less than 1. For example, a parameter value of 0.9 on H will provide, if the resolution is possible, an H matrix with about 90% of zero entries.

The rank is set here at 6. This choice is not totally arbitrary as a first trial, it is exactly the number of topics set by Vigeo-Eiris for the taxonomy of the scores.

The solution X provided by the algorithm satisfies:

X -W × t H 2 2 = 5.7 × 10 -2 X 2 2 .
In other words, the squared approximation error for rank c = 6 is 5.7 % of the squared norm of X.

Referring to the performances of PCA on chart 4, it is necessary to aggregate the first 32 factors in order to explain a fraction of 100 % -5.7 % = 94.3 % of the total variance. This illustrates one of the main strengths of NMF: relaxing the requirement of factor orthogonality provides in return a much stronger explanatory power.

Clustering observations and features

An interesting induced property of NMF is the ability to cluster variables and observations in a natural way.

Regarding the clustering of features, the principle is to link each of the f variables with the meta-variable (indexed between 1 and c) having the highest coefficient on this variable. Similarly, each observation is assigned the index of the meta-variable with the highest contribution to the observation.

This grouping scheme has the disadvantage of being dependent on the particular system of units retained for the meta-variables, i.e. a transformation by an arbitrary D matrix. In order to get around this disadvantage, the authors of [5] introduce the notion of leverage, which represents "the ability of a row or column to specifically exert an influence on a particular element without a corresponding increase in the influence on other components". Leverage is a dimensionless quantity between 0 and 1. The clustering scheme based on maximum leverage is more stable than that based simply on the maximum coefficient.

The most visible benefit that we obtain from the non-negative matrix is the replacement of the initial 77 features with a reduced subset of 6 meta-features; the initial features can themselves be written as positive linear combinations of these meta-features.

Figure 6 shows the distribution of the coefficients (the 'loadings') of the initial 77 features over the 6 new meta-features, sorted by decreasing order. We observe that as the rank of the meta-features increases, fewer and fewer of the original features have significant coefficients on them.

Introducing Separative NMF

One of the few weaknesses of NMF is its asymmetry of treatment between high and low values. By construction, standard NMF will group observations by looking for the meta-features having highest scores. The attribution of a cluster will be more discriminant for companies characterized by a high score on a meta-feature (e.g. high environmental scores), and less discriminant for companies characterized by a low score. But when it comes to ESG scores, knowing that a company has a low score on a meta-feature is just as important as the information provided by a high score.

In [START_REF] Fogel | Applications of a novel clustering approach using non-negative matrix factorization to environmental research in public health[END_REF], an interesting variant of NMF is proposed in order to go around this obstacle. The method is called PosNegNMF. The idea is to separate the original matrix and to rewrite it as the sum of three components:

X = X + -X -+ X 0 (3)
with the following notations:

• X 0 is the "baseline", or the matrix of columns-wise medians of X.

• X + = M ax(0, X -X 0 ), the matrix of positive high deviations from X 0 .

• X -= M ax(0, X 0 -X), the matrix of positive low deviations from X 0 .

Having done this, PosNegNMF will perform a single decomposition of the matrix containing positive and negative deviations. This is equivalent to two independent decompositions, yielding two H matrices with a common W matrix. More precisely:

X = def [X + |X -] = W × t [H + |H -] + (4) X = X + -X -+ X 0 = W × t (H + -H -) + X 0 + (5)
. The notation [A|B] represents the r x 2c horizontal concatenation of two r x c matrices A and B. The first equation describes the NMF decomposition of X. The second one reconstructs the original matrix as the sum of baseline, high and low deviations.

The new approach proposed here relies on the more general k(= 3)-dimensional algorithmic version of NMF, known as Non-negative Tensor Factorization (NTF, [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF]), to approximate the original matrix by a three-dimensional n × f × 2 tensor X:

X = def [X + |X -] = W ⊗ H ⊗ Q + (6)
Note that all elements of X are non-negative. Similarly to the NMF-based algorithm described above, this decomposition provides a canonical clustering based on the meta-features provided by H matrix. Let us emphasize that this H matrix is common to the two slices of the tensor X. In other words, what we obtain is not the same as PosNegNMF, which performs two NMF decompositions on X + and X -respectively. The H matrix describes common features having an influence on either higher or lower deviations of the scores from their median.

PosNegNMF yields an approximation X + -X -= W × t H with H = H + -H -. Thus, each column of H may have signed entries, which complicates interpretation since the features under study are non-negative scores. In contrast, our NTF-based approximation yields a non-negative H matrix in which each column describes common features that have an impact on either higher or lower deviations of the scores from their median, as determined by the sign of the difference Q[1, column] -Q[2, column]. The non-negativity of the H-loadings makes interpretation as straightforward as it is for the standard NMF. For this reason, we propose to call our approach Separative NMF, or S 2 NMF, since it aims at separating patterns of features with either positive or negative deviations.

Precision and sparsity

We will now evaluate the comparative advantage brought by S 2 NMF, by varying the level of control of parsimony of the variables. More precisely, we observe how the approximation error and the composition of the meta-features evolve as the model is constrained to work with the most significant features. In order to focus on the sole effect of the parsimony, we start from a fully populated data matrix. Missing data are replaced for each feature, by the average of non-missing values of this feature. Doing this prevents from any unwanted numerical side-effect that can appear in the presence of a sparsity constraint applied to a matrix with empty entries.

To perform this experiment, we proceed as follows:

• The number of components is set to 6.

• We run in parallel an ordinary NMF model with its S 2 NMF counterpart.

• Both models are first launched with a sparsity target of 0.9 on the matrix H of meta-features. To achieve such a sparse representation, each vector of H is approximated after each update of the numerical resolution by a close vector (in the sense of quadratic norm) with the required degree of sparsity [START_REF]Block coordinate descent for sparse nmf[END_REF], [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF]). By doing so, some features show zero coefficients on all meta-features in matrix H. This is the sparse version of the model.

• The features showing only null loadings at the previous step are dropped from the matrix X, and the model is run again on the reduced matrix with no sparsity constraint: this is the restricted version of the model.

Figure 7 summarizes the steps of this process. Table 2 presents the results of this experiment. The columns of the tables are defined as follows:

• Dimension is the number c of components, here set at 6 for all models.

• Features refers to the number of features effectively included in the factorization.

• H-sparsity is a tuning parameter of the NMF/S 2 NMF algorithm, representing a target level of sparsity, or a proportion of null values, in matrix H.

• The Error column is the relative squared error of the approximation.

We see that the ordinary NMF model collapses when the sparsity is constrained at the high level of 90%. The mean squared error of the factorization becomes greater than the variance of the original data. In order to recover an acceptable error level, we have to drop 41 features out of 77 that do not contribute to the composition of any meta-feature, leaving us with 36 significant features. By contrast, starting directly with a Separative NMF approach in sparse mode, effectively reduces the error to an acceptable level of 9%. Forty five features are flagged as irrelevant, leaving the dense model with a remaining 32 columns.

Clusterization efficiency

The advantage of S 2 NMF over NMF is more easily seen by looking at the ordered heatmaps presentations of the matrix. Figures 8 and9 show NMF-and S 2 NMF-based heatmaps of the restricted matrix. The restriction takes out the features that have been eliminated in the first step with a target sparsity set at 90 %.

It is visually clear that S 2 NMF has a greater ability than NMF to discern and separate clusters, at both ends of the score range (blue or red aggregates).

If we try to formalise this visual intuition, the mathematical translation is that the distribution of scores within clusters is more homogeneous than in the totality of the data. This homogeneity can be quantified thanks to the notion of entropy.

Let us recall its definition: for any probability distribution p on a finite set:

X = {x 1 , . . . , x n }
the entropy of the distribution is given by

E(p) = i=n i=1 p(x i ) × log 1 p(x i ) (7)
Information theory tells us that log(1/p(x)) is (up to an additive constant) the optimal number of bits required to encode the value x, given its probability of occurrence p(x). From this viewpoint, E(p) appears to be the expected encoding length of a random pick from X. The minimum possible entropy is 0, attained for a distribution concentrated at 100 % on a single value. The maximum entropy is realized for a uniform distribution on X, and its value is then log(n). Let us use a normalized entropy defined by

E * (p) = E(p) log(n) (8)
The efficiency of a NMF/S 2 NMF-based clustering can be assessed by looking into the diagonal blocks of the sorted matrix. We expect to observe less dispersion, hence less entropy, among the values a) taken by the features belonging to a cluster i, and b) restricted to those observations also belonging to cluster i. The rationale here is to group together observations having high scores on a given meta-feature (read in matrix W ), hence having likely high scores on the features having a high loading in this meta-feature (read in matrix H). The advantage of S 2 NMF vs NMF is in the ability to group together not only the high values, but also the low values which can contain important information as well.

"Clustering entropy" shown in Table 3 represents the relative entropy of the clustering induced by each model. It is calculated as follows:

• The normalized entropy of each diagonal cluster (i, i) is calculated. The i-th diagonal cluster is a submatrix restricted to features in cluster, i and to observations having their dominant score on feature i.

• The clustering entropy is the average of relative entropies measured on all diagonal clusters, weighted by their sizes.

• This average entropy on diagonal clusters can be compared with the normalized entropy of the whole matrix. The difference defines the "Entropy delta".

Entropy delta Restricted NMF -0.77% Restricted S 2 NMF -8.03% It is not a surprise to observe lower average entropies in the clusters, compared to the respective original matrices. This is an indirect proof of the efficiency of general NMF methods applied to this particular dataset. Moreover, we see that the loss in entropy is higher (-8.03% vs -0.77%) with S 2 NMF compared to NMF, and this is a numerical confirmation of the visual difference between the respective heatmaps.

Interpreting meta-features

Let us select as final model a 6-components S 2 NMF, after elimination of the features from an initial run with sparsity 0.9. By taking a closer look at the definition of the meta-features, some coherence appears in the most influential features, i.e. those having the highest (positive) coefficients. 3 show the 6 most influential features of each meta-feature. Interestingly, the topics of these influential features are interrelated. We can therefore assign a main theme to each meta-feature. For example, the theme related to feature V1 is mostly about Human Rights (hrts) and Human Resources (hr), while V2 is clearly concentrated on Environment (env). S 2 NMF has traded dimension for interpretability, by drastically lowering the dimension (from 77 to 6) at the price of getting more fuzzy, but still providing meaningful features.

Rank

Chart 10 shows the correspondence, or incidence matrix, between the 6 clusters identified by S 2 NMF and the 11 industrial sectors of the stocks. Each cell shows the proportion of the stocks in a given sector (in rows) classified in a given cluster (in columns). It is visually clear that the two clusterings are not mutually redundant. S 2 NMF provides a data-driven, endogenous way of grouping stocks with respect to their ESG scores. For example, we see that the Utilities sector has about half of its members classified on the cluster 1, whose dominant theme is Human Rights (table 3). Similarly, the cluster 2 (Environment) is only dominant for the Health Care sector, and poorly represented in Utilities and Materials.

Conclusion and future work

We have illustrated here the application of NMF to extra-financial score data. In contrast to principal component analysis, which focuses on risk reduction by creating mutually orthogonal factors, NMF allows the information contained in the data to be summarized faithfully using meta-variables, where the original variables appear as linear combinations with positive coefficients. These meta-variables are not constrained by mutual orthogonality; they induce a natural clustering of the original variables and observations, at the cost of much lower error than would be produced by a PCA of the same dimension. Moreover, the S 2 NMF variant presented here allows, by symmetrizing the role of low and high scores, to further improve the results in terms of error and cluster relevance.

An induced result is the interpretability of the meta-features thus created, as shown by the comparison of their definitions with the prior knowledge of the meaning of the variables, which was not taken into account by the decomposition.

Finally, we have shown that an initial application of NMF with a high level of sparsity allows for the removal of variables that are insufficiently relevant.

Future directions of research are:

• The determination of heuristics independent of any prior knowledge of the domain, to optimize the number of components retained in the factorization. Possible heuristics are combinations of sparsity, robustness, precision and interpretability metrics.

• Generalizing the use of S 2 NMF as a pre-processing of variables, in a features engineering step upstream of a supervised learning process. The aim here is twofold: to reduce computation time and, above all, to reduce the risk of overlearning. Replacing a large number of partially redundant features by a much smaller number of interpretable features containing globally the same amount of information is a good policy in machine learning, regardless of the domain.

• The clustering induced by S 2 NMF, being independent from sectors, can induce a smart strategy for the imputation of missing data.

• S 2 NMF is applied to 2-dimensional data, and involves the use of 3-dimensional tensors to perform the factorization. A seemingly interesting generalization is to transpose respectively PosNegNMF and S 2 NMF to 3-dimensional data (with the addition of a time axis). The same way S 2 NMF uses an additional dimension to factorize the duplicated matrix (as described in equation 3), the 3-dimensional factorization would be achieved with 4-dimensional tensors, and yield PosNegNTF and SNTF.

We thank Vigeo-Eiris for making some datasets available to the authors. We also thank Vincent Margot for useful discussions and remarks.

Figure 1 Figure 2 :

 12 Figure 1 Figure 2: Each cell represents the average fill rate of a given feature within a given sector.

Figure 3 :

 3 Figure 3: The average fill rate by feature.

Figure 4 :Figure 5 :

 45 Figure 4: Cumulative explained variance of the first n components or eigenvectors.

Figure 6 :

 6 Figure 6: Coefficients of the 6 meta-features sorted by decreasing order.

Figure 7 :

 7 Figure 7: A schematic representation of the various factorizations applied to the initial data. Rectangles with rounded corners represent input or output data, while rectangles with slanted edges represent algorithmic treatments.

Figure 8 :Figure 9 :

 89 Figure 8: The original matrix reordered according to NMF clusters, and restricted to the 37 features preserved by the sparsity condition.

Figure 10 :

 10 Figure 10: Incidence matrix between industrial sectors and S 2 NMF clusters.

Table 1

 1 illustrates the fill rate (number of shares for which the variable is filled in, divided by the number of shares). The fill rate is calculated as an average per industry sector.

		Average fill rate
	Consumer Staples	88.41 %
	Consumer Discretionary	85.52 %
	Materials	84.36 %
	Health Care	83.89 %
	Utilities	83.54 %
	Industrials	82.37 %
	Energy	79.22 %
	Real Estate	79.02 %
	Communication Services	78.03 %
	Information Technology	77.18 %
	Financials	76.21 %

Table 1 :

 1

Sectors sorted by decreasing average fill rate, the average being taken over all features. The global average fill rate is 81.1 %

Table 2 :

 2 Comparison of NMF and S 2 NMF in sparse and dense versions.

Table 3 :

 3 Relative entropies for NMF and S 2 NMF restricted models.

Table 3 :

 3 Assigning themes to meta-features based on their dominant features

Table


R Ravaud

G Lemarquand
email: guy.lemarquand@univ-lemans.fr

MODELLING AN IRONLESS LOUDSPEAKER BY USING THREE-DIMENSIONAL ANALYTICAL APPROACHES

This paper presents some improved three-dimensional expressions of the magnetic field created by tile permanent magnets uniformly and radially magnetized for the design of ironless loudspeaker structures. All the expressions determined have been reduced to compact forms.

We use these expressions for the optimization of ironless loudspeaker structures in which the radial field must be radially uniform. Indeed, as ring permanent magnets radially magnetized are rather difficult to manufacture, these magnets are replaced by assemblies of tile permanent magnets radially and uniformly magnetized. We present an example of ironless loudspeaker structure that has been optimized with our threedimensional approaches.

INTRODUCTION

Tile permanent magnets can be used for manufacturing ring permanent magnets radially magnetized by stacking them together. Such a way of manufacturing a ring is commonly realized because ring permanent magnets radially magnetized are difficult to fabricate. Consequently, this alternative method consists in using tile permanent magnets uniformly magnetized. However, they do not exhibit the same properties as tiles radially magnetized and can lower the quality of the radial field created by the permanent magnet structure.

Many authors have studied the field created by ring [START_REF] Durand | Magnetostatique[END_REF]- [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF], cylinder [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]- [START_REF] Selvaggi | Computation of the three-dimensional magnetic field from solid permanent-magnet bipolar cylinders by employing toroidal harmonics[END_REF], tile permanent magnets [START_REF] Furlani | Field analysis and optimization of ndfeb axial field permanent magnet motors[END_REF]- [START_REF] Ravaud | Magnetic field produced by a tile permanent magnet whose polarization is both uniform and tangential[END_REF] or by disk conductors [START_REF] Azzerboni | Magnetic field evaluation for disk conductors[END_REF]- [START_REF] Conway | Noncoaxial inductance calculations without the vector potential for axisymmetric coils and planar coils[END_REF]. All these analytical methods have enabled authors to study the stiffness or the force between permanent magnets for magnetic bearings [START_REF] Delamare | Classification and synthesis of permanent magnet bearing configurations[END_REF] [START_REF] Moser | Optimization of repulsive passive magnetic bearings[END_REF] or couplings [START_REF] Yonnet | Analytical calculation of permanent magnet couplings[END_REF] [START_REF] Elies | Analytical study of radial stability of permanent magnet synchronous couplings[END_REF]. More generally, optimizations can be carried out with analytical or semi-analytical approaches [START_REF] Perigo | General expressions for the magnetic flux density produced by axially magnetized toroidal permanent magnets[END_REF]- [START_REF] Zhilichev | Calculation of magnetic field of tubular permanent magnet assemblies in cylindrical bipolar coordinates[END_REF]. An application of the analytical calculation of the magnetic field produced by permanent magnets is the design of ironless loudspeakers [START_REF] Berkouk | Analytical calculation of ironless loudspeaker motors[END_REF]- [START_REF] Ravaud | Ironless loudspeakers with ferrofluid seals[END_REF].

This paper presents improved three-dimensional expressions of the magnetic field created by tile permanent magnets uniformly magnetized. The analytical expressions we present are not only more compact but also more simple to use than the ones determined in [START_REF] Ravaud | Magnetic field produced by a tile permanent magnet whose polarization is both uniform and tangential[END_REF]. Indeed, the contribution of the magnetic pole surface densities located on the straight faces of the tile permanent magnet can be described in a simpler expression than the one determined in [START_REF] Ravaud | Magnetic field produced by a tile permanent magnet whose polarization is both uniform and tangential[END_REF]. Even though the two problems are different between this paper and [START_REF] Ravaud | Magnetic field produced by a tile permanent magnet whose polarization is both uniform and tangential[END_REF], the form of the analytical expressions that must be determined analytically is the same. Then, these expressions are used in order to design an ironless loudspeaker structure. The first section presents the way of obtaining these analytical expressions by using the Coulombian model. Then, this paper discusses the wave rate of the radial field created by a stacked structure of tile permanent magnets uniformly magnetized. A prototype has been built to show the suitability of such an analytical three-dimensional approach. All the expressions obtained in this paper are available online [26].

MODELLING TILE PERMANENT MAGNETS WITH THE COULOMBIAN MODEL

This section presents the geometry studied and the model used for modelling the magnetic field created by tile permanent magnets uniformly magnetized. We explain how to find the surface densities located on the faces of a tile permanent magnet uniformly magnetized. The magnetic polarization of the tile permanent magnet studied is J = 1 T for the rest of this paper.

Modelling the tile permanent magnet using the Coulombian approach

The geometry considered and the related parameters are shown in Fig 1 . The coulombian model asserts that a permanent magnet whose polarization is J can be represented by using fictitious magnetic charges which are located inside the magnet and on its faces. Thus, there are both a magnetic pole volume density σ * v inside the magnet and a magnetic pole surface density σ * s located on its faces. These magnetic pole densities can be determined by using the following expressions: The magnetic polarization vector J is expressed as follows:

σ * v = -∇. J (1)
J = -J x u x -J y u y (3)
where J x = J cos( θ 2 +θ 1 2 ) and J y = J sin( θ 2 +θ 1 2 ). By using (1) with (3), we deduct that:

σ * v = 0 (4) This is a well-known result: a permanent magnet whose magnetic polarisation is uniform has not any magnetic pole volume densities.

By applying (2) to n 1 , n 2 , n 3 and n 4 (see Fig 2), we find:

+ J 4πµ 0 r 2 r 1 z 2 z 1 sin( θ 2 -θ 1 2 ) u(θ 2 ) u(θ 2 ) 3 drdz + J 4πµ 0 r 2 r 1 z 2 z 1 sin( θ 2 -θ 1 2 ) u(θ 1 ) u(θ 1 ) 3 drdz (6)
where µ 0 is the permeability of the vacuum, J is the magnetic polarization of the tile permanent magnet and u is expressed as follows:

u(θ i ) = (r -r j cos(θ -θ i ) u x -r j sin(θ -θ i ) u y + (z -z) u z (7)
where j = r 1 when dS = r 1 d θdz (the inner surface) or j = r 2 when dS = r 2 d θdz (the outer surface). The integration of ( 6) leads to the magnetic field components along the three defined axes: H r (r, θ, z), H θ (r, θ, z) and H z (r, θ, z) (expressed in A/m) which are given by ( 11), ( 14) and [START_REF] Moser | Optimization of repulsive passive magnetic bearings[END_REF]. We use the relations

x k = cos(θ -θ k ), x 2k = cos(2(θ -θ k )) y k = sin(θ k -θ) and t k = cot(θ -θ k ).
In addition, we define ξ i,j , η i,k and α i,j,k as follows:

ξ i,j = r 2 + r 2 i + (z -z j ) 2 -2rr i cos(θ -θ) (8)
η i,k = 2r 2 i -4rr i x k + r 2 (1 + x 2k ) (9)
α i,j,k = r 2 + r 2 i -2rr i x k + (z -z j ) 2 (10)
3.2. Radial Component H r (r, θ, z)
The radial component H r (r, θ, z) of the magnetic field created by one tile permanent magnet uniformly magnetized can be expressed as follows:

H r (r, θ, z) = J 4πµ 0 2 i=1 2 j=1 (-1) (i+j) ℘ (1,r) (i, j) + J 4πµ 0 2 i=1 2 j=1 2 k=1 (-1) (i+j) ℘ (2,r) (i, j, k) (11)
with

℘ (1,r) (i, j) = θ 2 θ 1 r i (z j -z)(-r + r i cos(θ -θ)) cos( θ 1 +θ 2 -2 θ 2 ) (r 2 + r 2 i -2rr i cos(θ -θ))ξ i,j d θ (12) ℘ (2,r) (i, j, k) = - r(rx k (x 2k -1) + 2r i y 2 k ) -r 2 (-1 + x 2k )η i,k arctan (z -z j )η i,k -r 2 (x 2k -1)α i,j,k -x k log [z -z j + α i,j,k ] (13) 3.3. Azimuthal Component H θ (r, θ, z)
The azimuthal component H θ (r, θ, z) of the magnetic field created by one tile permanent magnet uniformly magnetized can be expressed as follows:

H θ (r, θ, z) = J 4πµ 0 2 i=1 2 j=1 (-1) (i+j) ℘ (1,θ) (i, j) + J 4πµ 0 2 i=1 2 j=1 2 k=1 (-1) (i+j) ℘ (2,θ) (i, j, k) (14) with ℘ (1,θ) (i, j) = θ 2 θ 1 r 2 i (z -z j ) sin(θ -θ) cos θ 1 +θ 2 -2 θ 2 ξ i,j r 2 + r 2 i -2rr i cos(θ -θ) d θ (15)
and

℘ (2,θ) (i, j, k) = z -z j α i,j,k y k sin θ 2 -θ 1 2 - r i y k t k α i,j,k sin θ 2 -θ 1 2 arctan z -z j ry k + ry k x k t k α i,j,k sin θ 2 -θ 1 2 arctan z -z j ry k (16) 3.4. Axial Component H z (r, θ, z)
The axial component H z (r, θ, z) of the magnetic field created by one tile permanent magnet uniformly magnetized can be expressed as follows:

H z (r, θ, z) = J 4πµ 0 2 i=1 2 j=1 (-1) (i+j) ℘ (1,z) (i, j) z u r u u z 1 2 2 1 J = -J cos ( ) -J cos ( ) J = -Ju r u u a a 2 1 A B a = 2 x y u u x y u Figure 3.
Representation of two tile permanent magnets: the left one is radially magnetized, the right one is uniformly magnetized

+ J 4πµ 0 2 i=1 2 j=1 2 k=1 (-1) (i+j) ℘ (2,z) (i, j, k) (17)
with

℘ (1,z) (i, j) = θ 2 θ 1 -r i ξ i,j cos θ 1 + θ 2 -2 θ 2 d θ (18)
and

℘ (2,z) (i, j, k) = sin θ 2 -θ 1 2 log [r i -r cos(θ -θ k ) + ξ i,j ] (19)
3.5. Comparison of the magnetic field created by a tile permanent radially magnetized and a tile permanent magnet whose polarization is both uniform and radial Tile permanent magnets whose polarization is both uniform and radial are generally used instead of tile permanent magnets radially magnetized (Fig 3). However, the magnetic field they produce exhibits some differences and thus can lower the quality of the permanent magnet devices.

The three components of the magnetic field created by each tile permanent magnet (radially or uniformly magnetized) are represented in Fig. 4. The stacking of several tile permanent magnets uniformly magnetized can be used for creating the same radial field as the one created by a ring permanent magnet radially magnetized. Strictly speaking, the radial field undulates in front of the tile magnets, especially at the place where two tiles are stacked. We propose here to determine the wave rate of the radial field created by the ring of stacked tile permanent magnets versus the number of tiles used. The wave rate W R is defined by

W R = 1 2H rmoy (|Sup[H r ]| -|Inf [H r ]|) ( 20
)
where H rmoy is the mean radial field at a given radial distance and a given axial distance from the ring stacked with tile permanent magnets uniformly magnetized.

We study three configurations owing 8, 16 and 32 tile permanent magnets uniformly magnetized. Each tile has the same dimensions.

We take for all the illustrative examples the following dimensions : the inner radius r 1 is 0.025 m, the outer radius r 2 is 0.028 m, the height z 2z 1 is 0.003 m and J = 1 T, the observation radial distance is r = 0.024 m and the observation axial distance is z = 0.001 m. We represent in Fig 7 the radial field for the three configurations. Figure 7 clearly shows that the higher the number of tile permanent magnets used is, the lower the wave rate of the radial field is. We can extend this remark when we use N tile permanent magnets uniformly magnetized. Figure 8 represents the wave rate of the radial Representation of the wave rate of the radial field created by a ring of stacked tile permanent magnets uniformly magnetized versus the number of tiles used. For each tile, the inner radius is 0.025 m, the outer radius is 0.028 m, the height is 0.003 m field created by a ring of stacked tile permanent magnets versus the number of tiles used. Such a representation is useful because it gives interesting information about how many tile permanent magnets should be used if the wave rate must be inferior to a given percentage.

In addition, when N → ∞, we find the same radial field value as the one determined by Babic and Akyel [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF].

ILLUSTRATION: DESIGN OF AN IRONLESS LOUDSPEAKER STRUCTURE

This section presents an example of utilization of the expressions determined in this paper. To do so, let us consider a typical ironless structure made of tile permanent magnets stacked together for creating a radial field in the air gap as shown in Fig. 9. In this optimized configuration, this ironless loudspeaker is composed of three ring permanent magnets made of 16 tile permanent magnets uniformly and radially magnetized.

As stated previously, the radial field undulates versus θ. By using the analytical expressions determined in this paper, we can optimize the tile permanent magnet dimensions for obtaining the magnetic field whose undulation is the less important. However, the number of tile permanent magnets used must be as small as possible. Indeed, such ironless structures are built for being commercialised. As it is both difficult and expensive to stack a large number of tile permanent magnets together, we must look for an optimal configuration that creates a radial field whose wave undulation is as small as possible though the Representation of an ironless loudspeaker with three ring permanent magnets and two ferrofluid seals number of tile permanent magnets used is as small as possible too. The height of the tile permanent magnets is fixed in our configuration (h = 3 mm) and the radial width of each tile permanent magnet equals r 2 -r 1 = 3 mm. Therefore, the only parameter that can be otpimized is the angular width of each tile permanent magnet, that is, the number of tile permanent magnets used in the ironless loudspeaker structure. This optimization has been performed by using the threedimensional analytical expressions determined in this paper and the expression of the wave undulation defined in [START_REF] Perigo | General expressions for the magnetic flux density produced by axially magnetized toroidal permanent magnets[END_REF]. The wave undulation must be inferior to at least 10%. The radial field created H is used for two different purposes. First, it is used for moving the emissive piston shown in Fig 9 by using the Laplace's force defined as follows:

F L = I(t)d l ∧ (µ 0 H) (21)
where I(t) is the current in the coil that is located on the emissive piston and l is the length of the coil wire. The second objective of this radial field is to contribute of the creation of a ferrofluid seal. Indeed, the magnetic field created by such a structure can be used for creating two ferrofluid seals whose shapes influence directly the

p m (r, z) = µ 0 M H r (r, θ, z) 2 + H θ (r, θ, z) 2 + H z (r, θ, z) 2 (22)
Therefore, the final shape of the ferrofluid seal does not undulate as the radial field as shown in Fig 7 . However, the radial field is the most important component that contributes to the deformation of the seal shape. We present in Fig 10 an example of prototypes that has been built and optimized in our laboratory.

As shown in Fig 10, between two tile permanent magnets, we see the undulation of the magnetic field that stems mainly from the radial field undulation. This behaviour is predicted theoretically with the analytical expressions determined in section 2.

CONCLUSION

This paper has presented improved three-dimensional analytical expressions of the components of the magnetic field created by tile permanent magnets uniformly and radially magnetized. First, we have presented the expressions of the three components H r (r, θ, z), H θ (r, θ, s) and H z (r, θ, z) created by a tile uniformly magnetized in a compact form. Such expressions are important for practical purposes because tile permanent magnets uniformly magnets can be used to manufacture a entire ring radially magnetized. Indeed, they are easier to fabricate than arc-shaped permanent magnets radially magnetized. In addition, we have seen that the wave rate of the radial field created by a stacked magnet structure depends directly on the number of tile permanent magnets used. The more we use tile permanent magnets uniformly magnetized, the less the wave rate is important. This element of information can be very useful for manufacturers involved in the design of ironless loudspeakers. We have built an ironless loudspeaker prototype with an optimized number of tile permanent magnets. The effects observed experimentally are well predicted theoretically. As the expressions determined in this paper have a very low computational cost, they can be easily used for optimizing ironless structures.

Figure 1 . 1 Figure 2 .

 112 Figure 1. Representation of the configuration considered. The ring inner radius is r 1 , the ring outer one is r 2 , its height is z 2z 1 , its angular width is θ 2θ 1

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Representation of the the radial, axial and azimuthal components of the magnetic field created by a tile which is either radially magnetized (thin line) or uniformly magnetized (thick line). r in = 0.025 m, r out = 0.028 m, z bz a = 0.003 m, z = 0.0015 m, r = 0.024 m and θ 2θ 1 = π 3 rad

Figure 7 .

 7 Figure 7. Representation of the radial field created by a ring of stacked tile permanent magnets uniformly magnetized. The inner radius is 0.025 m, the outer radius is 0.028 m, the height is 0.003 m and J = 1 T. (Thick line : θ 2 -θ 1 = π 32 rad, thin line : θ 2 -θ 1 = π 16 rad, dashed line : θ 2θ 1 = π 8 rad)

Figure 8 .

 8 Figure 8. Representation of the wave rate of the radial field created by a ring of stacked tile permanent magnets uniformly magnetized versus the number of tiles used. For each tile, the inner radius is 0.025 m, the outer radius is 0.028 m, the height is 0.003 m

Figure 9 .

 9 Figure 9. Representation of an ironless loudspeaker with three ring permanent magnets and two ferrofluid seals

Figure 10 .

 10 Figure 10. Photography of the ironless loudspeaker structure: we have used 3 times 16 tile permanent magnets for creating two ferrofluid seals, therefore we have θ 2θ 1 = π 8 rad; r 1 = 0.02485 m, r 2 = 0.04485 m, z 2z 1 = 0.02 m