Yaoyao Guoying Dong

Fiona Zhao

Numerical and Experimental Investigation of the Joint Stiffness in Lattice Structures Fabricated by Additive Manufacturing

Keywords: Additive Manufacturing, Finite Element Analysis, Lattice Structure, Mechanical Property

Highlight

 A joint stiffening element is proposed to simulate lattice structures fabricated by Additive Manufacturing.  The influence of the joint on the stiffness of the lattice structure can be considered in the proposed model.  A parametric study is conducted to quantify the influence of the joint on a certain lattice structure and to find the optimal parameters of the joint stiffening element.  Three-point bend testing is conducted, and it shows that the accuracy is significantly improved by the proposed model for bending dominant lattice structures.  The computational cost of the proposed model is much lower than the tetrahedral element model.

Introduction

The lattice structure in this paper is defined as a type of cellular architecture that has interconnected struts and joints in a three-dimensional (3D) space [START_REF] Dong | A Survey of Modeling of Lattice Structures Fabricated by Additive Manufacturing[END_REF]. With the development of Additive Manufacturing (AM), the fabrication of mesoscale lattice structures has become feasible. The lattice structure can achieve unique physical properties by adjusting the distribution of materials in the design space, thus it has promising applications in the industry. The lattice structure has been widely used as a lightweight material to reduce the weight or increase the stiffness of the structure [START_REF] Tang | Lattice Structure Design and Optimization With Additive Manufacturing Constraints[END_REF][START_REF] Alzahrani | Design of truss-like cellular structures using relative density mapping method[END_REF]. It was also designed for energy absorption purpose. Tancogne-Dejean et al. [START_REF] Tancogne-Dejean | Additively-manufactured metallic micro-lattice materials for high specific energy absorption under static and dynamic loading[END_REF] found that the octet lattice material featured a constant plateau stress between initial yield and densification as well as an exceptionally high specific energy absorption. Furthermore, the lattice structure has also been used for heat dissipation due to its high surface-to-volume ratio [START_REF] Brooks | Design of conformal cooling layers with self-supporting lattices for additively manufactured tooling[END_REF].

The physical properties of lattice structures have been studied for decades. In the early stage, Gibson et al. [START_REF] Gibson | The mechanics of two-dimensional cellular materials[END_REF][START_REF] Gibson | The mechanics of three-dimensional cellular materials[END_REF] investigated the mechanical properties of two-dimensional (2D) and 3D cellular materials with beam theory. It was found that the effective Young's modulus was controlled by the relative density of the cellular material. After AM started to be used to fabricated lattice structures, more and more complicated lattice structures were manufactured and tested through experiments. Labeas and Sunaric [START_REF] Labeas | Investigation on the Static Response and Failure Process of Metallic Open Lattice Cellular Structures[END_REF] conducted compressive experiments to understand the static response and failure process of lattice structures made of Stainless Steel 316L. The result showed that the structural response was strongly influenced by the strut geometrical characteristics such as the aspect ratio, the unit-cell size, and the shape. Alsalla et al. [START_REF] Alsalla | Fracture toughness and tensile strength of 316L stainless steel cellular lattice structures manufactured using the selective laser melting technique[END_REF] investigated the fracture toughness of Stainless steel 316L lattice structures through tensile tests. It was found that the result was highly influenced by the building orientation. The horizontal struts were weaker than the struts built in the vertical direction. Yang et al. [START_REF] Yang | Non-stochastic Ti-6Al-4V foam structures with negative Poisson's ratio[END_REF] investigated the flexural property of auxetic lattice structures through bending tests. The result indicated that the flexural strength of auxetic lattice structures was higher than regular sandwich panel structures. Furthermore, impact loading tests [START_REF] Winter | Plate-impact loading of cellular structures formed by selective laser melting[END_REF], dynamic tests [START_REF] Salehian | Dynamic analysis of a lattice structure by homogenization: Experimental validation[END_REF], and fatigue test [START_REF] Jamshidinia | Fatigue properties of a dental implant produced by electron beam melting (EBM)[END_REF] have also been conducted to investigate the performance of lattice structures fabricated by different types of AM techniques.

Because the manufacturing cost and the time cost of most AM processes are still high, simulation methods are widely used to predict the mechanical performance of lattice structures, which accelerates the design process and saves the cost. Simulation methods are also important in the optimization process of lattice structures. Homogenization Method (HM) and Finite Element Method (FEM) are the most frequently used methods to simulate the mechanical property of lattice structures. Due to the manufacturing defects of as-build lattice structures, several methods have been used to improve the simulation method in order to more accurately predict the mechanical performance of lattice structures. Park et al. [START_REF] Park | Effective mechanical properties of lattice material fabricated by material extrusion additive manufacturing[END_REF] proposed a two-step homogenization method considering the manufacturing influence of AM. Firstly, the effective structural element parameters were determined by a voxel-based model that included the shape deviations in the as-fabricated struts. Secondly, these parameters were implemented in the homogenization model to calculate the effective mechanical property of the lattice structure. It was found that the result obtained from the two-step approach was closer to the experimental result compared to the result from the direct homogenization method. Therefore, the shape variation caused by the manufacturing process is a crucial factor that should be considered in the simulation process. Furthermore, the manufacturing defects should also be considered in FEM. Cahill et al. [START_REF] Cahill | Finite element predictions compared to experimental results for the effective modulus of bone tissue engineering scaffolds fabricated by selective laser sintering[END_REF] generated finite element models for lattice structures with both smooth surfaces and rough surfaces. The result revealed that the effective modulus of the lattice with smooth surfaces was 68% higher than that of the lattice with a rough surface. One way to improve the accuracy of the simulation result is to update the computer-aided design (CAD) model with the measurement from the fabricated sample [START_REF] Yang | Mechanical properties of 3D re-entrant honeycomb auxetic structures realized via additive manufacturing[END_REF]. Another way to consider the roughness of the lattice surfaces is to use beam elements [START_REF] Karamooz Ravari | Numerical investigation on mechanical properties of cellular lattice structures fabricated by fused deposition modeling[END_REF][START_REF] Campoli | Mechanical properties of open-cell metallic biomaterials manufactured using additive manufacturing[END_REF] with varying thickness to discretize the lattice strut. Furthermore, Micro-CT has been used to reconstruct the geometrical model of the as-built lattice structure. The reconstructed model can be used to generate the tetrahedral element model for FEM [START_REF] Suard | Towards stiffness prediction of cellular structures made by electron beam melting (EBM)[END_REF][START_REF] Sercombe | Failure modes in high strength and stiffness to weight scaffolds produced by Selective Laser Melting[END_REF].

Another aspect that needs to be considered in the simulation is the influence of the joint on the stiffness of the lattice structure. The concentration of the material in the nodal region will increase the cross-section area of the lattice strut. Labeas and Sunaric [START_REF] Labeas | Investigation on the Static Response and Failure Process of Metallic Open Lattice Cellular Structures[END_REF] increased the thickness of the beam element by 40% on 1/10 of the strut length at both strut ends to compensate for the underestimation of the stiffness. Smith et al. [START_REF] Smith | Finite element modelling of the compressive response of lattice structures manufactured using the selective laser melting technique[END_REF] used the same way to increase the stiffness of the strut at the joint area. Luxner et al. [START_REF] Luxner | Finite element modeling concepts and linear analyses of 3D regular open cell structures[END_REF] directly increased the stiffness of the beam element at the nodal region 1000 times than the stiffness of the strut. On the contrary, Park and Rosen [START_REF] Sang-In | Homogenization of Mechanical Properties for Additively Manufactured Periodic Lattice Structures Considering Joint Stiffening Effects[END_REF] proposed a homogenization procedure using semi-rigid joint frame elements for lattice structures. This semi-rigid joint method has been widely used in macro-scale frame structures in civil engineering [START_REF] Díaz | Optimum design of semi-rigid connections using metamodels[END_REF]. They mentioned that the frame structure could not achieve a perfectly rigid joint. Therefore, the strut was modeled using semi-rigid jointed frame elements that consist of a shear deformable frame region and two joint regions at both ends. The joints were modeled as point-connections in conventional frame elements. The elastic modulus obtained by the proposed model is less than the conventional frame model, which means the actual stiffness of the lattice structure was reduced by the joints.

