X Song

L Leotoing
email: lionel.leotoing@insa-rennes.fr

D Guines

E Ragneau

Effect of continuous strain path changes on forming limit strains at necking and fracture of DP600 sheets

Keywords: In-plane biaxial tensile test, FLCN, FLCF, non-linear strain path

The strain path may change in actual sheet metal forming processes, so the determination of formability of sheet metal should consider the non-linear strain path. For identifying the forming limit (FL) strains under non-linear strain path, a conventional two-step procedure with unloading is classically used to produce the strain path change, which results in no continuous measure of strain. The in-plane biaxial tensile test with a cruciform specimen is an interesting alternative to overcome the drawbacks of conventional method. The strain path change can be made without unloading during a single test. In this work, the experimental FL strains of DP600 sheets under two types of non-linear strain path are investigated and then compared with those under linear strain paths. The Oyane ductile fracture criterion is used in the finite element simulation to predict the experimental results.

Introduction

Most of the forming limit curves at necking (FLCNs) and the forming limit curves at fracture (FLCFs) have been experimentally determined using methods that produce proportional loading with insignificant changes in strain paths. However, in some innovative forming processes, such as single point incremental forming, complex strain paths with several changes have been identified by different authors. As an example, the figure 1 shows the evolution of the strain path during an incremental forming operation of a truncated aluminium cone with a wall angle of 60 o [START_REF] Eyckens | Strain evolution in the single point incremental forming process: digital image correlation measurement and finite element prediction[END_REF]. As can be seen, drastic changes in the strain path are observed during the forming operation, without unloading, and a very high level of deformation can be achieved, well beyond the conventional forming limit strains at necking. In addition, it has been demonstrated that strain path change has an influence on the FLCN and more specifically on the high level of forming limits (FL) [START_REF] Da Rocha | Prediction of the forming limit diagrams of anisotropic sheets in linear and non-linear loading[END_REF].

Because the strain path change obviously affects the FL strains at necking, no unified FLCN exists in the strain space to represent the formability of sheet metal [START_REF] Yoshida | Path-dependence of the forming limit stresses in a sheet metal[END_REF]. The strain paths in sheet metal forming are classically divided into two types: linear strain path and non-linear strain path. Traditionally, the identification of strain-based FLCN or FLCF is limited to the sheet metal undergoing linear strain paths. Three types of experimental FLCN (Figure 2) are concluded by Barata et al. [START_REF] Da Rocha | Prediction of the forming limit diagrams of anisotropic sheets in linear and non-linear loading[END_REF] to investigate the effect of strain path change on the FL strains at necking. For type 1, the FLCN is determined under proportional loading. Every point of the FLCN is defined by major and minor strains at necking. As shown in Figure 2 (a), the whole FLCN (ABCD) is produced by different strain paths from biaxial stretching (OA) to uniaxial tension (OD) through plane-strain tension (OC). The red dash line shows the FLCF, which is higher than the FLCN [START_REF] Song | Investigation of the forming limit strains at fracture of AA5086 sheets using an in-plane biaxial tensile test[END_REF]. For type 2, the FLCN is identified under two steps of proportional loading. Different prestrain levels under a constant strain ratio P1 are used, and then an abrupt change is produced towards the strain ratio P2 for every prestrain level. As shown in Figure 2 (b), the curve ABCD is identified by the uniaxial prestrain followed by biaxial stretching, and the curve AEFD is produced by different levels of biaxial prestrain followed by uniaxial tension. For type 3, the FLCN is also determined under two steps of proportional loading. The strain path and strain level are kept constant in the first step, and different strain paths are used in the second step. As shown in Figure 2 (c), the biaxial stretching is followed by various ratios of linear strain path for the curve ABCD, while the uniaxial tension is followed by different ratios of linear strain path for the curve EFGH.

Figure 1: Evolution of the major and minor surface strains for a 60 o wall angle cone during incremental forming process [START_REF] Eyckens | Strain evolution in the single point incremental forming process: digital image correlation measurement and finite element prediction[END_REF] The strain path change is always produced by a two-step procedure in conventional tests to identify the experimental FLCN under non-linear strain paths. In the first step, oversized specimen is produced by the tensile test, Marciniak or bulge test to obtain prestrains in uniaxial tension, plane-strain tension or biaxial stretching. In the second step, the oversized pre-strained sheet is deformed under new loading condition in order to achieve a strain path change. Ishigaki [START_REF] Ishigaki | Deformation analysis of large sized panels in the press shop[END_REF] at Toyota Motors Company have applied the multi-stage loading to increase the formability of sheet metal. The engineers firstly used four stages of prestrain in uniaxial tension to reach 37% of major strain. Whereafter, two stages of biaxial stretching were adopted which dramatically improved the formability. Graf and Hosford [START_REF] Hosford | The influence of strain-path changes on forming limit diagrams of Al 6111 T4[END_REF] have investigated the influence of strain path change on the FL at necking of aluminum alloy 6111-T4. Different prestrain levels in various loading conditions are considered for the FLCNs. The results show that the abrupt change of strain path can significantly affect the FL at necking. The prestrain in biaxial stretching decreases the level of FL strain if succeeded by uniaxial tension or planestrain tension. The prestrain in uniaxial tension increases the level of FL strain if succeeded by biaxial stretching and plane-strain tension, when the direction of principal strain is maintained but decreases if the direction of principal strains is changed after the prestrain. The prestrain in plane-strain tension leads to a small increase of the level of FLCN when the principal strain direction is maintained, but a substantial decrease is happened with the change of principal strain direction. Stoughton et al. [START_REF] Stoughton | Path independent forming limits in strain and stress spaces[END_REF] have analyzed the FLCNs for 2008-T4 sheet metal under different prestrains reported by Graf and Hosford [START_REF] Graf | Effect of changing strain paths on forming limit diagrams of Al 2008-T4[END_REF] in early years. Recently, Volk et al. [START_REF] Volk | Failure prediction for nonlinear strain paths in sheet metal forming[END_REF] have adopted the two-step procedure to plot experimental FLCN of a high strength steel HC300X under non-linear strain paths and six different prestrains from uniaxial to equibiaxial are considered. In addition, the experimental results of two-step procedure have been widely used for comparison with numerical results in many predictive researches [START_REF] Chow | Prediction of forming limit diagrams for AL6111-T4 under non-proportional loading[END_REF][START_REF] Uppaluri | An analytical approach for the prediction of forming limit curves subjected to combined strain paths[END_REF][START_REF] Cao | Prediction of localized thinning in sheet metal using a general anisotropic yield criterion[END_REF]. To represent the formability of metallic sheets under non-linear strain paths, the stress-based forming limit concept was proposed in the early 1980s [START_REF] Arrieux | Determination of an intrinsic forming limit stress diagram for isotropic sheets[END_REF]. It is found that the forming limit stress curve (FLSC) is almost path-independent [START_REF] Yoshida | Path-dependence of the forming limit stresses in a sheet metal[END_REF]. Many authors [START_REF] Stoughton | A general forming limit criterion for sheet metal forming[END_REF][START_REF] Stoughton | Sheet metal formability analysis for anisotropic materials under nonproportional loading[END_REF][START_REF] Butuc | An experimental and theoretical analysis on the application of stress-based forming limit criterion[END_REF] have investigated the FLSC but the main disadvantage of this approach is that the stress state cannot be measured experimentally. The FLSC is calculated from the measured forming limit strains using postulated constitutive behaviours, i.e. yield function and hardening rule. Some disadvantages can be found in the conventional two-step procedure for controlling the strain path. Firstly, the conventional method requires different experimental methods to obtain different strain paths. Secondly, the unloading between two steps is obligatory and no continuous measure of strain is produced, while non-linear strain paths are observed without any unloading in some forming processes (like single point incremental forming process [START_REF] Eyckens | Strain evolution in the single point incremental forming process: digital image correlation measurement and finite element prediction[END_REF]). Thirdly, multiple strain path changes are impossible. Only simple prestrains can be realized which makes difficult to study the formability with several changes of strain path. Using the in-plane biaxial tensile test with the cruciform specimen to control the strain path can be an interesting alternative to overcome the drawbacks of conventional methods. The potentiality of the in-plane biaxial tensile test with a cruciform specimen to study the effect of strain path change on the FLCN and to investigate the FLCF under linear strain paths was validated by the presenting authors [START_REF] Leotoing | Investigations of the effect of strain path changes on forming limit curves using an in-plane biaxial tensile test[END_REF][START_REF] Song | Investigation of the forming limit strains at fracture of AA5086 sheets using an in-plane biaxial tensile test[END_REF]. Firstly, the strain path during the test can be directly controlled by the motion of actuators along the two axes, which is sufficient to cover the whole forming limit diagram under linear and non-linear strain paths. Secondly, the unloading is not obligatory and the measurement of strain is continuous during strain path changes. Lastly, complex strain path with several changes can be envisaged for continuous process conditions (temperature, strain rate). For the in-plane biaxial tensile test, the critical point remains the design of the cruciform specimen. A high number of cruciform shapes have been proposed in literature, depending on the mechanical behavior to be studied: yield locus, hardening, forming limits ... Recently, a standardized shape, based on the works of Kuwabara et al. [START_REF] Kuwabara | Measurement and analysis of differential work hardening in coldrolled steel sheet under biaxial tension[END_REF] has been proposed. However, this shape is not suitable for forming limit investigations. Strains are very low in the central area and fracture takes place in the arms. A first geometry made from a metal sheet with an initial thickness of 4 mm was validated for the determination of FLCN and FLCF [START_REF] Song | Investigation of the forming limit strains at fracture of AA5086 sheets using an in-plane biaxial tensile test[END_REF] for AA5086. For linear strain paths, the correlation between results from conventional Marciniak tests and biaxial tensile tests is rather good for FLCN [START_REF] Leotoing | Cruciform shape benefits for experimental and numerical evaluation of sheet metal formability[END_REF]. Recently, a new shape with an initial thickness of 2 mm was optimized by the authors and validated for linear strain paths [START_REF] Song | Characterization of forming limits at fracture with an optimized cruciform specimen: Application to DP600 steel sheets[END_REF]. The initial thickness reduction from 4 mm to 2 mm permits to study metal sheets with a larger range of applications. For the two geometries, a thickness reduction is necessary in the central region and slits on arms reduce transverse rigidity and strain localization in the junction area of two perpendicular arms.

A predictive method can cover a wide range of linear and non-linear strain paths, which is easier and cheaper than experimental determination [START_REF] Chow | Prediction of forming limit diagrams for AL6111-T4 under non-proportional loading[END_REF]. In the literature, most of the studies for prediction of FL under non-linear strain path are investigated with the M-K model. Graf and Hosford [START_REF] Graf | Calculations of forming limit diagrams for changing strain paths[END_REF] have analyzed the effect of abrupt and gradual changes of strain path on FLCN of aluminum alloy 2008-T4 by using the M-K model. The FLCN moves to the right and down compared with the original FLCN when the first step corresponds to biaxial stretching, while it moves to the left and up when the uniaxial tension is adopted in the first step. Yoshida et al. [START_REF] Yoshida | Path-dependence of the forming limit stresses in a sheet metal[END_REF] have analyzed the influence of strain path change on the FL stresses of sheet metals by using the M-K model. Two types of loading are considered: one type has unloading between two steps of loading while the other one has no unloading. The results show the strain path has influence on the FL in strain space, while the FL in stress space depend on the strain path only for the second type of loading. Kuroda and Tvergaard [START_REF] Kuroda | Effect of strain path change on limits to ductility of anisotropic metal sheets[END_REF] have investigated the effect of strain path change on the FLCN by using the M-K model. The predictive results show that FLCN depends on the unloading between the two steps of loading. The elastic straining is very important for the dependence. When the unloading is adopted between the two steps, the stress point of yield surface moves in a different way, which leads to a different flow localization response. Chow et al. [START_REF] Chow | Prediction of forming limit diagrams for AL6111-T4 under non-proportional loading[END_REF] have extended an anisotropic damage model to predict the FLCNs of AA6111-T4 sheet under non-proportional loading. The prestrains in uniaxial tension, plane-strain tension and biaxial stretching are considered. The second loading direction is parallel or perpendicular to the first loading direction. The investigation of the damage model is made by comparing the calculated results of FLCN with the experimental results published by Graf and Hosford [START_REF] Hosford | The influence of strain-path changes on forming limit diagrams of Al 6111 T4[END_REF]. When the first and second loading direction are parallel, the correlation between the predictive FLCN and the experimental result is good. However, the predictive FLCN is lower than the experimental FLCN when the two loading directions are perpendicular.

Isik et al. [START_REF] Isik | Formability limits by fracture in sheet metal forming[END_REF] have analyzed the FLCF by a theoretical method. Based on the damage mechanics and the correlation between FLCF and fracture toughness, the authors have concluded that the FLCF is a straight line. On the contrary, the FLCN is not a material property and is strain path dependent. In addition, the linear profile of FLCF is validated experimentally under linear strain paths, while non-linear strain paths are not considered.

In this study, the in-plane biaxial tensile test is used to investigate the FL strains at necking and fracture of DP600 sheet metal under two types of non-linear strain paths. A ductile fracture criterion for predicting the onset of fracture under non-linear strain paths will also be discussed.

2 Experimental procedure and results

Cruciform specimen

As shown in Figure 3, an optimized cruciform specimen was proposed by the authors to identify the FL strains at fracture of DP600 sheets under linear strain paths [START_REF] Song | Characterization of forming limits at fracture with an optimized cruciform specimen: Application to DP600 steel sheets[END_REF]. This specimen is designed for sheets with an initial thickness of 2 mm and the shape is not specific to DP600 sheets since its performance was also validated for other materials, like aluminium alloys. For generating strain localization in the central area of the specimen (fracture zone), a circular arc profile is adopted in the thickness direction of the circular reduced zone and six slots with an optimized arrangement are used. The strain path in the central area of the specimen can be controlled directly by four independent actuators.

Figure 3: Geometry of the cruciform specimen [START_REF] Song | Characterization of forming limits at fracture with an optimized cruciform specimen: Application to DP600 steel sheets[END_REF] As shown in Table 1, different linear strain paths including equibiaxial stretching, plane-strain tension and uniaxial tension are produced by different velocity ratios of actuators acting on the arms of cruciform specimen.

Table 1: Different velocity ratios of actuators for producing linear strain paths

Velocity ratio R 1 0.05 (plane-strain tension) uniaxial tension Velocity along axis X (mm/s) 1 2 1

Velocity along axis Y (mm/s) 1 0.1 free As shown in Figure 4, two types of non-linear strain paths are considered in this work. For the type 2-A, two steps of loading without unloading are used, in which the first step corresponds to the uniaxial tension and the equibiaxial stretching is used for the second step. For the type 2-B, the first step corresponds to equibiaxial stretching and the plane-strain tension is adopted for the second step. As shown in Table 2, five levels of displacement in uniaxial tension for the type 2-A and six levels in equibiaxial stretching for the type 2-B are considered in the first step. One can see on figure 3 that the thickness reduction in the center is obtained by removing material from one side only. The opposite side of the specimen is plane. This plane surface is used to generate the speckle pattern. The non-symmetrical geometry of the specimen due to the thickness reduction of the central zone can lead to some out of plane specimen bending effects. Those effects have been numerically investigated and it is shown that the out of plane specimen displacement before the fracture is less than 1 mm (confirmed by experiments) which leads to negligible bending effects. A camera is perpendicularly positioned with respect to the plane surface which remains quasi-plane during the test. Consecutive images are recorded by a fast camera and an acquisition of 250 images/s is chosen. The true in-plane strain components are evaluated by the Digital Image Correlation (DIC) method. A 2D system is then used. Due to the large deformations, the reference image is updated at each step of displacement calculation by considering the previous image.

Forming limit strains at fracture under non-linear strain paths

The time-dependent method is relatively efficient for identifying the beginning of fracture. The FL strains at fracture of DP600 sheets under linear strain paths have been identified by the presenting authors [START_REF] Song | Characterization of forming limits at fracture with an optimized cruciform specimen: Application to DP600 steel sheets[END_REF]. This method is now applied for determining the beginning of fracture and the FL strains at fracture under the two types of non-linear strain path.

The identification of onset of fracture and FL strains at fracture for the pre-displacements of 2.0 mm in uniaxial tension and 4.0 mm in equibiaxial stretching are presented here as examples. Figure 5 shows the equivalent strain field of the central area at the time 0.004s before fracture with different pre-displacements. The evolution of the equivalent strain at the calculated area (2.4 × 2.4 mm size) for the two pre-displacements are plotted in Figure 6. It can be found that the macroscopic crack (Figure 7 and Figure 8) occurs with an abrupt change of equivalent strain in a very short time (0.004 s) under each type of non-linear strain paths. The equivalent strain is calculated by assuming volume consistency and is defined by :