Even though the influence of the joint of the stiffness of lattice structures has been considered in previous studies. It is not clear how much the joint will affect the stiffness of the strut. In this paper, the stiffness of the joint in the lattice structure will be investigated quantitatively. Then, a joint stiffening element is proposed to precisely estimate the influence of the joint on the stiffness of lattice structures. It can be used to build FEM models to simulate the mechanical performance of entire lattice structures. In Section 2, the influence of the joint on the stiffness is discussed first. Then, a joint stiffening element is proposed to model the lattice strut with the influence of joints. In Section 3, a parametric study is carried out to find the optimal coefficients for the joint stiffening element. A case study is presented in Section 4 to demonstrate the process of simulating lattice structures with the joint stiffening element. The as-fabricated material property is obtained by tensile specimens to eliminate the imperfection caused in the manufacturing process. The result from the proposed simulation model and the traditional models are compared with the experimental result. Finally, the conclusions and future research are summarized in Section 5.

The proposed joint stiffening element 2.1 Effects of joints on the stiffness of lattice structures

The joint in the lattice structure has the influence on the structural stiffness for two reasons. The first reason is that the size of the joint in the CAD model is relatively large compared to the strut thickness. The second reason is that the material tends to concentrate in the joint area during the fabrication by AM. In this paper, only the joint geometry designed by the CAD model is considered. A comparison of the joint sizes of two lattice structures is shown in Figure 1. The joint in Figure 1(a) is only the intersection of the cylindrical struts. But in Figure 1(b), there is a large joint used to connect the surrounding struts. In the geometrical modeling process of the lattice structure, there are lots of struts that need to be connected by the Boolean operation. Generally, joints are used to trim the intersection of lattice struts so that the difficulty of the Boolean operation can be reduced. However, it makes the joint size larger than the thickness of the strut. For the bending dominant lattice topologies, the stress will concentrate at the end of the strut. The large joint will affect the stiffness at the end of the strut. If the beam element is used to simulation the mechanical performance of the lattice structure, the total stiffness of the structure tends to be underestimated. Because in the traditional beam element model, the stiffness is uniform on the beam element.

Figure 1 The comparison of the joint sizes with different types of intersection

A unit cell of Vintiles lattice structure is simulated by both the solid element model and the beam element model to investigate the stiffness underestimation. The solid element is generated by the CAD model of the lattice structure. The geometry of the joints and the lattice struts are both meshed with the tetrahedral element. Therefore, in the solid element model, the influence of the joint on the stiffness of the lattice structure is considered. However, the beam element model uses Timoshenko beam elements to mesh the wireframe of the lattice structure. All the beam elements have the uniform thickness. The geometry of the joint is not considered in the beam element model. Some parameters of both simulation models are summarized in Table 1. Convergence tests are conducted for both the solid element model and the beam element model as shown in Figure 2. It is found that if the number of the solid element is small, the result is 42.5% higher than the converged result. Therefore, a large number of elements is required for an accurate simulation result of the 3D solid element model. However, when the lattice structure has a lot of unit cells, it will be extremely difficult to generate the 3D mesh of the structure. The meshing algorithm of the commercial software may fail in the mesh generation. This is one reason that the beam element model is preferred to simulate complicated lattice structure. The comparison of the converged simulation results is shown in Figure 2(c). The total CPU time of the solid element model is 19.2s (119754 elements) while that of the beam element model is 0.3s (720 elements). However, it should be noted that this is only a single unit cell. If the lattice structure has plenty of unit cells, the computational cost of the solid element model will be significantly high. For example, if the lattice structure has 10 unit cells along each axis, the total CPU time of the solid element will be more than 5 hours. But the time of the beam element model will be only 5 minutes. Therefore, beam element is the ideal element type to simulate complex lattice structures. The current problem of the beam element is that the stiffness of joints in the lattice structure is underestimated. It can be seen from Figure 2(c) that under the same loading condition, the effective elastic modulus of the solid element model is 1755MPa while that of the beam element model is only 1133MPa. This proves that the stiffness of the beam element model is much lower than that of the solid element model.

Therefore, the beam element model should be modified so that the influence of the joint can be considered in the simulation process. and are the deflections and the rotations at two ends of the element. and are the deflections and the rotations at the internal nodes of the element. Because the cross-section of the joint has a larger area than that of the lattice strut, the flexural stiffness of the joint is greater than the flexural stiffness of the strut. Therefore, a flexural coefficient is proposed to compensate for the stiffness of the joint. The force-displacement equation for the Timoshenko beam element is derived from the method in the literature [START_REF] Reddy | An introduction to the finite element method[END_REF]:

(

(

where is the Young's modulus, is the shear modulus, is the area of the strut cross section, is the moment of inertia, is the shear correction coefficient, and is the length of the beam element.

is the stiffness matrix of the Timoshenko beam element. is the flexural displacement and the rotation on the nodes, and is the transverse force and moment on the end nodes of the beam element.

To represent the length of the joint segment and the strut segment, a joint length coefficient is proposed. and can be calculated as: (6) [START_REF] Gibson | The mechanics of three-dimensional cellular materials[END_REF] where is the total length of the joint stiffening element. By substituting and into Eq. ( 2), the stiffness matrices of the joint segment and the strut segment can be obtained as:

(8) (9)
The stiffness matrix of the joint segment is multiplied by so that the stiffening effect of the joint can be considered in the simulation of lattice structures. The global stiffness matrix can be assembled from and . The total force-displacement equation for the joint stiffening element can be written as [START_REF] Yang | Non-stochastic Ti-6Al-4V foam structures with negative Poisson's ratio[END_REF] Since the degrees of freedom at the second node and the third node do not correspond to a degree of freedom of any other element, they can be eliminated to obtain the element stiffness matrix that corresponds to the degrees of freedom at the first and fourth nodes only. Static condensation is employed to reduce the number of element degrees of freedom [START_REF] Bathe | Finite element procedures[END_REF][START_REF] Qu | Model order reduction techniques with applications in finite element analysis[END_REF]. The elimination of the degrees of freedom at the second and third nodes is carried out using Gauss elimination. Eq. ( 10) can be partitioned as:

(11) (12)
For simplicity, Eq. ( 11) and Eq. ( 12) can be written as:

(13) (14)
The internal displacement and rotation can be calculated as: [START_REF] Cahill | Finite element predictions compared to experimental results for the effective modulus of bone tissue engineering scaffolds fabricated by selective laser sintering[END_REF] By substituting Eq. ( 15) into Eq. ( 14), the total force-displacement equation can be written as: [START_REF] Yang | Mechanical properties of 3D re-entrant honeycomb auxetic structures realized via additive manufacturing[END_REF] where [START_REF] Karamooz Ravari | Numerical investigation on mechanical properties of cellular lattice structures fabricated by fused deposition modeling[END_REF] The four-nodes joint stiffening element can be condensed into a two-node element by Eq. [START_REF] Karamooz Ravari | Numerical investigation on mechanical properties of cellular lattice structures fabricated by fused deposition modeling[END_REF].

is the flexural stiffness matrix of the joint stiffening element. There are two parameters and that control the joint property. The stiffness coefficient determines the amplification factor of the flexural stiffness at the joint. The length coefficient indicates the length of the joint area compared to the total length of the strut. The calculation of and will be discussed in Section 3.

Apart from the flexural stiffness, the axial stiffness of the lattice strut is also influenced by the joint in the lattice structure. The axial stiffness matrix of the joint stiffening element is derived as: [START_REF] Campoli | Mechanical properties of open-cell metallic biomaterials manufactured using additive manufacturing[END_REF] (19) [START_REF] Sercombe | Failure modes in high strength and stiffness to weight scaffolds produced by Selective Laser Melting[END_REF] where ( 21) [START_REF] Luxner | Finite element modeling concepts and linear analyses of 3D regular open cell structures[END_REF] where is the length coefficient for the joint segment under the axial load, is the axial stiffness coefficient of the joint segment. Because the formulation of the axial stiffness matrix is not related to the moment of inertia, which is different from the flexural stiffness matrix. The influence of the joint on the axial part and the flexural part are different. Therefore, new parameters and are used for the axial stiffness matrix.

Similar to the axial stiffness matrix, the torsional stiffness matrix of the joint stiffening element is derived as:

(23) (24) (25)
where is the torsion constant. Because both the moment of inertia and the torsion constant are proportional to the fourth power of the strut radius. The influence of the joint on the flexural stiffness and the torsional stiffness are similar. Therefore, the length coefficient and the stiffness coefficient for the flexural load and the torsional load are the same. The total forcedisplacement equation can be assembled from Eq. ( 16), Eq. ( 18), and Eq. ( 23) as:

(26)

Calculation of parameters

In the proposed joint stiffening element model, there are four coefficients evaluating the influence of the joint. They are the flexural stiffness coefficient ; the axial stiffness coefficient ; the length coefficient and . The determination of these parameters will be introduced in this section. The first step is to quantify the influence of the joint on the mechanical performance of the lattice strut by the tetrahedral element model. A convergence test is conducted to ensure the accuracy of the result. Secondly, the joint stiffening element model of the lattice strut is used to obtain a group of results with a range of parameters. Finally, one pair of parameters in the joint stiffening element model will be selected whose result is closest to the result in the tetrahedral element model.

Tetrahedral element model of the lattice strut

In order to get the tetrahedral element model of lattice strut, one lattice strut is cut from the lattice structure as shown in Figure 4. Other lattice struts connected to the joint are removed so that only one strut is left with joints at the two ends. The single lattice strut model is then meshed by 4-node tetrahedral elements. Different mesh sizes are used to generate multiple models for the convergence test. Because only a single strut is simulated, the computational cost is acceptable when the result converges. Examples of the numerical model are shown in Figure 5. Both numerical models are fixed at the left end. The numerical model in Figure 5(a) is forced with a unit transverse displacement to simulate the flexural behavior of the lattice strut. The numerical model in Figure 5(b) is forced with a unit axial displacement to simulation the tensile behavior of the lattice strut. To indicate the influence of the joint on the stiffness of the lattice strut, two beam element models are also generated with the same boundary condition as that in tetrahedral element models. The simulation results of the tetrahedral element models and the beam element models are shown in Figure 6. The total length of the strut in the simulation is 10. The relationship between the displacement and the axial position are demonstrated by the curves in Figure 6. It is found that there is a significant difference between the solid element model and the beam element model. The slope of the red curve is smaller than that of the blue curve at the two ends of the strut. This is because the stiffness at the two ends of the strut is much higher than the middle of the strut. The deformation should concentrate in the middle of the strut. The solid element model can reflect the actual stiffness of the joint while the beam element model does not consider the joint. Therefore, the strut is very strong at both ends in the solid element model. It causes the deformation at the ends of red curves (solid element model) smaller than that in the blue curves (beam element model).

Determination of the optimal parameters

The comparison of the results in Figure 6 shows that the stiffness at the end of the beam element should be increased to match the actual stiffness of the joint. The next step is to determine the stiffness and the length coefficients of the joint stiffening element. Firstly, a range of coefficients is selected to generate the joint stiffening element model of the lattice strut. Then, a group of results can be obtained from the joint stiffening element models. To evaluate the accuracy of selected coefficients, n points with position on the lattice strut will be chosen. The error between the solid element model and the joint stiffening element model is calculated by [START_REF] Qu | Model order reduction techniques with applications in finite element analysis[END_REF] where is displacement of the solid element model at the position , is the displacement of the joint stiffening element model. If the error is small, it means the result of the joint stiffening element is close to the solid element result.

A series of stiffness and length coefficients are generated for the flexural load and the axial load. Each combination of the stiffness and length coefficients has been used to calculate the deformation of the lattice strut. Then, the error between the solid element model and the joint stiffening element model are calculated from Eq. [START_REF] Qu | Model order reduction techniques with applications in finite element analysis[END_REF]. The relationships between the error and the coefficients are shown in Figure 7. In this example, is from 1 to 100, is from 0.1 to 0.2, is from 1 to 20, and is from 0.05 to 0.2. Each parameter has 50 items within its range. When the combination of the stiffness coefficient and the length coefficient gives the minimum error, they will be selected as the optimal parameters for the joint stiffening element.

Validation of selected parameters

After the determination of the optimal coefficients, these coefficients are imported into the joint stiffening element. The result of the optimized joint stiffening element is compared to the solid element model as shown in Figure 8. It is found that the results of both models coincide, which means the displacement in the joint stiffening element is nearly the same as the displacement of the solid element model. Therefore, the joint stiffening element can accurately reflect the influence of the joint on the stiffness of the lattice strut. Because the beam element does not consider the influence of the joint, it underestimates the stiffness of the lattice structure. It is more accurate to simulate the lattice structure with the proposed model than using the beam element model alone.

Case study

To validate the proposed joint stiffening element, two lattice structures with the Cubic-center topology and the Vintiles topology are designed for the three-point-bend testing as shown in Figure 9. Both lattice structures have cells along x, y, and z-direction. The dimension of the unit cell is 10mm in all three axes. The relative density of both structures is 15%. The Cubic-center topology is stretching dominant while the Vintiles topology is bending dominant. It is assumed that the joint will have more influence on the bending dominant topology than that of the stretching dominant topology. Because the stress is concentrated at the end of the beam element if it is under bending load. The underestimation of the joint stiffness will be magnified by the stress concentration. This assumption will also be validated in this case study.

As-fabricated material property

The lattice structures are fabricated by Fused Deposition Modeling (FDM) process with ABS materials. In order to obtain the bulk material properties more accurately, several tensile parts are designed as shown in Figure 9(c). Two specimens are designed so that the testing strut can be put parallel, inclined, and perpendicular to the build platform. The lattice samples and the tensile specimens are fabricated together in one batch as shown in Figure 10(a). The asfabricated samples are shown in Figure 10(b). Because there are several factors in the FDM process, such as the toolpath and the layer adhesion that may have the influence on the material property. The as-fabricated material tends to be different from the bulk material. The influence from the manufacturing process should be eliminated as much as possible so that the joint stiffening effect can be isolated. Otherwise, the manufacturing influence may reduce the stiffness of the lattice structure, which may counteract the joint stiffening effect and lead to wrong results. Therefore, the tensile specimens are designed with struts that have the diameter close to the strut thickness in the lattice structure. All the parameters used in the fabrication are the same for both the tensile specimen and the lattice structure. The tool path in the lattice strut is similar to the tool path in the specimen. The strength of the adhesion between layers can be investigated by the vertical specimen. The strain rate of the tensile test is 1mm/min. The load cell is 5kN. A strain gauge is used to measure the strain of the specimen. The averaged testing results and the standard deviations of the tensile samples are shown in Figure 11. It is found that the elastic modulus of the horizontal specimen and the inclined specimen are very close. But the elastic modulus of the vertical specimen is 20% lower than the other two specimens. It is concluded that the adhesion between layers will reduce the elastic modulus of as-fabricated parts. Then, the as-fabricated material property will be imported into the simulation model. The lattice strut with different orientation will be assigned with different properties. This can minimize the inaccuracy that comes from the material property. The elastic modulus and the flexural modulus are also tested by ASTM D638 [START_REF]Standard Test Method for Tensile Properties of Plastics[END_REF] and ASTM D790 [START_REF]Standard Test Methods for Flexural Properties of Unreinforced and Reinforced Plastics and Electrical Insulating Materials[END_REF] for comparison as shown in Figure 11. Because the tool path has influence on the mechanical property of the FDM fabricated parts, the elastic modulus of these samples is different. Only the proposed test samples have the same tool path with that of the lattice strut. It is more reasonable to use the elastic modulus of the proposed tensile specimens in the simulation process.

Determination of optimal coefficients

The next step is to determine the coefficients for the joint stiffening element models of both lattice structures. Firstly, the strut models of these two lattice structures are generated with tetrahedral elements similar to the models shown in Figure 5. The unit flexural displacement load and the unit axial displacement load are applied to the strut models. Secondly, a series of flexural coefficients and axial coefficients are imported into the joint stiffening element models. Then, Eq. ( 27) is used to calculate the error between the tetrahedral element model and the joint stiffening element model. When the error gets the minimum value, the corresponding coefficients are selected as the optimal parameters. The optimal coefficients for the flexural load and the axial load are summarized in Table 2. The optimal coefficients are used in the joint stiffening element model to simulate the threepoint bending performance of the designed lattice structures. To compare the result of the proposed models, the designed lattice structures are also meshed by Timoshenko beam elements and tetrahedral elements. The Timoshenko beam element model is generated by the lattice wireframe, which means it doesn't consider the influence of the joint on the stiffness of the structure. By contrast, the tetrahedral element is generated based on the geometrical model of the lattice structure, which takes the joints into consideration.

Results and Discussion

Physical tests are used to validate the simulation results. The three-point-bend testing is conducted on the TestResources 313 universal tester as shown in Figure 12. The loaddisplacement curves of both lattice structures are shown in Figure 13. Each of them is fabricated and tested three times to reduce the uncertainty of the FDM machine. The average stiffness and the standard deviation of stiffness is calculated for both lattice structures. It can be found in the experiment that the stretching dominant lattice structure tends to exhibit a higher stiffness than the bending dominant lattice when they have the same relative density. But deformation of the stretching dominant lattice at the ultimate load is smaller than that of the bending dominant lattice. The load-displacement curves in Figure 13 show that the Cubic-center lattice structure fails earlier than the Vintiles lattice structure. This is because the Cubic-center topology is stretching dominant while the Vintiles topology is bending dominant. The load on the strut of stretching dominant lattices is mainly the axial force. But in the bending dominant lattice structure, most of the struts are under flexural load. In order to compare experimental results with simulation results, the average stiffness and the standard deviation of the as-fabricated lattice structure is calculated from the loaddisplacement curves. The results of the stiffness obtained by both the simulation method and the experimental method are summarized in Figure 14. The tetrahedral mesh in the simulation model is generated by Altair HyperWorks. All the simulation models except the proposed model are solved by Abaqus. Four-node tetrahedral element (C3D4) is used in the solid element model. Two-node beam element (B31) is used in the beam element model and the models in the literature. A self-contained beam element solver is written in Matlab to solve the proposed model. The parameters used in the joint stiffening element are integrated into the Matlab code. The parameters in Table 2 are the inputs of the Matlab code.

For the Cubic-center lattice structure, the difference between the beam element model and the solid element model is relatively small. However, the stiffness of the Vintiles lattice structure predicted by the beam element model is much lower than that predicted by the solid element model. It means that the joint stiffening effect is more obvious on the bending dominant topology than the stretching dominant topology. It can be explained by the formulation of the stiffness matrix of the beam element. The bending part of the stiffness matrix includes the moment of inertia, which is proportional to the fourth square of the radius. The axial part of the stiffness only includes the cross-section area, which is proportional to the square of the radius. If the radius at the joint area is increased, the bending part is more positively influenced. The deformation of the struts in Vintiles lattice is mainly contributed by the bending deformation. But the deformation of the struts in Cubic-center lattice is mainly contributed by the axial deformation. Therefore, the Vintiles lattice structure is more positively influenced by the joint. The second reason is that the stress in the strut of Vintiles lattice concentrates at the joint area. But the stress is more uniformly distributed in the strut of Cubic-center lattice structure. Since the cross-section area of the joint is larger than that of the strut. It can be considered as adding materials in the stress concentration area. Therefore, it has more influence on the stiffness of the Vintiles lattice structure.

The result also shows that the joint stiffening element model can increase the estimated stiffness compared to the beam element model. The result of the Vintiles lattice structure is doubled when the stiffness of the joint is considered in the proposed model. Furthermore, it is found that, compared to the results of the beam element and tetrahedral element models, the results in the proposed model are closer to the experimental results for both topologies, especially for the Vintiles topology. It can be concluded that the proposed model is able to predict the mechanical performance of lattice structures more accurately. The joint stiffening effect must be considered in bending dominant topologies. Otherwise, the stiffness of the structure will be seriously underestimated by the beam element model.

The proposed model is also compared to existing methods in the literature that consider the increase of the joint stiffness as shown in Figure 14. One method is that 1/10 of the strut at the end are modeled with a 40% higher cross-section area [START_REF] Labeas | Investigation on the Static Response and Failure Process of Metallic Open Lattice Cellular Structures[END_REF]. Another method is to increase the stiffness of the beam element 1000 times inside a sphere centered at the joint that has the same radius as the strut [START_REF] Luxner | Finite element modeling concepts and linear analyses of 3D regular open cell structures[END_REF]. Even though these two methods can improve the stiffness of the lattice structure, they have a large deviation from the experimental result. It is found that the first method underestimates the stiffness while the second method overestimates the stiffness. This is because the parameters in their methods are determined by experience. But the proposed model has a systemic method to optimize parameters used in the joint stiffening element, which lead to a more accurate result in the proposed model. Furthermore, both models in the literature require a larger quantity of the beam element than the proposed model. The quantity of mesh is around 10 times than the proposed model.

The advantage of the proposed model compared to the solid element model is that the computational cost of the proposed model is much lower than the solid element model. Because the geometry of the lattice structure is very complicated with a lot of small features. It will lead to much difficulty in the generation of the 3D tetrahedron mesh. The size of each tetrahedron element needs to be very small to improve the mesh quality, which means the total number of the is very large. It can be found in Figure 14 that the coarse tetrahedral mesh model of the Vintiles lattice structure does not obtain an accurate result. Even though the accuracy of the solid element model of the Vintiles lattice structure can be improved by using a finer mesh. The computational cost is much higher than that of the other simulation models.

The comparison of the quantity of the mesh between the solid element model and the proposed element model is shown in Table 3. A large number of meshes means the degrees of freedom of the structure is high. It takes more memory and time to solve the linear system of the structure. Instead, the beam element can mesh the lattice structure by using the wireframe information.

Each lattice strut can be modeled with one beam element. This will significantly reduce the computational cost in the simulation of the mechanical performance of lattice structures. The computational time to solve the proposed model and solid element models are also summarized in Table 3. The fine mesh model takes around one hour to solve the FEA equations while the proposed model only takes a second. It should be noted that the time of the meshing process has not been included yet. The time to generate the 3D solid element is also much longer than that to generate the joint stiffening element model.

In this paper, the joint stiffening element is an improved beam element that can compensate for the joint stiffening effect. It also has a much lower computational cost than the solid element model. It can accelerate the design process of lattice structures. It can also give designers more accurate results when simulating the mechanical performance of lattice structures.

Conclusions

In this paper, a concept called joint stiffening element is proposed to investigate the influence of the joint on the stiffness property of lattice structures. The joint stiffening element is divided into three segments. The segments at two ends represent the joint area in the lattice strut. The middle segment represents the actual strut length of the lattice structure. Four parameters are defined in the proposed element to evaluate the stiffening effect of the joint on a single lattice strut. Then, the stiffness matrix of the proposed element is derived by rearranging the displacement vector and the load vector. A parametric study has been conducted to determine the optimal parameters of the joint stiffening element for a specific lattice structure. The unit axial load and transverse load are added on the lattice strut to investigate the influence of the joint on the tensile and flexural properties.

To validate the proposed model, the lattice structures with Cubic-center and Vintiles topologies have been fabricated for the three-point-bend testing. Tensile specimens are tested to obtain the as-fabricated material property in order to reduce the influence caused by the manufacturing process. It is found that the elastic modulus of the vertical strut is 20% lower than that of the horizontal strut and inclined strut. The lattice structure is simulated using the beam element model, the solid element model, the proposed model, and the models in literature. It is found that the influence of the joint is more obvious on the bending dominant topology than the stretching dominant topology. For the Vintiles lattice structure, the stiffness predicted by the beam element model is only 52.7% of the experimental result. However, the proposed model that considers the joint influence can accurately estimate the stiffness of the structure. For the Cubic-center lattice structure, the difference between the result of the beam element model and the experimental result is only 12.5%. This is because most struts in Cubiccenter lattice structure are under axial load instead of bending load, in which case the joint stiffening effect is not that significant. Nevertheless, the result obtained by the proposed model is still closer to the experimental result than the traditional beam element model. Future research will focus on the feasibility of the joint stiffening element on the heterogeneous lattice structures. Because the joint size is nonuniform in heterogeneous lattice structures, a relationship between the parameter and the joint size should be established. The application of the joint stiffening element in the nonlinear analysis will also be investigated to simulate the plastic deformation of lattice structures.

Figure 2

 2 Figure 2 The simulation results with the convergence test to show the influence of the joint2.2 The proposed concept of the Joint stiffening elementTo consider the joint stiffening effect on the lattice structure, a joint stiffening element is proposed. It is modified from the Timoshenko beam element. The conceptual configuration of the joint stiffening element is shown in Figure3. The joint stiffening element consists of three segments. Two segments at the end of the element represent the joint area in the lattice structure. The middle segment of the joint stiffening element represents the actual length of the lattice strut.and are the length of the joint stiffening segment and the strut segment, respectively. and are the deflections and the rotations at two ends of the element. and are the deflections and the rotations at the internal nodes of the element. Because the cross-section of the joint has a larger area than that of the lattice strut, the flexural stiffness of the joint is greater than the flexural stiffness of the strut. Therefore, a flexural coefficient is proposed to compensate for the stiffness of the joint.

Figure 3

 3 Figure 3 The conceptual configuration of the joint stiffening element.

Figure 4

 4 Figure 4 The geometrical model of the joint stiffening element inside a lattice structure.

Figure 5

 5 Figure 5 Numerical model and the boundary conditions of the single lattice strut.

Figure 6

 6 Figure 6 Comparison of the simulation results of the solid element model and the beam element model without considering the joint stiffening effect.

Figure 7

 7 Figure 7 The relationship between the error and the coefficients

Figure 8

 8 Figure 8 The comparison of the optimized joint stiffening element result and the solid element result

Figure 9

 9 Figure 9 The geometrical model of the three-point bend testing, (a) Cubic-center lattice structure, (b) Vintiles lattice structure, and (c) tensile specimens used to obtain the elastic modulus of the bulk material.

Figure 10 (

 10 Figure 10 (a) Printing set-up for the lattice structures and the tensile specimens. (b) as-built samples.

Figure 11

 11 Figure 11 The elastic modulus of the tensile specimens and the flexural modulus of the bending specimens.

Figure 12

 12 Figure 12 Mechanical test on (a) the tensile specimen and (b) the three-point bending lattice structure.

Figure 13

 13 Figure 13 The load-displacement curves of the Cubic-center and Vintiles lattice structures, σ is the standard deviation of the stiffness

Table 1

 1 The parameters in the simulation models

	Topology	Cell size	Number of	Strut	Young's	Poisson's
			cells on each	thickness	modulus	ratio
			axis
	Vintiles	50mm	1	8mm	150GPa	0.3

Table 2

 2 Optimal coefficients for the designed lattice structures

		Flexural coefficient	Axial coefficient
	Cubic-center	31.5	0.155	2.16	0.162
	Vintiles	4.6	0.213	2.36	0.183

Table 3

 3 The comparison of the computational costs between the solid element model and the proposed model

		Cubic-center			Vintiles
		Element	DOFs	Time(s) Element	DOFs	Time(s)
	3D element (coarse)	369440	393087	61	252149	275109	36
	3D element (fine)	15529031 11057385	3397 15486480 10719162	3686
	Proposed model	2040	2796	0.25	4356	14724	0.56

Acknowledgment

This research work is supported by National Sciences and Engineering Research Council of Canada Discovery Grant RGPIN436055-2013.


Guilherme Dean
email: guidean@unicamp.br

Miguel Couceiro
email: miguel.couceiro@loria.fr

Leonardo Tomazeli Duarte
email: leonardo.duarte@fca.unicamp.br

A STATISTICAL APPROACH TO DETECT SENSITIVE FEATURES IN A GROUP FAIRNESS SETTING A PREPRINT

The use of machine learning models in decision support systems with high societal impact raised concerns about unfair (disparate) results for different groups of people. When evaluating such unfair decisions, one generally relies on predefined groups that are determined by a set of features that are considered sensitive. However, such an approach is subjective and does not guarantee that these features are the only ones to be considered as sensitive nor that they entail unfair (disparate) outcomes.

In this paper, we propose a preprocessing step to address the task of automatically recognizing sensitive features that does not require a trained model to verify unfair results. Our proposal is based on the Hilbert-Schmidt independence criterion, which measures the statistical dependence of variable distributions. We hypothesize that if the dependence between the label vector and a candidate is high for a sensitive feature, then the information provided by this feature will entail disparate performance measures between groups. Our empirical results attest our hypothesis and show that several features considered as sensitive in the literature do not necessarily entail disparate (unfair) results.

1 These methods combine different fairness interventions such as in-processing and post-processing.

Introduction

The use of Machine Learning (ML) methods to deal with real-world problems has grown exponentially in the last decades [START_REF] Brink | Real-world machine learning[END_REF][START_REF] Iqbal | Machine learning: Algorithms, real-world applications and research directions[END_REF]. In particular, ML is currently integrated in every decision support (DS) system with critical societal impacts. This fact has raised increasing concerns about the (un)fairness of such 'artificial deciders' in several contexts [START_REF] Ninareh | A survey on bias and fairness in machine learning[END_REF], such as in recidivism risk prediction [START_REF] Angwin | Machine Bias -ProPublica[END_REF], facial recognition systems [START_REF] Inioluwa | Actionable auditing: Investigating the impact of publicly naming biased performance results of commercial AI products[END_REF], platforms for job applications [START_REF] Hangartner | Monitoring hiring discrimination through online recruitment platforms[END_REF], hate speech and abusive language detection [START_REF] Davidson | Racial bias in hate speech and abusive language detection datasets[END_REF]. This motivated several efforts towards both the detection and mitigation of unfairness in ML and DSs. The former are mostly centered in the detection of undesirable bias, while the later seek to reduce such biases. These methods fall into four main categories, namely, pre-processing [START_REF] Calders | Building classifiers with independency constraints[END_REF][START_REF] Roh | Fairbatch: Batch selection for model fairness[END_REF][START_REF] Dean Pelegrina | Analysis of trade-offs in fair principal component analysis based on multi-objective optimization[END_REF][START_REF] Dean | A novel approach for Fair Principal Component Analysis based on eigendecomposition[END_REF], in-processing [START_REF] Hu Zhang | Mitigating unwanted biases with adversarial learning[END_REF][START_REF] Agarwal | A reductions approach to fair classification[END_REF][START_REF] Iosifidis | Adafair: Cumulative fairness adaptive boosting[END_REF], post-processing [START_REF] Hardt | Equality of opportunity in supervised learning[END_REF][START_REF] Fish | A confidence-based approach for balancing fairness and accuracy[END_REF][START_REF] Bhargava | Limeout: An ensemble approach to improve process fairness[END_REF], and hybrid-processing 1 [START_REF] Zemel | Learning fair representations[END_REF][START_REF] Iosifidis | Fae: A fairness-aware ensemble framework[END_REF][START_REF] Alves | and Advanced (a) Adult income. (b) COMPAS. (c) LSAC. (d) Taiwanese credit default. Figure 5: Relation between NOCCO (with linear kernel) and fairness measures for all datasets. We compare 5a (Adult income), 5b (COMPAS)[END_REF].

Fairness of decision systems is commonly assessed via metrics that rely on the outcomes2 of decision models, and inspired by anti-discrimination laws under which decision policies or practices can be declared as discriminatory based on their effects on people belonging to certain sensitive demographic groups [START_REF] Grgic-Hlaca | The case for process fairness in learning: Feature selection for fair decision making[END_REF], e.g., gender, race, age and sexual orientation. Still, fairness is a social construct [START_REF] Abigail | Measurement and fairness[END_REF] and an ethical concept [START_REF] Tsamados | The ethics of algorithms: key problems and solutions[END_REF], and its definition remains prone to subjectivity.

To overcome subjectivity, several fairness notions have been recently proposed to objectively measure and to capture different aspects of fairness. These include group based notions (e.g., [START_REF] Dwork | Fairness through awareness[END_REF][START_REF] Hardt | Equality of opportunity in supervised learning[END_REF]), individual-based notions (e.g., [START_REF] Dwork | Fairness through awareness[END_REF][START_REF] Biega | Equity of attention: Amortizing individual fairness in rankings[END_REF]), and causal and counterfactual based notions (e.g., [START_REF] Pearl | Causality[END_REF][START_REF] Kusner | Counterfactual fairness[END_REF]). Even though these metrics would replace human subjectivity by objective metrics, they are still dependent on subjective decisions. Indeed, most of the above metrics rely on the choice of particular features (variables) with respect to which fairness will be measured. For instance, in the case of group based metrics, fairness is mostly measured with respect to disparate performance and decision results among different subpopulations. For example, when predicting recidivism risk [START_REF] Angwin | Machine Bias -ProPublica[END_REF], disparate results among the white and black subpopulations, or among men and women. In most of the fairness analysis in ML and DS problems, it is assumed that the population is divided into privileged and unprivileged groups on which disparate results are expected. In other words, the sensitive features that discriminate these groups are known a priori. Examples of commonly considered sensitive features include race, gender, age, sexual orientation, etc. See, e.g., [START_REF] Le Quy | A survey on datasets for fairnessaware machine learning[END_REF] for a set of ML problems and the sensitive features frequently considered in the literature.

However, the choice of sensitive features also remains subjective and, as they are chosen a priori, it is dataset agnostic and does not necessarily entail outcome disparities in a given use case scenario. Indeed, consider an ML problem in which both gender and race are among the set of features. Depending on problem addressed, one may have uneven results with respect to race but similar results regardless the gender. In this situation, one may only be aware of race as a sensitive feature when training the ML model. As further illustration, consider one of the above mentioned platforms for job applications. Ethically, physical disabilities of candidates should not constitute a discriminating (i.e., it should be considered sensitive) of candidates. However, if the job is indeed of physical nature, then some physical abilities may be required and, therefore, this information should not be considered as sensitive.

In this paper, we address this critical issue of identifying sensitive feature or subfeatures (i.e., corresponding protected groups or subgroups), and we propose to defining them as features that entail disparate (unfair) outcomes. However, in order to verify disparities among different (sub)groups, one usually needs to train the ML model and calculate the pertaining fairness metrics. This step may be computationally costly depending on the ML model is adopted. Thus this raises the question: how to detect possible disparities before the training step?

To tackle this issue, we propose a statistical approach based on the Hilbert-Schmidt independence criterion (HSIC) [START_REF] Gretton | Measuring statistical dependence with hilbert-schmidt norms[END_REF], which is commonly used in data science [START_REF] Wang | Learning with hilbert-schmidt independence criterion: A review and new perspectives[END_REF] to measure the dependence between two matrices (or vectors). The HSIC has been used as a feature selection method [START_REF] Song | Supervised feature felection via dependence estimation[END_REF][START_REF] Song | Feature selection via dependence maximization[END_REF], for instance, in classification or regression tasks. The HSIC is usually used under the assumption that "good features" are those that maximize the dependence degree between these features and the vector of labels or values, and one selects the subset of features with the highest dependence degrees, as they are more relevant to the output. In [START_REF] Barshan | Supervised principal component analysis: Visualization, classification and regression on subspaces and submanifolds[END_REF], the authors used the empirical estimate of HSIC in order to derive a supervised version of Principal Component Analysis (PCA). The idea was to perform dimensionality reduction while maximizing the dependence between the projected data and the vector of outcomes. HSIC has also been used when defining the loss-function for regression of classification task [START_REF] Wang | Identification of membrane protein types via multivariate information fusion with hilbert-schmidt independence criterion[END_REF][START_REF] Greenfeld | Robust learning with the hilbert-schmidt independence criterion[END_REF]. In this context, a noteworthy application addresses fairness in automatic decisions by learning models that reduce disparities between sensitive groups [START_REF] Pérez-Suay | Fair kernel learning[END_REF][START_REF] Li | Kernel dependence regularizers and gaussian processes with applications to algorithmic fairness[END_REF].

Our hypothesis in this work is that, for a given feature, if 1. this feature brings information that splits people into different groups, and 2. the HSIC3 between this feature and vector of labels is high, then there will be disparate outcomes among the (sub)groups discriminated by this feature. The intuition is that the features with high HSIC are important to the classification task as they carry key information to discriminate (sub)classes. Hence, it is likely to observe disparate outcomes among the discriminated (sub)groups.

Note that not all features that split the population into two groups can be considered as sensitive. For example, the charge degree in recidivism risk prediction [START_REF] Angwin | Machine Bias -ProPublica[END_REF] may split the population between those who committed a felony or a misdemeanor, but it is not assumed as a sensitive one. Moreover, sensitive features may split the population into more than two groups. For instance, as race, one may have whites, blacks and Asians. For each subfeature (e.g., race whites, race blacks and race Asians), a high HSIC will indicate (i) if the feature is sensitive and (ii) which group division may lead to unfair (disparate) results.

As we will see, HSIC indicates both the sensitive features and the associate groups of individuals that may lead to unfair outcomes. More precisely, high HSIC values between features (or subfeatures, in the case of categorical data) and the label vector may entail disparate results. We attest our hypothesis with empirical evaluations on several datasets frequently used in the literature, and that relate high HSIC values to high group fairness measures. This HSIC based preprocessing approach thus constitutes a noteworthy tool researchers and practitioners to automatically detect those features that should be considered sensitive in an ML or DSs tasks.

It is important to stress that our approach deviates from previous attempts such as [START_REF] Grgic-Hlaca | The case for process fairness in learning: Feature selection for fair decision making[END_REF][START_REF] Grgić-Hlača | Beyond distributive fairness in algorithmic decision making: Feature selection for procedurally fair learning[END_REF] in two major aspects, namely, by seeking an automatic approach to choosing which features should be considered sensitive, and by claiming that the choice of sensitive features should take into account both the task and use case scenario. In particular, we show that certain features often considered as sensitive from ethical and social standpoints do not result in outcome disparaties, and thus their use in ML and DSs should not be prohibited. In fact, our approach also enables proxy detection as we will discuss in Section 5.

The rest of this paper is organized as follows. In Section 2, we outline the notations used in this paper. Section 3 discusses the Hilbert-Schmidt independence criterion and its normalized version. The proposed approach to automatically detect sensitive features and sensitive groups is presented in Section 4. In Section 5, we conduct the numerical experiments and discuss the obtained results. Finally, in Section 6, we present our conclusions and future perspectives.

Main contributions. Here, we simply highlight the main contributions of the paper.

• We address the problem of detecting sensitive features in group fairness settings, and propose a statistical approach based on the HSIC that does not require a trained model to verify disparate outcomes.

• We present preliminary empirical results on four well known datasets that support our hypothesis. The higher is the HSIC-based dependence measure between a feature and the labels, the higher are the disparate results entailed by using the information provided by such a feature.

• The analysis of our results also shows that several features taken as sensitive in the literature do not necessarily entail disparate (unfair) results.

• We also raise some perspectives for future work. Our proposed approach can be easily extended to decide whether a combination features (e.g., the use of features describing gender and race) can lead to disparate results. Moreover, we will investigate how to adapt our framework to multiclass classification problems and to individual fairness settings.

Notation

Throughout this paper, we assume an underlying binary classification problem such that each m-dimensional sample is assigned to a class y ∈ {-1, 1}. The set of n samples and the associated vector of classes are represented by X and y, respectively. Each column of X, defined by X (j) , describes a feature G j , j = 1, . . . , m.

Very often some features in ML problems are categorical. However, as most of the learning algorithms require numbers as input features, we consider in this paper the one-hot encoding method to transform the categorical features into binary ones. Mathematically, a vector X (j) describing a categorical feature G j with feature space {G j,1 , G j,2 , . . . , G j,q } can be replaced by 4 q binary vectors X (j,k) , k = 1, . . . , q, that indicate which category each sample belongs to. We illustrate this scenario with the 3-dimensional dataset

X = male 2 30-60 f emale 0 < 30 f emale 1 > 60 , (1)
whose features are G 1 : gender with possible values male or f emale, G 2 : number ncrimes of committed crimes in the last 5 years, and G 3 : age that can be less than 30, between 30 and 60 or greater than 60 years old.

After encoding the categorical features, one achieves the extended 5-dimensional dataset

X = 1 0 2 0 1 0 0 1 0 1 0 0 0 1 1 0 0 1 , (2)
where G 1,1 = gender male , G 1,2 = gender f emale , G 3,1 = age ≤30 , G 3,2 = age 30-60 and G 3,3 = age ≥60 are the novel subfeatures. Note that, after conducting one-hot encoding for all categorical features, one increases the dataset dimension from m to m.

Hilbert-Schmidt independence criterion

There are several measures adopted in feature selection tasks to evaluate the relation between inputs and outputs [START_REF] Sotiris | Feature selection for machine learning classification problems: A recent overview[END_REF].

As examples, one may cite the correlation [START_REF] Hall | Correlation-based feature selection for machine learning[END_REF] and mutual information [START_REF] Vergara | A review of feature selection methods based on mutual information[END_REF] between variables. Methods based on correlation are generally straightforward to be deployed as this measure is easy to be calculated. However, correlation only provides second-order information between variables and, therefore, does not necessarily imply dependence. On the other hand, mutual information brings the knowledge about the dependence degree between variables. The price to be paid here is that, in order to calculate the mutual information, one needs to estimate the probability density function of such variables. This task is very impractical in most cases.

Another measure that has been used in recent works in the literature is the Hilbert-Schmidt independence criterion [START_REF] Wang | Learning with hilbert-schmidt independence criterion: A review and new perspectives[END_REF]. The HSIC measures the dependence degree between finite number of observations5 (x i , y i ), i = 1, . . . , n. Assume K x and K y as the kernel matrices of x and y, respectively. In our analysis, we consider the radial basis function (RBF) kernel, defined by

K RBF (z i , z i ) = e -1 n (zi-z i ) 2
, and the linear kernel, defined by K linear (z i , z i ) = z i z i . However, other kernels can also be considered (see [START_REF] Schölkopf | Learning with kernels: Support vector machines, regularization, optimization, and beyond[END_REF] for details). An empirical calculation of HSIC is given by

HSIC(x, y) = tr (K x HK y H) (n -1)

2

,

where H = I -n -1 ee T is the centering matrix and e is a n-dimensional column vector of ones. Note that, even if the dataset x is centered, this is not necessarily true for the kernel K x . In order to centralize this kernel, one applies the centering matrix H on both sides of K x . Recall that

tr (K x HK y H) = tr (HK x HK y ) = tr Kx K y ,
where Kx is the centered kernel. Moreover, when we multiply K x on the left (resp. right) by H, one removes its columns (resp. rows) mean.

As highlighted in [START_REF] Li | Kernel dependence regularizers and gaussian processes with applications to algorithmic fairness[END_REF], the HSIC calculation is sensitive to the scale of the observations and, therefore, an appropriate normalization should be conducted in order to compare relative dependence degrees. We thus consider a normalized version of HSIC called NOCCO (NOrmalized Cross-Covariance Operator) [START_REF] Fukumizu | Kernel measures of conditional dependence[END_REF]. It is defined as follows:

N OCCO(x, y) = tr (R x R y ) , where • R x = HK x H (HK x H + n I n ) -1 , • R y = HK y H (HK y H + n I n ) -1 ,
• is a regularization parameter (e.g., 10 -6 ), and

• I n is a n × n identity matrix.
It is important to highlight that the use of NOCCO has two main advantages. Firstly, as the calculation is based on the trace of matrix products, it is easily calculated (we do not need to estimate, for instance, probability density functions). Secondly, the use of kernels brings more information than second-order moment to measure the relation between variables. Therefore, it is useful to estimate the dependence degree between them.

Proposed approach

The purpose of this work is to automatically detect sensitive features and the associated impacted groups. Our hypothesis is that if a feature that provides information that splits groups of people (e.g., if the individual is white or black) is important to increase the performance of a machine learning model (measure by means of NOCCO), then this feature may create disparities between such groups in terms of performance measures.

Consider the extended dataset X, with columns X(1) , . . . , X(m) , obtained after encoding categorical features. As mentioned in Section 3, HSIC can be used to measure the dependence between two vectors. For each numerical feature and encoded features (i.e., for all columns of X), we calculate NOCCO with respect to the vector of outcomes y. For numerical features, the dependence degree with respect to the vector of labels is the NOCCO itself. In the case of categorical features, if this feature is a binary one, then the NOCCO for both groups will be the same and, therefore, we only need to calculate itonce. Otherwise, if the categorical features have three or more categories, we assume as the dependence degree the maximum NOCCO among the associated encoded features. This procedure leads to m measures of dependence

d j = tr R Gj R y , j = 1, . . . , m, (3)
where R Gj = HK Gj H HK Gj H + n I n -1 and K Gj is the kernel matrix associated with feature G j . In the example of Section 2, if NOCCO for age ¡30 is greater than the NOCCO for both age 30-60 and age ¿60, we define d 3 as the dependence measure between X(4) and y.

Based on all d j , j = 1, . . . , m, we may evaluate whether a feature should be considered as a sensitive one. In this paper, we assume an automatic procedure that takes into account a predefined threshold t. Algorithm 1 presents a pseudo-code of our proposal. For each feature G j , if d j ≥ t, G j is considered sensitive and the G j,k that leads to the higher NOCCO is the harmed group. Otherwise, we do not consider G j as sensitive. In our empirical setting, we take t to be the median value among all d 1 , . . . , d m . Note that this can be easily adapted to integrate other values of t.

Algorithm 1 (Automatic approach for sensitive features detection)

Input: X and y.

Output: Set of sensitive features S = {G j , G j , . . .} and sensitive groups S g = {G j ,k , G j ,k , . . .}. 1: Encoding categorical features: Encode the categorical features and obtain the extended dataset X ← encode (X).

2: Calculate the kernel matrix of y: K y = kernel (y) 3: Calculate the NOCCO for all features and subfeatures: for j ∈ {1, . . . , m} do if G j is either a numerical or binary feature then Calculate the kernel matrix associated with G j : K Gj = kernel (G j ) Calculate the dependence measure:

d j = tr R Gj R y else
Calculate the maximum NOCCO for the categorical variable: for k ∈ {1, . . . , q} do Calculate the kernel matrix associated with G j,k :

K G j,k = kernel (G j,k )
Calculate the dependence measure: d k j = tr R G j,k R y end for Define the dependence measure for feature G j : d j ← max d 1 j , . . . , d q j end if end for 4: Create the set of sensitive features and sensitive groups: S = ∅ and S g = ∅ 5: Calculate the median value among all d j , j = 1, . . . , m: t = median (d 1 , . . . , d m ) for j ∈ {1, . . . , m} do if G j is candidate for sensitive feature and d j ≥ t: then Update the set of sensitive features: S ← {S, G j } Update the set of sensitive groups: S g ← {S g , G j,k * } such that d k * j = max d 1 j , . . . , d q j end if end for

Empirical Evaluation

In this section, we present the experimental setup to evaluate our proposal on several real datasets. We identify possible sensitive features as well as the associated sensitive groups. Moreover, we provide a comparison between the obtained NOCCO values and the fairness measures with respect to those features with high NOCCO values.

Datasets

In order to assess our proposal, we consider four datasets frequently used in the literature: Adult income 6 , COMPAS recidivism risk [START_REF] Angwin | Machine Bias -ProPublica[END_REF], Law School Admission Council (LSAC) [START_REF] Wightman | LSAC national longitudinal bar passage study[END_REF] and Taiwanese credit default [START_REF] Yeh | The comparisons of data mining techniques for the predictive accuracy of probability of default of credit card clients[END_REF]. A brief description and how we preprocessed some features are provided in the sequel. See [START_REF] Le Quy | A survey on datasets for fairnessaware machine learning[END_REF] for further details as well as the list of features frequently considered as sensitive ones in the literature.

• Adult income dataset: In this dataset, the aim is to predict whether a person makes over 50K a year based on the following features: age, workclass, educational-num (educational degree), marital-status, occupation, relationship (husband, not in family, other relative, own child, unmarried or wife), race (Indian-Eskimo, Asian-Pacific Islander, black, white or other), gender (male or female), capital-gain, capital-loss, hours-per-week and native-country. After removing samples with missing values, we achieved 45222 samples. We also rearranged some categorical features: age = {< 25, 25-60, > 60}, workclass = {private, non-private}, maritalstatus = {married, never-married, other}, occupation = {office, heavy-work, service, other} and nativecountry = {US, non-US}. One typically assumes age, race and gender as sensitive features.

• COMPAS recidivism risk dataset: This dataset has the goal of classifying individuals as a potential criminal recidivist. There are 6167 samples and 8 features, namely sex (male or female), age cat (age category -less than 25, between 25 and 45 or greater than 45), race, juv fel count (number of juvenile felony), juv misd count (number of juvenile misdemeanor), juv other count (number of others infractions), priors count (number of priors) and c charge degree (charge degree -felony or misdemeanor). We here rearranged race {Caucasian, African-American, Other}. Race and sex are assumed as sensitive features.

• Law School Admission Council (LSAC) dataset: In this dataset, there are 23726 candidates described by 11 features: decile1b (decile given the grades in the 1st year), decile3 (decile given the grades in the 3rd year), lsat (score), ugpa (undergraduate GPA), zfygpa (1st year law school GPA), zgpa (cumulative law school GPA), fulltime (full-time or part-time work), fam inc (family income), male (whether the student is male or female), race and tier (tier of the law school). The goal is to predict whether a candidate would pass the bar exam. One assumes gender and race as sensitive features.

• Taiwanese credit default dataset: In this dataset, the aim is to predict customers' default payments in a Taiwanese institution. The 23 features are the following: limit bal (amount of given credit), sex (male or female), education (graduate school, university, high school or others), marriage (married, single or others), age (less than 35 or at least 35), repayment status (six features, from April to September), amount of bill statement (six features, from April to September) and amount paid (six features, from April to September). We assumed sex, education and marriage as sensitive features.

Fairness metrics

In order to attest our hypothesis that the NOCCO value can be associated with fairness measures (and, hence, used as a preprocessing approach to detect sensitive features), we compare the obtained values of fairness measures with respect to features with high NOCCO values. These measures are based on the following performances extracted from the trained ML model:

• True positive (T P ): # of instances correctly classified as class 1.

• True negative (T N ): # of instances correctly classified as class -1.

• False positive (F P ): # of instances wrongly classified as class 1.

• False positive (F N ): # of instances wrongly classified as class -1.

Aiming at evaluating some fairness measures when splitting the dataset into different groups of people, we restrict the aforementioned performance measures to such groups. Suppose, for example, that feature G j contains information that splits people into groups {G j,1 , . . . , G j,q }. As the NOCCO for group G j,k will indicate the dependence degree between such a group and the vector of labels, when evaluating fairness, we consider the disparities when people is divided between those belonging to group G j,k and those that does not. Therefore, we restrict the performance measures such that F P G j,k means the number of instances belonging to group G j,k wrongly classified as class 1 and T N -G j,k is the number of instances not belonging to group G j,k correctly classified as class -1 (we referred to -G j,k as the instances not belonging to group G j,k ).

Given the performance measures conditioned on groups, we may define the fairness measures considered in our experiments:

• Predictive equality (P E): A classifier satisfies predictive equality if both groups have equal false positive rate. In other words, we should have low vales of

f P E = F P G j,k F P G j,k + T N G j,k - F P -G j,k F P -G j,k + T N -G j,k
.

• Equal opportunity (EP ): A classifier satisfies equal opportunity if both groups have equal true positive rate.

In other words, we should have low vales of

f EP = T P G j,k T P G j,k + F N G j,k - T P -G j,k T P -G j,k + F N -G j,k .
• Equalized odds (EO): A classifier satisfies equalized odds if both groups have equal true positive and false positive rates. In other words, we should have low vales of

f EO = T P G j,k T P G j,k + F N G j,k - T P -G j,k T P -G j,k + F N -G j,k + F P G j,k F P G j,k + T N G j,k - F P -G j,k F P -G j,k + T N -G j,k
.

• Overall accuracy equality (OAE): A classifier satisfies overall accuracy equality if both groups have equal correct classification for both classes -1 and 1. In other words, we should have low vales of

f OAE = T P G j,k + T N G j,k n G j,k - T P -G j,k + T N -G j,k n -G j,k ,
where n G j,k and n -G j,k are the number of instances that belong and do not belong to group G j,k .

Result analysis

In our experiments, we first calculated the NOCCO for each feature and encoded feature. Thereafter, we calculate the median NOCCO and define the features that could be considered as sensitive ones. Aiming at attesting the usefulness of NOCCO as a measure to detect sensitive features, we comparing the NOCCO values with a set of fairness measures. For this purpose, we trained a Random Forest classifier using k-fold cross validation (with 10 folds). For each fold and encoded features, we calculated performances and fairness measures (equal opportunity, predictive parity, equalized odds and overall accuracy equality, by taking the absolute difference). For categorical features with more than two groups, we calculated these measures by assuming a one-vs-all strategy.

We present the obtained results for the Adult income, COMPAS, LSAC and Taiwanese credit default datasets in Figures 1,2, 3 and 4, respectively. In NOCCO calculation, we assumed the RBF kernel. Note that, for all datasets and all fairness measures (see Figures 1b,2b, 3b and 4b), the higher is the NOCCO, the higher is the fairness measures (i.e., the higher are the unfair results).

Results on the Adult income dataset In the results with the Adult income dataset, we obtained a median NOCCO equals to 0.0642. Features whose NOCCO is greater than this threshold are marital-status, capital-gain, relationship, educational-num, hours-per-week and capital-loss (see Figure 1a). We consider that splitting people based on either education-num or hours-per-week do not create unfair concerns. However, the other two features can be seen as sensitive ones. Although age and gender are frequently pointed as sensitive features in this dataset, marital-status and relationship appeared as higher candidates. The higher NOCCO was obtained when splitting people between married and non-married. Moreover, NOCCO is high when we divide the instances in those who are a husband in a relation and those who are not. As in this dataset, husband in a relation is a proxy to know that the person is a man, this information is directly associated with gender. Therefore, not only gender should be consider as a sensitive feature when looking at men/women, but we should also be aware of the information provided by relationship. Another interesting aspect in this dataset is that race led to a very small NOCCO. This suggested that, although frequently considered as a sensitive feature, in this problem, it will not lead to disparate results. We attested this conclusion in the fairness measures presented in Figure 1b. The fairness measures are very small when splitting people by race. Results on the COMPAS dataset In COMPAS dataset, the literature assumes gender and race as sensitive features.

In the NOCCO values presented in Figure 2a and with t = 0.0390, only race (when dividing people between African-American and other ones) and age category (greater than 45) could be considered as sensitive features. Although sex is frequently considered as sensitive, the obtained NOCCO value indicates that, in this dataset, it does not entail disparities between male and female. See in Figure 2b that fairness measures when splitting people by sex are much lower in comparison with race or age category. Another interesting remark presented in Figure 2b is that the the relation between NOCCO and most of the fairness measures are nearly linear.

Results on the LSAC dataset

The NOCCO values for the LSAC dataset are presented in Figure 3a. Only race and sex is considered as candidates for sensitive features in this dataset. With a median value of 0.0635, race should be considered as a sensitive feature and male should not. Indeed, as can be attested in Figure 3b, while splitting people by race creates unfair results (we obtained high values of fairness measures), a division by gender does not (the fairness measures are very small).

Results the Taiwanese credit default dataset Another interesting result was obtained with the Taiwanese credit default dataset. One may see in Figure 4a that features generally considered as sensitive ones (sex, education and marriage) have very small NOCCO values. Features associated with payment status have a much higher dependence with class labels than the considered sensitive features. Therefore, with a median NOCCO value of 0.2789, no feature could be defined as sensitive. This result is attested by the small values (with most of them lower than 0.10) of the fairness metrics presented in Figure 4. Indeed, if one compares these fairness metrics with the ones achieved with the previous datasets, one may assume that they are small.

Consistency with linear kernel

In order to verify the consistency of the previous results by assuming another kernel, we evaluate the relation between linear kernal-based NOCCO and the fairness measures. The results for the Adult income, COMPAS, LSAC and Taiwanese credit default datasets are presented in Figures 5a,5b, 5c and 5d, respectively. Clearly, all the results with the linear kernel are consistent with the ones obtained by using the RBF kernel. The only difference can be seen in Adult income dataset with respect to the detected sensitive group for race. Instead of race White (as detected by using the RBF kernel), we here assigned race Black as the sensitive group. The reason for this disparity lies in the number of samples for each race category. As the number of samples whose categories are different from white or black is less than 5%, when using one-vs-all strategy, the results for race White and race Black are practically the same. Therefore, we can safely consider that there was not a disparity between the obtained results.

Conclusions and perspectives

In this paper, we proposed an automatic approach to detect sensitive features based on a normalized version of the Hilbert-Schmidt independence criterion. As our method only requires the dataset and labels, we are able to detect sensitive features without the need of training a machine learning model and evaluating the performance measures.

The numerical experiments on several datasets and different kernels attested that the higher is the dependence between a sensitive feature and the vector of labels, the higher are the outcome disparities with respect to groups discriminated by such a feature. Therefore, the statistical dependence calculated from the dataset and labels proved to be a measure that can guide researchers and practitioners in detecting sensitive features, especially in situations and case scenarios where the pertaining fairness notions seem hard to identify.

We highlight that an automatic procedure to detect sensitive features is of importance in practical problems. Generally, one subjectively defines which features are sensitive. However, we showed that features frequently considered as sensitive in the literature may not entail disparate results. This finding can save effort when developing a machine learning model that takes into account information from predetermined sensitive features in order to mitigate disparate results provided by them.

This work opens interesting perspectives for future works. Although we used the NOCCO to calculate the dependence degree between labels and vectors of features that describe groups in the population (e.g., if the person is black or not), this measure can also be applied to coalitions of features (e.g., subgroups). In this case, we aim at verifying the dependence degree between tow or more sensitive information (e.g., if the person is a black woman with less than 30 years old) and the labels. This may show that combining information provided by two or more features may increase the chance of achieving disparate results. As another perspective, we intend to extend the proposed approach to deal with multiclass classification problems. In such scenarios, for each available class, we may detect which features should be considered as sensitive.

  (a) Dependence measure (NOCCO). (b) Relation between NOCCO and fairness measures.

Figure 1 :

 1 Figure 1: Results for the Adult income dataset (with RBF kernel).

  (a) Dependence measure (NOCCO). (b) Relation between NOCCO and fairness measures.

Figure 2 :

 2 Figure 2: Results for the COMPAS dataset (with RBF kernel).

  (a) Dependence measure (NOCCO). (b) Relation between NOCCO and fairness measures.

Figure 3 :

 3 Figure 3: Results for the LSAC dataset (with RBF kernel).

  (a) Dependence measure (NOCCO). (b) Relation between NOCCO and fairness measures.

Figure 4 :

 4 Figure 4: Results for the Taiwanese credit default dataset (with RBF kernel).

Some works also focus on procedural fairness, i.e., of the decision making processes (means) that lead to the outcomes[START_REF] Grgić-Hlača | Beyond distributive fairness in algorithmic decision making: Feature selection for procedurally fair learning[END_REF]. However, this concern will remain outside the scope of this paper.

In fact, as will be further discussed in this paper, we adopted a normalized version of HSIC.

One generally replaces by q -1 binary features in order to avoid a redundant column. However, in our analysis, the use of either q or q -1 binary features lead to the same results.

We here consider x = {xi} n i=1 and y = {yi} n i=1 as vectors. However, HSIC can be also used to calculate the dependence degree between matrices.

https://archive.ics.uci.edu/ml/datasets/adult

Acknowledgements

Work supported by São Paulo Research Foundation (FAPESP) under the grants #2020/09838-0 (BI0S -Brazilian Institute of Data Science), #2020/10572-5 and #2021/11086-0. The research of the second named author was partially supported by TAILOR, a project funded by EU Horizon 2020 research and innovation programme under GA No 952215, and the Inria Project Lab "Hybrid Approaches for Interpretable AI" (HyAIAI).