ε = √ 2 3 ∑ i,j ε 2 ij (1)
The volume consistency is debatable after the onset of necking due to the growth of voids in the material but the evolution of equivalent strain is only used to detect the time just before rupture. With the identified time, the major and minor strains of the calculated area at the moment just before the occurrence of macroscopic fracture are used as the FL strains at fracture. For other pre-displacements in uniaxial tension and in equibiaxial stretching, the same method is adopted to identify the FL strains at fracture. Figure 9 shows the FL strains at fracture under different types of non-linear strain path. In this figure, the solid markers and dashed lines represent the FL strains at fracture and the strain paths, respectively. As shown in Figure 9 (a), different pre-displacements in uniaxial tension correspond to different levels of prestrain (see Table 3). After the pre-displacement, the equibiaxial stretching is used in the As shown in Figure 9 (b), the abrupt change of strain path from equibiaxial stretching to plane strain is obvious for each pre-displacement. Different levels of prestrain in equibiaxial stretching are shown in Table 4. The major strain at fracture decreases with increasing the prestrain from 8% to 38%. As can be seen in Figure 9 (a), the instability shown in Figure 5 (a) corresponds to a zero extension line. For all the tested paths, before rupture, the minor strain no longer increases and remains constant. Only the major strain grows and a near plane strain state is observed although an equibiaxial stretching is applied. Figure 10 shows the forming limit strains at fracture under linear and non-linear strain paths (Type 2-A and Type 2-B). As described in previous study [START_REF] Song | Characterization of forming limits at fracture with an optimized cruciform specimen: Application to DP600 steel sheets[END_REF], the forming limit strains at fracture under linear strain paths is fitted by a straight line with k=-0.39 and A=0.63 , while those under non-linear strain paths is fitted by a straight line with k=-0.49 and A=0.68. For linear strain paths, the experimental results are in accordance with the ones of Gruben et al. [START_REF] Gruben | Determination of ductile fracture parameters of a dual-phase steel by optical measurements[END_REF], especially for uniaxial and equibiaxial loadings. For plane strain condition, the forming limits are lower than the ones of Figure 10(a). In the work of Gruben et al., Nakajima and Marciniak tests were performed. Recently, forming limits at rupture were given for DP600 by Cheng et al. [START_REF] Cheng | A modified lou-huh model for characterization of ductile fracture of dp590 sheet[END_REF] by using Nakajima tests. In this last work, forming limits at fracture for equibiaxial loading are very close to the ones of Figure 10(a). For the other strain paths, forming limits in [START_REF] Cheng | A modified lou-huh model for characterization of ductile fracture of dp590 sheet[END_REF] are systematically lower. As shown in figure 10, the difference between the FLCFs for linear strain paths and non-linear strain paths is small. In other words, the strain path change has very small effect on the forming limit strains at fracture.

Forming limit strains at necking under non-linear strain paths

The critical ratio method was successfully used to identify the FL strains at necking of DP600 sheet metal under linear strain paths [START_REF] Song | Characterization of forming limits at fracture with an optimized cruciform specimen: Application to DP600 steel sheets[END_REF][START_REF] Zhang | A comparative study of different necking criteria for numerical and experimental prediction of flcs[END_REF]. The same method is adopted to determine the FL strains at necking under the two types of non-linear strain paths. In order to apply the critical ratio method, the evolutions of equivalent strain measured respectively in zone A and B (Figure 5) are plotted (Figure 11). In this figure, the onset of necking with the pre-displacement of 2.0 mm in uniaxial tension and 4.0 mm in equibiaxial stretching is very clear. When a localized necking appears in the central zone, the equivalent strain in zone A increases abruptly due to the plastic instability, while the equivalent strain in zone B stops any increase and shows a saturation. The moment of the saturation of equivalent strain in zone B is used to define the onset of localized necking. The major and minor strains in zone A associated to this moment define the FL strains at necking. To identify precisely this moment, the ratio of equivalent strain increment between zone A and B is calculated and compared to a critical value. The increment is calculated between each image. The critical value of the equivalent strain increment ratio depends on the position of zone B. For the selected zone B (Figure 5), a value of 8 permits to detect any saturation and then the onset of the plastic instability. Figure 12 depicts the comparison of FL strains at necking under linear and non-linear strain paths. For the type 2-A of non-linear strain path, the FLCN shifts up with respect to the original FLCN under linear strain paths, and a small increase of formability at necking is observed in Figure 12 (a). This phenomenon is also reported by Leotoing et al. [START_REF] Leotoing | Investigations of the effect of strain path changes on forming limit curves using an in-plane biaxial tensile test[END_REF] for an aluminium alloy sheet and Kuroda et al. [START_REF] Kuroda | Effect of strain path change on limits to ductility of anisotropic metal sheets[END_REF] for a cold-rolled lowcarbon steel, when the level of prestrain in uniaxial tension is not high. In their researches, a high level of prestrain in uniaxial tension can reduce the formability at necking, which results in an abrupt decrease for the FLCN. In other words, the level of prestrain in uniaxial tension can either improve or reduce the formability at necking. However, as shown in Figure 12 (a), no reduction of formability at necking is observed under a high level of prestrain (major strain beyond 30%) in uniaxial tension. Therefore, it can be concluded that for the DP600 sheet metal, the prestrain in uniaxial tension improves the formability at necking under the type 2-A of non-linear strain path. As shown in Figure 12 (b), the FLCN under type 2-B shifts to the down comparing with the original FLCN under linear strain paths. The same results are obtained by Graf and Hosford [START_REF] Hosford | The influence of strain-path changes on forming limit diagrams of Al 6111 T4[END_REF]. For small prestrain (8% or 10%) in equibiaxial stretching, there is no obvious difference between the FLCNs under linear strain paths and non-linear strain paths.

3 Numerical simulation

Numerical model

As shown in Figure 13, a predictive model based on the optimized cruciform geometry has been proposed by the finite element (FE) method. The same displacement ratios than the experimental ones are applied along the two perpendicular directions to follow different strain paths. Displacements are applied at the ends of the two arms of the specimen at a distance of 80 mm far from the central point. For linear strain paths, the same displacement is imposed on each arm. For non-linear strain paths, the displacement values are given in table 2 for the step 1. The FE simulations are produced by the commercial software ABAQUS (implicit solver). Because the cruciform specimen is symmetrical, one-quarter is adopted to save the calculation time. Linear tetrahedral solid elements are used for the mesh. A refined mesh with the size of 0.25 mm is defined in the central area. The identification of DP600 sheet behaviour under in-plane equibiaxial loadings has been performed in previous work [START_REF] Liu | Identification of strain rate-dependent mechanical behaviour of DP600 under in-plane biaxial loadings[END_REF] for a strain rate ranging from quasi-static to intermediate strain rate. Based on this work, the material constants and material models used in the present work are briefly recalled hereafter. For the isotropic elasticity, Young's modulus of 200 GPa and Poisson's ratio of 0.3 are considered. For the plastic behaviour, the associated normal flow rule is assumed and Hill48 yield criterion for plane stress condition is used. 5, the parameters of Hill48 yield criterion for DP600 have been calculated from three anisotropic coefficients proposed by Ozturk et al. [START_REF] Ozturk | Effects of anisotropic yield functions on prediction of forming limit diagrams of DP600 advanced high strength steel[END_REF].

The general strain rate dependent hardening law on the basis of Ludwick formulation [START_REF] Liu | Identification of strain rate-dependent mechanical behaviour of DP600 under in-plane biaxial loadings[END_REF] has been used: The parameters σ 0 , K , n , m 1 and m 2 are given in Table 6.

σ = σ 0 εm1 + Kε p n εm2 (2)
Table 6: Identified parameters of the rate-dependent hardening law for Hill48 yield criterion [START_REF] Liu | Identification of strain rate-dependent mechanical behaviour of DP600 under in-plane biaxial loadings[END_REF] σ 0 (MPa) K (MPa) n m 1 m 2 339.2 839.7 0.3864 0.0052 0.0158

Ductile fracture criterion

The Oyane criterion [START_REF] Oyane | Criteria for ductile fracture and their applications[END_REF] has been used to predict the FLCF of DP600 sheets under linear strain paths by the presenting authors [START_REF] Song | Characterization of forming limits at fracture with an optimized cruciform specimen: Application to DP600 steel sheets[END_REF]. It should be noted that this criterion is based on an uncoupled ductile model which does not account for the progressive damage of the material. If damage mechanics certainly affect the onset of necking and consequently the fracture strains, uncoupled criteria are easier to implement and to calibrate than microscopic failure criterion. In this study, the Oyane criterion will be used to predict the FLCF under non-linear strain paths. The hydrostatic pressure σ h and the equivalent stress σ are considered in this criterion. Two parameters A and C need to be identified.

∫ εf 0 (1 + σ h Aσ ) dε p = C (3)

Numerical results

Figure 14 shows the numerical FLCFs under different strain paths calculated with the Oyane ductile fracture criterion. The parameters A and C were identified by minimizing the gap between the experimental fitting line and the fracture forming line for linear strain paths (A=-1 and C=0.45) (Figure 14(a)). The same values of A and C are used to predict the numerical FLCFs under non-linear strain paths (Figure 14(b)). The numerical FLCFs under linear and non-linear strain paths (type 2-A and 2-B) are almost overlapped. This conclusion is in accordance with the experimental results since the same fracture line can be defined for both linear and non-linear strain paths. Besides, it is possible to calibrate the Oyane criterion parameters in order to predict accurately the experimental results under linear and non-linear strain paths. The criterion is validated under a certain range of stress triaxiality σ h σ (from 1/3 in uniaxial tension to 2/3 in equibiaxial tension). To discuss the validity of the criterion for other strain states, we think that low stress triaxiality values could be obtained with the in-plane biaxial test by applying a tension along one direction and a compression along the other one.

Conclusions

An optimized cruciform specimen has been used with the in-plane biaxial tensile test to identify the FL strains at fracture and necking of DP600 sheets under non-linear strain paths. Two types of non-linear strain paths without unloading are used in the tests. Some conclusions are obtained: (1) The FL strains at fracture of DP600 sheets under non-linear strain paths can be fitted by a straight line. The difference between the experimental FLCFs under linear and nonlinear strain paths is very small, which indicates that the strain path change has very small effect on the experimental FLCF.

(2) The prestrain in uniaxial tension followed by equibiaxial stretching can increase the level of FL strain at necking of DP600 sheets, while the prestrain in equibiaxial stretching followed by plane-strain tension can decrease it. (3) The numerical FLCFs calculated by the Oyane criterion predict well the experimental FL strains at fracture under linear and non-linear strain paths. The strain path change has almost no effect on the prediction of FLCF. Calibration of the Oyane criterion can be performed from simple tests, under proportional loadings.

Figure 2 :

 2 Figure 2: Three types of experimental FLCN

Figure 4 :

 4 Figure4: Two types of non-linear strain path[START_REF] Song | Identification of forming limits at fracture of dp600 sheet metal under linear and unloaded non-linear strain paths[END_REF]

Figure 5 :

 5 Figure 5: Equivalent strain field of the central area at the time 0.004s before the fracture with different predisplacements

Figure 6 :

 6 Figure 6: Identification of the onset of fracture under different types of non-linear strain paths

Figure 7 :

 7 Figure 7: Specimen surfaces for DP600 sheets with the pre-displacement of 2.0 mm in uniaxial tension

Figure 8 :Figure 9 :

 89 Figure 8: Specimen surfaces for DP600 sheets with the pre-displacement of 4.0 mm in equibiaxial stretching

Figure 10 :

 10 Figure 10: Comparison of forming limit strains at fracture under different strain paths

Figure 11 :

 11 Figure 11: Identification of necking for DP600 sheets under different types of non-linear strain path

Figure 12 :

 12 Figure 12: Comparison of forming limit strains at necking different strain paths

Figure 13 :

 13 Figure 13: 3D Mesh of the cruciform specimen

  (a) Linear strain paths (b) Non-linear strain paths

Figure 14 :

 14 Figure 14: Prediction of forming limits at fracture DP600

Table 3 :

 3

	Major strain	7%	9%	13%	17%	27%

Major strains at fracture under different pre-displacements in uniaxial tension Pre-displacement 2.0 mm 2.5 mm 3.0 mm 3.5mm 4.0 mm

Table 4 :

 4 Major strains at fracture under different pre-displacements in equibiaxial stretching

	Pre-displacement 2.0 mm 2.5 mm 3.0 mm 3.5mm 4.0 mm 4.5 mm
	Major strain	8%	10%	15%	22%	30%	38%

Table 5 :

 5 Lankford's coefficient and Hill48 yield parameters

	r 0	r 45	r 90	F	G	H	L	M	N
	0.89 0.85 0.420 0.529 0.471 1.500 1.500 1.282


C Krier
email: krier@dice.ucl.ac.be

F Rossi

D François
email: francois@inma.ucl.ac.be

M Verleysen
email: verleysen@dice.ucl.ac.be

A DATA-DRIVEN FUNCTIONAL PROJECTION APPROACH FOR THE SELECTION OF FEATURE RANGES IN SPECTRA WITH ICA OR CLUSTER ANALYSIS

Keywords: features extraction, functional projection, Independent Component Analysis, clustering

Prediction problems from spectra are largely encountered in chemometry. In addition to accurate predictions, it is often needed to extract information about which wavelengths in the spectra contribute in an effective way to the quality of the prediction. This implies to select wavelengths (or wavelength intervals), a problem associated to variable selection. In this paper, it is shown how this problem may be tackled in the specific case of smooth (for example infrared) spectra. The functional character of the spectra (their smoothness) is taken into account through a functional variable projection procedure. Contrarily to standard approaches, the projection is performed on a basis that is driven by the spectra themselves, in order to best fit their characteristics. The methodology is illustrated by two examples of functional projection, using Independent Component Analysis and functional variable clustering, respectively. The performances on two standard infrared spectra benchmarks are illustrated.

Introduction

Predicting a dependent variable from the measure of a spectrum is a frequently encountered problem in chemometrics. For example, one might try to predict the concentration in sugar, nitrogen or alcohol or the percentage of fat content, from infrared spectra measured on juice, grass, wine or meat samples, respectively. Accurate predictions may help avoiding heavy and costly chemical measurements, as spectra measures are generally much more affordable and can be automated easily.

However, accurate predictions are not sufficient in many applications. Indeed it is often required to have an insight on which chemical components are responsible for the prediction, in order to interpret the model and further develop models and products. As basic chemical components are associated to limited spectra ranges, finding which parts (ranges) of the spectra are responsible for the prediction is important in many applications too.

For this reason, prediction models that use only parts of the spectra (after variable selection) will be preferred, when interpretation is sought. Not all prediction models are able to force the fact that only parts of the spectrum are used. For example, the traditional PLS uses loadings that are combinations of all spectral variables. Even if some of them have a larger weight in the model than others, in most situations it can hardly be concluded that only some parts of the spectra participate to the prediction model.

Selecting ranges in spectra may be seen as a variable selection problem where the variables are the wavelengths in the spectrum. The variable selection problem at hand is specific in the sense that its dimensionality is high. The dimension is indeed the number N of wavelengths available in spectra, which can easily reach several hundreds or thousands in modern, high-resolution spectrometers. Obviously, with such high resolution, consecutive wavelengths are highly correlated. This leads to the well-known colinearity problem, and all other difficulties related to the high dimension and the dependencies between coordinates (complex models, difficulties of convergence in the latter, overfitting, curse of dimensionality, etc.).

In addition, selecting individual variables is not a prefect solution though: because of the colinearity between consecutive variables, there is often no more reason to select an individual wavelength than the preceding or next one. The selection of groups of consecutive variables (or wavelength ranges) seems thus less arbitrary and can in addition reduce the computational burden of the selection process. Selecting groups of consecutive wavelengths may be achieved through the so-called functional modeling approach: contrarily to vector approaches (which are by far the most common ones in spectrometric modeling), the order (or indexes) of variables is taken into account. In other words, the information contained in the fact that an infrared spectrum is a smooth function of the wavelengths is not lost. Remind that in vector approaches, such as PLS, the numbering of variables is irrelevant: any swapping or permutation will rigorously lead to the same prediction results; the smoothness of spectra, although an important information, is not exploited.

The principle of functional approaches is to project the spectra on a basis of smooth functions, and then to use the coordinates of each spectrum on the basis as new variables [1]. In this way, the number of variables is reduced in the projection, and the functional nature of spectra is exploited through the smoothness of the basis functions. For example, B-splines are often used as basis functions [2], [3], as they offer several advantages: the projection is easy to compute, the basis functions cover approximately equal ranges, and taking into account the derivative of spectra is easily achieved [4]. Projection on wavelets [5], [6] also preserves the functional nature of spectra, but do not have the equal range property.

In both cases however, the basis is chosen a priori, without taking the spectra into account. In other words, the same basis is used for any problem (an exception is the adjustment of some parameters, such as the number of B-splines). Obviously, this solution might not be optimal: "better" bases could be found if the spectra themselves were taken into account. For example, it may happen that some parts of the spectra are highly informative, necessitating an increased number of basis functions in that part, while others are flat and/or not correlated to the dependent variable, thus necessitating less or no basis functions. It is thus preferable to choose the basis, or the groups of consecutive wavelengths, in a spectra-driven way, for increased prediction performances and interpretation.

Selecting appropriate groups of wavelengths according to the specificities of the spectra is the idea developed in [START_REF]Conference on Artificial Neural Networks[END_REF]. In this paper, an information-theoretic criterion is used to group together (to cluster) variables that carry similar information content. As the clustering is driven by the value of the spectra, the size of the groups is automatically adjusted to optimize the prediction performances. However, the approach in [START_REF]Conference on Artificial Neural Networks[END_REF] is not functional, in the sense that the indexes of the spectral variables are not used. The procedure may result in the selection of consecutive variables, as a consequence of their colinearity, but it may not be the case too.

In this paper, we develop a framework for the functional analysis of smooth spectra, through their projection on bases that are spectra-driven. In that sense, we combine the advantages of functional analysis (taking a priori into account the smooth nature of spectra) and of data-driven variable grouping (using a basis adapted to the data, for increased prediction performances and interpretation). Two examples of this data-driven functional approach are detailed. The first one decomposes the spectra into independent components that will be further used as basis functions. Independent Component Analysis (ICA) [START_REF]Wiley Series on Adaptative Learning Systems for Signal Processing[END_REF] is used for the decomposition. Because of the independence between basis functions, it will be shown that they occupy a limited range of the spectral variables, providing interpretability. The second one borrows the variable clustering idea from [START_REF]Conference on Artificial Neural Networks[END_REF], replacing a group of variables by a "mean" characteristic, though exploiting the functional nature of the data.

The two proposed approaches are intended to select ranges of spectrometric variables, and not only isolated variables. They differ in the fact that the clustering methodology is a crisp method, in the sense that each spectral variable is assigned to a single cluster, while this is not exactly the case for the ICA approach, where each spectral variable has a certain "weight" in each cluster (but close to zero for most of them). It is expected that selecting ranges of spectral variables will help interpretability, in particular to identify basic components that are contained in the samples measured by spectrometry. This does not mean however that the goal of the method is to extract the spectra themselves of the basic component; this would be a much more ambitious goal, and would require further independence hypotheses between components that are not found in practice. This paper is organized as follows. Section 2 reviews recent methods to select variables in spectrometric modeling. The proposed methodology is then presented in Section 3 at a general level. Details relating to the two projection basis envisaged in this paper are given in the next section. Section 5 explains how to select the projected variables in order to build a prediction model. Finally, Section 6 presents and discusses the experimental results on two Infra-Red spectroscopy datasets

Theory: state-of-the-art

This section presents state-of-the-art methods to select variables in spectrometric modeling. The first part is about functional approaches with a priori fixed bases, and the second part addresses a variable clustering, but non-functional, solution.

Functional projection

In general, the Infra-Red (IR) spectra can be seen as smooth functions sampled at different wavelengths corresponding to the spectral variables. This allows considering functional analysis as a tool to model and interpret such data. B-splines are commonly used as a basis to project the spectra into a lower dimensional space [2]. These functions present indeed a "local" behavior, in the sense that each spline corresponds to a specific range of spectral variables; this specificity preserves the interpretability of the results. Moreover, the groups of spectral variables defined by the splines contain only consecutives spectral variables, which present often a high degree of colinearity. Splines, and thus ranges of spectral variables, can then be selected through their coefficients. This selection is often based on the mutual information between the coefficients and the parameter to predict [START_REF] Rossi | [END_REF].

The choice of this kind of functional basis is however made a priori, without taking the information included in the spectral variables into account, which is certainly not optimal. This means also that different problems will be handled with the same functional basis. In addition, all ranges of wavelengths defined by the splines share the same size, and are not adapted to the shape (peaks and flats) of the spectra.

Clustering

Another way of grouping similar spectral variables is the clustering. Van Dijk and Van Hulle propose in [START_REF]Conference on Artificial Neural Networks[END_REF] to use a hierarchical clustering with a similarity notion adapted to spectral variables and clusters of variables. This similarity is based on the mutual information between the spectral variables.

Contrarily to the functional projection, clusters can have different sizes, adapted to the particularities of the spectra. The information included in the spectral variables is also exploited when determining the clustering in the similarity measure. This clustering is however not a functional projection approach, because spectral variables grouped in a cluster are not necessarily consecutive in the spectra; two variables corresponding to different wavelengths can indeed be grouped in a cluster, while another variable related to an intermediate wavelength is included in a second one. This phenomenon, avoided by the projection on Bsplines, may be problematic from the interpretability point of view.

Methodology

In the present paper, the prediction problem is solved in a three-step procedure. First, spectra are replaced by their projection on a functional basis. Two types of bases are considered here: a basis formed by the independent components (computed by ICA) of the spectra, and clusters of spectral variables. A classical preprocessing such as derivation, Standard Normal Variate (SNV) [10] or Multiplicative Scatter Correction (MSC) [11] , can be applied before the projection.

In the first case the new variables are the coefficients of the projection on the independent components.

In the second case the new variables are mean values taken over the clusters. The projection of the spectra is unsupervised for both the ICA and clustering approaches (only the number of clusters in this second case is chosen in a supervised manner), in the sense that the projection is not designed explicitly to optimize the quality of the prediction of the dependent variable.

Second, the new variables are selected by a supervised criterion. Indeed the number of variables resulting from the first step, which makes no or little use of the dependent variable, will generally be too large for an efficient use in the prediction model. Therefore, the second step consists in selecting which of the resulting variables will be used in the prediction model.

The proposed feature selection methods are linear in the sense that the final selected features are specific linear combination of the original features, both in the case of the ICA and the case of variable clustering. A linear model built on those variables is also a linear model build on the original variables and is thus likely to have lower performances than a PLSR model which is optimal among such linear models.

Consequently, a nonlinear model has to be considered in order to preserve the quality of the prediction.

As a nonlinear prediction model will be used, the criterion used to select the variables must be nonlinear too. The Mutual Information (MI) will be used, as a nonparametric nonlinear criterion to evaluate the usefulness of variables, in a greedy procedure (forward search followed by an exhaustive search).

Third, the prediction model itself will be built on the variables selected after step two. The nonlinear model used in the experimental section is a traditional Least-Squares Support Vector Machines (LS-SVM) [12]. Other models such as Multi-Layer Perceptrons (MLP), Support Vector Machines (SVM) and Radial-Basis Function Network (RBFN) [START_REF] Moody | Connectionist Models Summer School (CA) 1988[END_REF] could be used too.

Sections 4 and 5 will detail the functional projection and variables selection steps, respectively. Section 6 will present some experimental results.

Functional projection

This section shows two possible ways to reduce the number of variables in spectra. The functional character is preserved in both cases. The first solution consists in projecting the spectra on a basis found by Independent Component Analysis. The second solution consists in clustering the variables under a functional constraint. The first one is that the sources are statistically independent; this leads to the so-called Independent Component Analysis (ICA). It is also necessary to make some hypotheses on the mixture. Though other options are possible, the hypothesis that will be made here is that the mixing process is linear and Numerous methods based on some measure of the independence between the estimations S ˆare able to find an appropriate matrix W. Note that the ICA solution suffers from two indeterminacies, a possible permutation and a scaling of the sources. The resulting WA product will therefore not necessarily be the identity matrix; it will be a matrix with one non-zero element per line and per column.

In our context of infrared spectra analysis, the observations z i are the spectra of the analyzed products, and the sources s i are supposed to be the spectra of their independent basic chemical components. The implicit reasonable hypothesis of the above model is thus that the spectrum of a mixing of basic components is a linear combination (with coefficients corresponding to the proportions) of the component spectra.

In the above model each column of A corresponds to a single source. The column is thus considered as the set of coefficients (for each analyzed spectrum) corresponding to this source, in other words the coefficients of the projection on this source. The second step mentioned in Section 3 and detailed in Section 5 will therefore consist in selecting columns of A.

A last question remains to be discussed in the context of the ICA projection: the choice of the number of independent components. Because of the intrinsic formulation of the ICA problem, this number cannot exceed the number of observed spectra. However, a lower number is preferred, both to facilitate the subsequent variable selection, and to avoid numerical problems in the ICA solution. A simple solution to the choice of the number of independent components is to reconstruct the spectra (from Z = A S) with a limited number of sources, and to choose the lower number which leads to a small reconstruction error (the threshold may be fixed according to the level of noise in the data).

Clustering

The spectral variable clustering method described in this section is based on the combination of the variable clustering algorithm proposed in [START_REF]Conference on Artificial Neural Networks[END_REF] with functional constraints that facilitate the interpretation of clusters of variables.

When a smooth function is discretized over a fine grid, the resulting values are generally very similar on sub-intervals of the original range with "reasonable size": smoothness implies that f(w 1 ) and f(w 2 ) are very close when w 1 and w 2 are. As a consequence, spectral variables corresponding to close wavelengths take nearly identical values. It is therefore tempting to reduce the resolution of the spectra by downsampling them over wavelength ranges in which they are almost constant (on each individual spectrum). This can be done by clustering identical or similar variables together, and then by replacing each cluster of variables by the average of the variables it contains. This procedure is similar to a standard k-means algorithm, applied here to variables (features) instead of observations. The prototype summarizing a cluster is therefore a vector whose number of elements is the number of spectra: it corresponds to a new variable.

However, the goal of this procedure is to reduce the number of features while maintaining interpretation possibilities. It is therefore important to design the clustering procedure such that it produces clusters that fulfill the latter requirement. This is done first by requesting clusters to contain only consecutive variables: each cluster is associated this way to a wavelength range that is much easier to interpret than a set of arbitrary and possibly unrelated spectral variables. It should be noted that this requirement is not artificial: while a spectral variable is naturally more similar to a spectral variable with a close wavelength, the very smooth nature of some spectra leads to "long distance" similarities, as illustrated on Figure 1. This Figure represents the absolute value of the correlation between a spectral variable and all the other variables for the Tecator dataset (see Section 6 for details on this dataset). While the first three maximal correlations correspond to the variable itself and its two neighbors (leading to the wavelength range [906.6, 910.6]), the fourth maximally correlated variable corresponds to the wavelength 942.9 nm.

Clustering variables without functional constraint could thus lead to disconnected hardly interpretable clusters. The second element used to enforce easy interpretation is the use of a simple similarity criterion between variables, namely the absolute correlation, i.e.

(

) ( ) ( ) ( ) ( ) This measure has several advantages over the mutual information used in [START_REF]Conference on Artificial Neural Networks[END_REF]. It is both easier and faster to evaluate via standard estimators. However, the most important aspect is that the correlation measures only simple (linear) dependencies between variables, whereas the mutual information evaluates complex nonlinear relationships. When it is used to build a predictive model, this is a valuable characteristic, but in the case of variables clustering, it can lead to clusters that are quite difficult to interpret. For instance, if 2 , the mutual information between the variables will be high; they will tend to be considered as similar and hence put in the same cluster. However, the actual values of X 2 and X 1 will be quite different, a confusing fact for the experts. Moreover, the possible complex nonlinear relationship between variables in a cluster turns the design of a synthesis (for example mean) variable into a difficult problem. If variables in a cluster are highly correlated, they are also highly correlated to the mean variable of the cluster, while if they have high mutual information, they have no particular reason to have a high mutual information with the mean variable. Indeed, a high mutual information between two variables corresponds to the existence of an approximate bijective function mapping one variable to the other: even if there are exact one-to-one mappings between several variables in a cluster, there might be no such mapping between those variables and their mean.

j i j i j i j i X Var X Var X E X E X X E X X - = Σ ) , (
X 2 =(X 1 )
In practice, the proposed method is based on an agglomerative (bottom-up) full linkage hierarchical clustering algorithm. It proceeds as a standard hierarchical clustering method, with the additional constraints of merging only adjacent clusters. For N original spectral variables, the algorithm works as follows:

1. Create N clusters, c 1 ={X 1 }, ..., c N ={X N }.
2. Compute the similarity between each cluster and its following cluster, i.e., between c i and c i+1 , as the absolute correlation between the variables they contain (i.e., X i and X i+1 ).

Repeat N-1 times:

a. Find the pair of consecutive clusters that are the most similar.

b. Merge the clusters.

c. Update the similarities between consecutive clusters as follows: the similarity between c i and c i+1 is the minimum of the absolute correlation between any pair of variables, the first one being in c i and the second in c i+1 .

At each step of the algorithm, the number of clusters is reduced by one, starting from N and going down to 1. Clusters have a hierarchical structure, summarized by a dendrogram. For each number of clusters between 1 and N, the original variables are replaced by one variable per cluster, defined as the mean of the variables contained in the cluster. As clusters correspond to sub-interval of the original wavelength range, the resulting new variables can be considered as forming a piecewise constant approximation of the original spectra. In this sense, the proposed methodology is related to the B-spline based functional representation used in [START_REF] Rossi | [END_REF] and recalled in Section 2. The major difference is that [START_REF] Rossi | [END_REF] used B-splines of order 4 or 5, with regularly spaced knots. The piecewise constant approximation build here corresponds to B-splines of order 1, with non-regular spacing between knots.

The algorithm described above constructs iteratively N successive clusterings (N-2 non-trivial ones). In order to choose one of them, the expected performances of the complete model (variable clustering + variable selection + nonlinear prediction) would be an ideal selection criterion. In theory this is achievable.

In practice however, a robust evaluation of each nonlinear model would require computationally-demanding 15 techniques such as cross-validation, leading to unaffordable computation times when N is large. An approach that does not make use of the full model is thus preferred. A linear measure (the correlation) between features being used to build the hierarchical clustering, it is suggested to use a linear model to evaluate each of the N clustering results. More precisely, for each M between 1 and N, the cross-validated performances of a linear regression model built on the clustering result with M classes are calculated. The optimal M is chosen as the one that maximizes this performance indicator.

Feature selection

New variables are extracted from the spectra after a projection or clustering procedure as detailed in the previous section. However, in the projection case the output to predict (the dependent variable) is not taken into account. In the clustering approach, it is not fully taken into account either: the similarity criterion used to group the original features together makes no use of the output to predict, the latter being used only to choose among the N clustering results. However, in this case, only a linear predictive model is used.

Therefore, in both cases, it is expected that the result of the projection or clustering step will contain a prohibitive number of variables (although reduced compared to the number of original variables in the spectra) to be entered directly in a nonlinear prediction model. Indeed the projection or clustering step groups together features that are sufficiently correlated, but does not provide a way to eliminate those (or the resulting clusters) that carry no information to predict the output. A further feature selection may thus be considered, at least if the selection criterion used for this step makes use both of the output variable, and of the nonlinear potential of the prediction model; the goal here is to select only the variables that are relevant for the prediction.

A selection procedure consists in two basic blocs: the relevance criterion and the selection procedure.

The solution described in the following sections is similar to the ones used in our previous works [START_REF] Rossi | [END_REF], [START_REF] Rossi | [END_REF].

Relevance criterion

First, a relevance criterion has to be chosen: how to measure if a variable, or a group of variables, is pertinent for the prediction (in the present context, a group of variables is not anymore a group of spectral variables, but rather a group of new variables obtained after functional processing). The best choice for the criterion is certainly to estimate the performances of the model itself. However, let us remind that we are in a context where the number M of variables is high. Testing all 2 M -1 models will obviously lead to a prohibitive computational cost. In addition with a greedy (forward) procedure as detailed below, the model performance criterion is not adequate. Indeed the forward procedure will begin with a single variable, then two, etc., and there is no reason that the model itself will perform well on a too small number of variables.

Instead of an approach based on prediction performances, a simpler one is preferred, where the criterion is evaluated independently from the model. The correlation could have been used, but it suffers from two drawbacks: it measures only linear relations, and is most usually limited to individual variables, not groups of. Therefore, the Mutual Information (MI) is preferred. The MI between X and Y is defined as

( ) ( ) ( ) ( ) ( ) ∫ = dxdy y x y x log y x , MI Y X Y X Y X μ μ μ μ , , , , Y X (1)
where μ X , μ Y and μ XY define the probability density function of variables X, Y and the joint X,Y variable respectively. Note that X and Y in the above definition may be multidimensional, which answers to the above mentioned requirement to handle groups of variables. In our context, X will denote one or a group of variables obtained after projection or clustering, and Y is the dependent variable.

The mutual information must of course be estimated on the learning set. The evaluation of the mutual information by equation ( 1) requires the knowledge of the densities. In practice, their estimation on the learning set is unreliable when the number of variables in X increases, due to the curse of dimensionality.

A better estimator has been proposed in [15] to handle multi-dimensional variables. This estimator is based on the search for nearest neighbors among the learning set in the spaces spanned by X and Y, therefore avoiding the direct estimation of the densities; it will be used in the experimental section of this paper.

Forward selection

Once the criterion is chosen, a procedure to select variables has to be defined. Remind that the variables considered here are those resulting from the first step, i.e. the projection or feature clustering. Each of them represents a cluster of initial variables that are consecutive in the spectra. The further selection described in this section does not force to select clusters that are consecutive. The idea is that several ranges of spectral variables may be of interest for the prediction problem (see Section 6 for examples); it is thus important to let the method select non-consecutive clusters.

As explained above, an exhaustive search on all combinations of variables among the set of available ones would be computationally prohibitive. A forward selection is preferred. It consists in iteratively selecting variables at each step of the procedure, without questioning a choice previously made. The first selected variable is the one which maximizes the MI with the dependent variable Y. The second one is chosen as to maximize the MI between Y and the group of two variables formed by the previously selected one and the new one to choose. The procedure is repeated by adding one variable at each step.

In practice, the estimated value of the MI might decrease at some step of the forward selection. This is however only a consequence of the bias and variance of the estimator. In theory indeed, the true value of the MI can only increase when variables are added in X. Therefore the maximum of MI cannot be used as stopping criterion (although such choice in sometimes found in scientific publications). Rather, the forward procedure is iterated until a "reasonable" number of variables is selected. By "reasonable", we mean a number P for which 1) the curse of dimensionality remains limited in the MI estimator, and 2) the following exhaustive search remains computationally affordable. In practice, the value of P may for example be limited to 7 or 8.

Finally, an exhaustive search on all remaining variables is performed. In Section 5.1 it was argued that building 2 M -1 non-linear models would have been prohibitive when M is large. Now than M is replaced by P, where P equals for example 7 or 8, building 2 P -1 models becomes affordable, even when nonlinear models are considered. Therefore 2 P -1 nonlinear models are built, their performance evaluated (for example by cross-validation), and the one which leads to a minimum of the estimated error is selected. In the next Section about experiments, LS-SVM models are used for this purpose. The details of the learning procedure for these models are given in [4].

In practice, the forward selection proceeds as follows:

1. Define V 1 as the X i that maximizes MI(X i ; Y).

For k from 2 to P,

define V k as the X i that maximizes MI(V 1 ,…,V k-1 , X i );Y).
3. For each subset L of {V 1 ,…,V P }, build a LS-SVM on L, evaluate the best performances of the LS-SVM by cross-validation.

4. Select the best subset as the one that maximizes the cross-validation performances.

Experiments

This Section introduces first the datasets from the food industry used in the paper to evaluate the proposed methods. The next part details the experimental methodology. The results are then presented and discussed.

Datasets

Two datasets are used to illustrate and evaluate the functional approaches described in this paper. The first database is Tecator [START_REF]Tecator meat sample dataset[END_REF]. It consists of 215 near-infrared spectra of meat samples, recorded on a Tecator Infratec Food and Feed Analyzer in the 850-1050 nm wavelength range. Each spectrum is discretized into 100 spectral variables. The spectral variables are the absorbance, defined by log(1/T), where T is the measured transmittance. All spectra have been normalized according to the SNV method (mean equal to zero and variance equal to 1). The spectra are used to predict the fat content of the meat samples and are divided into a learning and a test set. The learning set defined for the experiments contains 172 spectra, which leaves 43 spectra for the test.

The second dataset, Wine [START_REF] Prof | [END_REF], is related to the alcohol concentration in wine samples. In this case, spectra consist of absorbance measures recorded in the mid-infrared range at 256 different wavenumbers between 4000 and 400 cm -1 ; the number of spectral variables is therefore 256. The learning and test sets contain respectively 94 and 30 spectra. Three spectra (numbers 34, 35 and 84) are however considered as outliers and consequently removed from the learning set.

Experimental methodology

The experimental methodology is summarized in Figure 2. For each database, the functional approaches described in Sections 4 and 5 are compared to two more "classical" approaches in similar experimental conditions: the fist one is a non-linear model (LS-SVM) and the second one is a standard PLSR. These experiments are conducted on the original data to provide additional reference performances.

All experiments are carried out on identical training and test sets. In order to choose the metaparameters for the LS-SVM (i.e. the regularization parameter and standard deviations of the kernel), a 3fold cross-validation procedure is carried out. Again, the same subsets of spectra are used for all compared methodologies. These different approaches are implemented in Matlab® version 7.1, using the LS-SVM Toolbox (LS-SVMLab1.5) [START_REF] Pelckmans | LS-SVMlab : a Matlab/C toolbox for Least Squares Support Vector Machines[END_REF]. In the case PLSR, processing has been done with the R statistical software [START_REF]R: A Language and Environment for Statistical Computing[END_REF], using the pls package [START_REF] Mevik | [END_REF]. In the case of the two proposed approaches, the spectra are projected on the functional basis for the two datasets. The projected variables (respectively the independent components and the representatives of the clusters) are then selected by mutual information and a LS-SVM is built on the selected projected variables.

The ICA is performed by the FastICA version 2.5 toolbox [START_REF]FastICA for Matlab 7.x and 6.x[END_REF]. The options chosen for the analysis are:

"symmetric" (estimate all the independent components in parallel) and "stabilization" (stabilized version of FatsICA). The non-linear function used for the experiments is g(u) = u 3 . The clustering procedure leads to 27 clusters for Tecator and 30 for Wine. Concerning the selection of the projected variables by MI, the forward procedure is stopped at 7 variables, which means 2 7 -1 LS-SVM models built on all possible combinations of variables. As suggested in [15], the mutual information estimator takes 6 neighbors into account; this neighborhood is defined according to the Euclidian norm.

The performances of the predictive models are evaluated on a test set Ω, by the Normalized Mean Squared Error (NMSE). The NMSE is defined as follows:

( )

, ) ( 1 2 Y Var Y Y Q NMSE i i ∑ Ω Ω Ω - =
where i Y is the parameter to predict, i Y ˆ the corresponding prediction and Q Ω the number of data in Ω. The NMSE is thus always non-negative; a NMSE is greater than 1 means that the prediction performances are worst than if all predicted values where equal to the mean of Y [START_REF] Caelen | IEEE International Symposium on Computer-Based Medical[END_REF] (actually, the NMSE is equal to

1 Q Q Ω Ω - )
in this extreme case.

The most time-consuming step is the forward selection carried out after the clustering or the projection of the spectra on the independent components; these last steps are very fast (a few second). The whole procedure takes between 30 minutes and an hour.

Results

This section presents the results of each considered methodology in terms of prediction performances, for the two datasets.

Tecator

The ICA leads to a decomposition of 12 independent components for the Tecator dataset. In the case of this dataset, it must be mentioned that the number of the independent components has been limited to 12, because large instabilities in the results of the ICA algorithm were observed when it was tried to extract more independent components. However with 12 independent components, the reconstruction error of the spectra was around 35%, which is not good. This might be seen as a severe limitation in this case: to reach around 1% of reconstruction error, around 30 independent components would have been needed, inducing the instability problems. Nevertheless, it must be mentioned that even if the reconstruction error of the spectra is large, good prediction results might be obtained. Indeed, it might happen that the lost information in the reconstruction is not relevant for the prediction. In other words, for example, the prediction model will use only small ranges in the spectra, while completely discarding other parts that consequently do not have to be reconstructed correctly. It is therefore important to check, a posteriori, if the quantitative results obtained with the method are comparable (even if not equal or better) to those obtained by methods using all variables, regardless of the reconstruction error. This will be done below in Table 1.

Among the 172 spectra used for learning, two subsets of 57 spectra and a subset of 58 spectra are used for the 3-fold cross-validation procedure. Table 1 presents the NMSE on the test set (43 spectra) for the three considered approaches, as well as the number of variables involved in the LS-SVM (in the case of ICA and clustering, this means the number of independent components and clusters, respectively). The best prediction performances are achieved by the clustering methodology, closely followed by the "classical" model. The ICA approach give acceptable but worse results, while PLSR gives the worst results. It should be noted that all results, even the worst ones, correspond to good prediction performances. This means that interpreting the selected variables should be meaningful. In both case the functional projection allows one to reduce the number of variables significantly.

Figure 3 shows in grey the ranges of wavelengths selected by the clustering methodology; the vertical lines correspond to the limits of the clusters. The selected independent components provided by the ICA approach are indicated in Figure 4. As expected, they contain a limited range of spectral variables, while of course, as detailed in the introduction of this paper, they do not strictly correspond to pure spectra of pure components. The two proposed approaches select wavelengths around 930 nm, which corresponds to a bump in the shape of the spectra. This result is in agreement with the results e.g. in [START_REF] Rossi | [END_REF]. In addition, the large size of the cluster covering this wavelengths range tends to show that the whole bump in the shape of the spectra plays the same role in the prediction of the meat samples fat content. Both ICA and clustering methodologies select wavelengths in the neighborhood of 880-910 nm, which is also in conformity with what can be found in the literature. However, the first independent component is quite different from the wavelengths usually selected. This is probably due to the change of the baseline around 1040 nm.

Method

As a comparison, Figures 5 and6 give the loadings of the PLSR. While they show some interesting trends, they are obviously far more difficult to interpret than the ICA components (not to mention the wavelength cluster). This is mainly due to the fact that most of the coefficients of the loadings are non zero.

There are of course some wavelength intervals with zero influence on some specific loadings (e.g., [975,1050] for the third loading on Figure 5), but none of them is as focused as ICA components or wavelength clusters.

On the Tecator dataset, the functional approaches lead to good prediction quality based on localized wavelength ranges, something that can not be achieved by PLSR. Similar results were obtained in [START_REF] Rossi | [END_REF] by selecting individual spectral variables but at the expense of a much higher computation cost (roughly 500 time more processing time, i.e. more than 10 day of calculation compared to half an hour on identical hardware).

Wine

In the case of the Wine database, the ICA leads to 30 independent components. For this example, the number of independent components has been chosen only according to a reconstruction error argument, as no problem of ICA instability was observed. With 30 components, the reconstruction error goes below 1%, and does not significantly reduce if more components are used. On the contrary, using less than 30 The case of ICA is less clear as its predictive performances are not very good. Nevertheless, from the five selected components, only the third one does not correspond to a part of the spectra which is selected by the clustering approach, while the second component is a mixed between some wavenumbers around 3800 cm -1 and a large bump centered on 2600 cm -1 . Moreover, parts of the spectra defined by the second and fifth independent components are close to wavenumbers selected in [START_REF] Benoudjit | [END_REF], [24], in which they gave good predictive performances. Therefore, despite bad predictive performances, the components remain associated to interesting part of the spectra. A possible explanation of the differences between the clustering approach and ICA on this dataset is the very thin aspect of the cluster retained by the first method (which is a consequence of the high number of clusters originally retained by the linear selection method). As it is partly supervised, it can detect that a rather high precision is needed in some parts of the spectra in order to achieve good predictive performances. On the contrary, ICA components appear to have a rather large support (compared to the clusters) which are then incompatible with high quality prediction. This explanation is partly confirmed by the spiky aspect of the PLS loadings. On the contrary, the Tecator dataset corresponds to smoother loadings and to wider variable clusters. As already mentioned, this illustrates the need to compare the results to state-of-the-art methods using all variables (and not looking for interpretability).

Method

Conclusion

Reducing the number of variables taken into account in a prediction model for spectrometric data is not only important to avoid problems related to the curse of dimensionality (such as overfitting or difficulties of convergence), but also helps to increase the interpretability. Identifying the parts of the spectra that play a role is indeed a crucial issue for the practitioners.

Existing methods projecting the spectra into a lower-dimensional space present several drawbacks: the projection of the spectra on a basis such as B-splines is not adapted to the shape of the spectra and does not exploit the information on the parameter to predict that can be found in the spectral variables. In contrast, the clustering approach developed by Van Dijk and Van Hulle in [START_REF]Conference on Artificial Neural Networks[END_REF] is data-driven, but does not exploit the functional character of the spectra. As a result, clusters are allowed to contain non-consecutive spectral variables, which may be a problem from the interpretability point of view.

In this paper, we suggest a methodology to address the variable projection problem while avoiding the drawbacks pinpointed above. This approach is based on a functional clustering or the projection of the spectra on independent components. The projected variables are then selected according to the mutual information between them and the parameter of interest.

Both methods result in limited ranges of spectra used for the prediction model. These ranges may be exploited by practitioners to understand on a chemical point of view which pure components might be mostly responsible for the performances of the prediction model. It is not expected however to retrieve the spectra themselves of the pure components, as this would require further independence hypotheses that are not met in practice.

It must be mentioned that the goal of the methods presented in this paper is qualitative, in the sense that finding limited ranges of spectral variables is the objective. Quantitatively, the performances of the models may vary considerably. This is due to the fact that the procedure is designed as a compromise between computational burden and performances. In some situations, it may occur that state-of-the-art prediction performances would require more variables, or more clusters than affordable by the method. When using the methods presented in this paper, it is therefore of good practice to compare their quantitative results (prediction performances) with state-of-the-art (linear and nonlinear) methods that do not look for interpretability. The closest the performances are from the latter, the best confidence can be put in the quality of the resulting interpretation.

Experiments on the Tecator and Wine databases show that the models built on the projected variables can give accurate predictions. The projections lead to a significant reduction of the number of variables to 36 take into account in the models. Moreover, the selected wavelengths ranges correspond mainly to wavelengths identified as meaningful in the literature.

Clustering and ICA are given here as examples of functional preprocessing (projection), aimed at considerably reducing the number of variables used for the prediction, while exploiting the functional character of spectra for increasing the interpretability of the results. Further work will consider other types of functional projection, and how to choose between the different projection possibilities to obtain a good compromise between prediction performances and interpretability.

4. 1 .

 1 ICAIndependent Component Analysis (ICA)[START_REF]Wiley Series on Adaptative Learning Systems for Signal Processing[END_REF] is a signal processing technique aiming at separating measured signals into their "sources" components. Sources are considered as the basis components which are mixed to produce the measured signals. The general problem of recovering sources from the measured signals is called Blind Source Separation (BSS). It is blind because the mixture process (from the sources to the measured signals) is unknown and unobserved. To be solved, this problem required further hypotheses.

  instantaneous. Let us define S the kxn matrix of source signals s i , Z the pxn matrix of observed signals z i , and A the p×k mixing matrix, where k is the number of independent sources, p the number of observed signals, and n the number of observations per signal. Then the model may be written as . S A Z = The ICA problem reduces in finding estimations S ˆ of the sources by combining the observed signals

Figure 1 :

 1 Figure 1: Absolute Correlation between a spectral variable and all the other variables for the Tecator dataset.

  .

Figure 2 :

 2 Figure 2: Summary of data processing methods.

Figure 3 :

 3 Figure 3: Tecator spectra and spectral variables selected by the clustering approach (in grey).

Figure 4 :

 4 Figure 4: Selected independent components for Tecator.

Figure 5 :

 5 Figure 5: The six first selected loadings for PLSR in the case of Tecator.

Figure 6 :

 6 Figure 6: The five last selected loadings for PLSR in the case of Tecator.

Figure 7 :

 7 Figure 7: Wine spectra and spectral variables selected by the clustering approach (in grey).

Figure 8 :

 8 Figure 8: Selected independent components for Wine.

Figure 9 :

 9 Figure 9: The four first selected loadings for PLSR in the case of Wine.

Figure 10 :

 10 Figure 10: The four last selected loadings for PLSR in the case of Wine.

Table 1 :

 1 Prediction performances for Tecator.

		Number of (latent) variables	NMSE test
	PLSR on raw data	11	0.0170
	LS-SVM	100	0.0086
	ICA + MI + LS-SVM	7	0.0135
	Clustering + MI + LS-SVM	5	0.0077

Table 2 :

 2 Prediction performances for Wine.

		Number of (latent) variables	NMSE test
	PLSR on raw data	8	0.0058
	LS-SVM	256	0.0042
	ICA + MI + LS-SVM	5	0.3777
	Clustering + MI + LS-SVM	3	0.0086

Acknowledgments

The authors thank the anonymous referees for their valuable suggestions that helped to improve this paper. C. Krier is funded by a Belgian FRIA grant. Parts of this research results from the Belgian Program on Interuniversity Attraction Poles, initiated by the Belgian Federal Science Policy Office. The scientific responsibility rests with its authors.

1

independent components makes the reconstruction error increase significantly (10% already with 27 components).

Among the 91 spectra used for learning, two subsets of 30 spectra and one subset of 31 spectra are used for the 3-fold cross-validation procedure. The NMSE on the test set (30 spectra) and the number of variables used by the LS-SVM are detailed in Table 2, together with the result of PLSR. The lines of matrix A are normalized (reduction and centering) before building the LS-SVM. The best performances in prediction are obtained by the LS-SVM model build on all spectral variables, followed by PLSR and by the clustering approach. Those three methods give good prediction performances. As for the Tecator dataset, this means that the wavenumber ranges selected by the clustering approach should be meaningful. On the contrary, the ICA approach has rather bad performances: this gives an indication that care must be exercised when interpreting ICA components in this case. As for the Tecator dataset, the functional approach reduces the number of variables used by the nonlinear model. Figure 7 indicates in grey the ranges of wavenumbers selected by the clustering methodology. The vertical lines delimitate the clusters. Figure 8 shows the selected independent component. As expected, the wavenumber ranges selected by the clustering approach are very satisfactory in the context of alcohol concentration prediction: they correspond to the absorption range of the O-H bond present in alcohol (around 3600 cm -1 ). While PLS loadings lead to similar conclusion (see Figures 9 and10), this is far less clear in their case. Indeed most of the loadings associate high coefficients to wavenumbers around 3600 cm - 1 , but they also give an as large weight to wavenumbers around 1500 cm -1 . More generally, while the situation is better for PLS than in the case of the Tecator dataset, loadings remain spread over the whole spectral range and are therefore more difficult to interpret than cluster of wavenumbers.