Ky Dang Van
email: dangvan@lms.polytechnique.fr.

Modelling of damage induced by contacts between solids

Keywords: Damage, Multi-axial fatigue, Fretting, Moving loads, Shakedown, Finite element analysis Endommagement, Fatigue multiaxiale, Fretting, Adaptation, Éléments finis

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

Introduction

Contact between solids may induce different types of fatigue defects which are very difficult to model. For instance, the prediction of damage and wear induced by rolling-sliding contact is of primary importance for the railway industry. Therefore, one has to evaluate the mechanical state induced by railway traffic generating complex stress and strain distributions. Moreover, in many cases, the occurrence of damage involves plastic strain, and the determination of the limit response under repeated loadings requires special computational methods in order to make the calculations tractable within a reasonable time. Finally rolling contact loading generates complex three dimensional varying stress and strain states which cannot be easily reduced to an equivalent stress or strain for the fatigue prediction. As a consequence, fatigue phenomena induced by rolling sliding contact in rails are still difficult to model and no fully satisfactory approaches are yet available.

Another important example of contact damage is the fretting-fatigue defined as a surface damage induced by small amplitude oscillatory displacements between metal components in contact. Depending on the imposed force or the displacement amplitudes, damage is either wear or crack nucleation. In industrial structures, fretting can be at the origin of catastrophic failures which are particularly difficult to predict. In order to study this phenomenon, different authors, [START_REF] Vingsbo | On fretting maps[END_REF][START_REF] Vincent | Testing methods in fretting fatigue: a critical appraisal, standardization of fretting fatigue test methods and equipment[END_REF] have established test methodologies using different laboratory experimental set up. The results are generally presented as fretting maps, and are useful for the qualitative understanding of damage formation. They show that fatigue crack initiation occurs in the 'mixed stick-slip' regime. Within this regime, relations between the tangential load T , the displacement d and the number of cycles to failure N evolve in a characteristic way: the first loops differ from the stabilised loops in which T and d are proportional. Unfortunately, these parameters are not intrinsic and the obtained results cannot be used for the prediction of fretting on different structures. Another methodology applicable to any kind of structure and loading is needed. For this purpose one must first evaluate the local stress and strain resulting from the cyclic contact loading and second use a multiaxial fatigue criterion in order to predict the crack initiation. For the first task, only numerical calculation is possible. For fretting results from oscillating contact loadings, different elastoplastic computational methods can be used. The incremental method, based on a step by step integration, can be employed to calculate the possible stabilised state, but it imposes cumbersome and lengthy calculations. In order to ovoid this drawback some new methods are proposed.

In conclusion, the general methodology proposed here for designing structures (see Fig. 1) which may undergo damage arising from contact between solids is summarized in the scheme represented above. In order to be able to perform the calculations, the local friction coefficient µ, must be known. However, this parameter is generally not precisely determined. For instance, in classical mechanics µ is often taken between 0.1 and 0.3. This range of values for µ is generally used by railway engineers. However, direct measurements developed recently [START_REF] Fouvry | A multiaxial fatigue analysis of fretting contact taking into account the size effects, in: Fretting-Fatigue Current Technology and Practices[END_REF] lead to much higher values, sometimes as high as 0.9 or 1. Moreover, this coefficient may depend of the load and the rolling-sliding regime (traction or braking) as it is shown in some studies [START_REF] Mandel | Résistance au roulement d'un cylindre indéformable sur un massif parfaitement plastique[END_REF]. Evaluation of the local value of this parameter is still a open question. In the following, we shall just consider different possible values of this parameter in order to discuss the influence of the friction coefficient on the different observed types of contact damage.

Evaluation of the stabilised mechanical states due to moving contacts

The prediction of the risk of damage occurrence on a solid submitted to repeated moving contact requires first to evaluate the intrinsic thermo-mechanical parameters and their evolution in the vicinity of the contact zone. As the loading is in general high, plastic flow occurs and the corresponding asymptotic state is either elastic shakedown, plastic shakedown or ratchet. Depending on the load level, the determination of the eventual stabilised state necessitates sometimes thousands of cycles. Therefore, the use of classical incremental methods is not recommended, and we have proposed for these applications other approaches [START_REF] Van | Steady-state flow in classical elastoplasticity: application to repeated rolling and sliding contact[END_REF].

Stationary methods are efficient for the evaluation of stress and strain in rail: (i) the pass by pass stationary method (PPSM) is suitable for studying the build-up cycle by cycle of the plastic strain and the residual stress distributions; (ii) the direct stationary method (DSM) permits the direct determination of the stabilised state in the case of shakedown (elastic or plastic), and the characterisation of the ratcheting (non-convergence of the algorithm).

The (iii) direct cyclic method is an alternative for studying the same kind of problems as previously mentioned. However, in the case of arbitrary cyclic loadings resulting, for example, from oscillatory contacts, only this last method is applicable.

Another promising method, which will not be discussed in this article, was proposed by M. Peigney and C. Stolz [START_REF] Peigney | An optimal control approach of inelastic structures under cyclic loading[END_REF]. It is based on optimal control techniques.

Recall of the stationary method

This method was already presented in previous papers [START_REF] Van | Steady-state flow in classical elastoplasticity: application to repeated rolling and sliding contact[END_REF][START_REF] Van | On some recent trends in modelling of contact fatigue and wear in rail[END_REF]. It applies to a structure (a rail in the presented applications) subjected to contact loads moving with a velocity V in a fixed direction e z (the longitudinal axis of the rail). Assuming a steady state in a reference frame moving with the loads, the time derivative of any material quantities becomes:

Ḃ =-B,z (1)
Writing the equations governing the problem in the moving frame related to the load, one obtains modified balance equations:

div σ = ρV 2 u ,zz
and constitutive laws:

σ = Lε e + σ 0 A k = Z : α k ε = ε e + ε p -ε p ,z = Λ ∂f ∂σ ,f ( σ , A k ) 0,Λ 0,Λ f( σ ,A k ) = 0 α k,z = Λ ∂f ∂A k
where Λ denotes the plastic multiplier, f the yield function, σ 0 the initial stress, α k are internal parameters and A k their associated forces linked by the tensor Z, L is the fourth order elastic coefficient tensor.

As described in Eq. (1), time derivatives have been replaced by space derivatives so that the integration of the previous set of equations is performed in space, along the direction of motion of the loads. In this method, neither the load nor the structure have to be translated.

Two numerical procedures for the determination of the stabilised state are derived. The first one is the pass by pass stationary method (PPSM) for the calculation of a single pass. Repeating the procedure for the successive passes, the method computes the stabilised state and its characteristics (number of cycles before reaching the stabilised state, residual stresses, plastic deformations, other internal parameters...). The stabilised state can be ratcheting, and in such a case, the ratcheting rate corresponds to surface displacement or plastic strains caused by a pass.

In the second procedure called 'Direct stationary method' (DSM), the possible stabilised state due to repeated moving load is searched directly. The stress field should be periodic, the plastic strain field and the field of other internal parameters also in the absence of ratcheting. The ratcheting is indicated by a non-convergence of the method.

These methods permit computation of strains and stresses in the rail due to repeated moving contacts.

Recall of the 'Direct cyclic method'

This method proves its efficiency in the case of arbitrary cyclic loading resulting from small oscillatory contacts. The main idea is to take advantage of the periodicity of the loading. The stabilised mechanical response is either a shakedown state or is ratcheting. Knowledge of the exact evolution of the cycles is of limited interest. The direct cyclic method (DSM) allows the determination of the asymptotic response of the structure without following the loading step by step.

Let us denotes by F d (t), t ∈[0,T], the external loading path of period T corresponding to either fixed or moving contact. The interval [0,T] is discretized in n intervals 0 <t 1 < •••<t n = T . In order to seek directly the asymptotic elastoplasticity solution, we proposed the following iterative procedure. For each internal elastoplastic iteration j at time t i , we first seek the stress tensor σ js i which is statically admissible with the loading F d (t i ) and the strain tensor ε j i , which is cinematically admissible with the prescribed displacement U d (t i ) corresponding to the instant t i . These two conditions are equivalent to the resolution of equilibrium equations:

KU j i = Q i + Q p j -1 i in which Q i , Q p j -1 i
, and K represent the nodal external forces corresponding to F d (t i ) the vector of loads due to plastic strain ε p j -1 i (assumed known) and the stiffness matrix of the linear elastic structure respectively.

In a second stage, we determine the stress tensor σ j p i+1 plastically admissible. For that, at each time, the plastic strain ε j p i+1 is obtained by projecting the strain increment ε j i+1ε j i on the elastic convex domain defined by the initial state (σ

j p i ,ε j i ,ε p j i ) at time t i which is known. σ j p
i+1 is then calculated from the constitutive law written in incremental form:

ε j i+1 -ε j i = K σ j p i+1 -σ j p i + ε p j i+1 -ε p j i
The stabilised state is obtained when the plastic strain ε p j 1 at the beginning of the cycle is the same as at the end of the cycle ε p j n . Otherwise, the same steps are performed in the iteration j + 1 by initialising ε

p j +1 1 ε p j +1 1 to ε p j n ε p j
n and the process is repeated until elastic or plastic shakedown is achieved. Ratcheting is revealed by the non-convergence of the procedure.

Comparison between this method and incremental method presented in [START_REF] Maouche | On a new method of evaluation of the inelastic state due to moving contacts[END_REF], proves the validity and the efficiency of the DCM method.

Fatigue treatment

Prediction of damages induced by moving contact loadings (rolling contact, or structures undergoing fretting ...) necessitates the use of multi-axial fatigue criteria able to deal with complex multi-axial situations. In fact, the local stress and strain fields are completely tri-dimensional with no fixed directions during their evolution. It is important to describe in a fatigue criterion the main features of the loading paths at any point of the structure. For this purpose a unified fatigue approach was proposed based on multi-scale approach and shakedown assumption.This new point of view on the fatigue strength and the fatigue damage gives a general framework to fatigue modelling, for high cycle as well as for low cycle. It allows analyzing, in an efficient manner, all types of industrial fatigue problems, including contact.

It is well known that the fatigue phenomenon in metals is associated with the movement of dislocations and creation of persistent slip bands in metal grains. Conducting classical laboratory tests the following features can be observed: in the high cycle fatigue (HCF) regime, the stress amplitude is low, the test specimen remains globally elastic and plasticity is localized in some well oriented grains; in the low cycle fatigue (LCF) regime, plasticity occurs at the structural level as well as the grain level. However, the fatigue mechanisms are the same in both cases and the different fatigue regimes differ essentially in the spatial extension of the plasticity. In order to propose a unified view for fatigue modelling, we shall suppose that the fatigue initiation corresponds to the limit of shakedown possibility of the material in both scales. In what follows we shall distinguish the mesoscopic scale corresponding to the grain and the macroscopic scale used by engineers which supposes implicitly the concept of representative volume element (RVE) defines as sufficiently small volume so that one can distinguish the microscopic heterogeneities and sufficiently large to be representative of the macroscopic behavior. σ ′ , ε ′ , ε ′ P , ρ ′ are respectively the stress strain, plastic strain and residual stress tensors at the mesoscopic scale, and σ , ε, ε P , ρ are the corresponding macroscopic quantities. Engineers calculate the mechanical parameters of a structure not at a point but over a volume that defines the scale which is used. In this volume centered on that point, stress and strain σ and ε are assumed constant. However, this volume contains grains which are heterogeneous. In HCF, at the mesoscopic level, plasticity occurs in certain number of grains and generates heterogeneous plastic strain and local residual stress states distribution in the RVE. Since the first cracks initiate in these positions, σ ′ (t) should be evaluated here. This problem corresponds to a localisation problem which cannot be solved without some additional assumptions.

The relation valid at any time t of the loading between macroscopic and mesoscopic stresses is:

σ ′ = L.σ + ρ ′
where L is the elastic localisation tensor and ρ ′ is the local residual stress tensor linearly related to plastic strains. The additive term thus depends of the loading path. In order to calculate σ ′ , we suppose that shakedown occurs before fatigue. More precisely, we assume that fatigue corresponds to the limit elastic shakedown possibility of the structure at the mesoscopic and macroscopic scale. Thank to this assumption, one can apply theoretical results such as the Melan-Koiter static shakedown theorem or its extension to generalised standard materials; see Q.S. Nguyen [START_REF] Nguyen | Stabilité et mécanique non linéaire[END_REF]. (We recall that generalised standard materials are suited for description of material behaviour of metals.) This last result, which gives a sufficient condition for elastic shakedown, can be stated as follows:

"If there exist a time θ and a fixed generalized self equilibrated stress field [ρ(x), A(x)] and a safety coefficient m>1 such that for any point x of the structure and ∀t>θ, f [m(σ el (x, t) + ρ(x)), mA(x)] < 0, the structure will shakedown elastically." f(σ , A) represents the material yield function, A correspond to 'forces' associated to hardening parameters. Thus, (ρ(x), A(x)) are generalized 'residual stresses'. σ el is the stress response of the structure under the same external loading, but under the assumption that the constitutive material has a pure elastic behavior.

This theorem is however difficult to apply, since, the residual stress field ρ must be self equilibrated, a condition which is difficult to fulfill. Another proposal due to Mandel et al. [START_REF] Mandel | Adaptation d'une structure élastoplastique à écrouissage cinématique[END_REF] opens a possibility to overcome this difficulty; the plasticity criterion can be written J 2 [(σα)] <k 2 (P eq ); α the kinematic strain hardening parameter is the center and k( P eq ) the isotropic parameter is the radius of the actual convex of plasticity. Generally, k increases with equivalent plastic strain P eq until some maximum value k * and then decreases. Since σ = σ el + ρ, one obtains:

J 2 (σ el -z) <k 2 (P eq )
where z = αρ. z is not necessarily self-equilibrated since α has not to fulfill this condition.

For cyclic loadings, σ el and ρ are functions of time t ; however, z(t) must be inside or on the convex (hypersphere) with center σ el (t). The trajectory of σ el (t) is known at any point x of the structure because it is the solution of an elastic problem.

Shakedown may occur only if for t>θthe intersection of these convex is not empty. The limit case of possible elastic shakedown occurs when the intersection domain is reduced to a point, for a value corresponding to the maximum of k 2 (p eq ) = k * . The dissipated energy density required for reaching this limit state corresponds to the local bound for possible shakedown. Beyond this bound failure will occur.

We use these different results at different scales in order to propose fatigue modeling. We shall make the following fundamental assumption: fatigue corresponds to the limit of elastic shakedown possibility of the material of the structure at the mesoscopic and macroscopic scale. We then examine the consequences of this hypothesis in order to derive a unified fatigue modeling and the associate computational methods for high cycle and low cycle fatigue.

High cycle fatigue

Characteristics of this fatigue regime are: (i) the loadings are not very important, (ii) the structure is globally elastic; and (iii) fatigue failure initiates only in some grains. It is difficult to evaluate directly the plastic strains and the dissipated energy in the grains that will fail. It is more judicious to proceed indirectly by an inverse way: since the fatigue limit corresponds to the shakedown limit, one can evaluate the local stress at the stabilized state if these grains are at or near the fatigue limit. Starting from the knowledge of σ (x, t) at a point x of the structure, we derive the local stress σ ′ by constructing the smallest hypersphere containing σ (x, t) the center of which corresponds to σ * = α *ρ * ≈-ρ * , because it can be shown that |α * |/|ρ * |≪1.

Thus we have an estimation of the local stress tensor σ ′ (t) = σ (t) + ρ * in grains which may fail by fatigue. Different fatigue criteria can be formulated based on the stabilized stress tensor:

-for instance a linear combination of shear τ and concomitant hydrostatic tension p (Dang Van criterion):

τ + ap < b -Papadopoulos criterion:
k * + ap max <b can be interpreted in the same way: the radius of the elastic domain at the limit state of possible shakedown depends of the maximum hydrostatic tension in the loading cycle.

The different computational steps for the application of these criteria can be found, for instance, in Ref. [START_REF] Van | Fatigue analysis by the multiscale approach, in: High Cycle Metal Fatigue, From Theory to Applications[END_REF].

Low cycle fatigue

In this regime macroscopic plastic strain is important and then difference between ε p and ε ′ p decreases so that there is few difference between macroscopic and mesoscopic parameters. However we propose to choose the dissipated energy W because it is a scalar parameter easy to compute (particularly in presence of thermo-mechanical loadings) and which leads to good predictions on structures. The fatigue criterion can be simply expressed by a relationship, where N is the number of cycles required for crack initiation:

W = cN b
In order to improve the predictivity of the previous relation when high triaxiality is present, it is possible to introduce the hydrostatic pressure in the criterion:

(W + ap max ) = cN b

Application to analysis of some rail damage

Different types of fatigue cracks may appear in the rail head, in the vicinity of the rolling track. These defects are of different natures. In the French railways they depend on the time period. From the 1960s to the early 1970s, the 'kidney shape cracks' are the most frequently encountered. They initiate in the rail head quite far from the rolling surface (around 7 mm) and propagate downwards in the rolling direction. Nowadays, squats and head checks are encountered on straight tracks or shallow curves. Contrary to the kidney shape cracks they initiate very near from the surface (some hundreds of microns to mm) on quite a long distance in the deformed zones of the running track. When they reach a certain size below the surface, they can propagate downwards to form progressively transverse cracks which propagate quite fast. That is the reason why their prevention is necessary. Head checks are angle cracks similar to squats which appear on the high rail in curves and crossing rails. The reasons of this difference between the two types of observed damage (in the depth and near the surface) are, on the one hand the nature of the rail steel (in the 1960s, the rail heads contained many inclusions due to the use of the Thomas process steels) and, on the other hand, the increasing of load and power per axle. The snake skin defect is typically found in some parts of the RER network (Réseau Express Régional) of Paris. It occurs near some railway stations with very high density and heavy traffic in the acceleration and deceleration region of the rails. These defects correspond to surface micro-crackings, approximatively parallel to the rolling-sliding direction.

To understand the circumstance of occurrence of these defects, studies were undertaken using the proposed methodology:

-First we compute strains and stresses in the rail due to repeated moving contacts, using the stationary algorithms presented previously, in order to determine the nature of the stabilised state in relation to the traffic and the type of rail steel (elastic shakedown, plastic shakedown or ratcheting). -Then we use the multiaxial high cycle fatigue criterion to predict the fatigue behaviour of rails.

Kidney shape cracking

The 2-D calculations were performed to evaluate the stabilized state and the residual stress pattern in the rail symmetry surface at different depths. In this way, we obtained the local stress cycle for different traffic conditions. Application of the fatigue criterion allows us to predict the locus of the crack initiation in relation with the load parameters. For instance, in the case of a 2D Hertzian contact and full sliding, the results obtained are summarized on Fig. 2. This type of fatigue corresponds to the region 1 of the map: initiation in depth, high cycle fatigue in the elastic and elastic shakedown zones.

Region 2 where friction coefficient becomes important corresponds to surface damage. The method we propose is applicable for any kind of rolling, with and without sliding regime.

Analysis of squat and head check

A complete 3-D rolling contact simulation was performed to simulate the phenomena. In this case, the microscopic observations of the orientations of the fibers in the running zone show that the superficial contact effects have to be taken into account. In these observations one can distinguish different fibers orientations (corresponding to severe local superficial plastic flow) depending on the position on the rail head surface and squat damage appears in a precise place. Dynamic contact on wheel-rail system calculations performed by INRETS, based on Kalker method, shows complex evolution of the rail-wheel contact with different shapes of the contact area, pressure and spin distribution. To simulate the phenomena, we choose to reproduce a simplified sequence, consistent with the real cinematic. The PPSM method was used to make calculations for the evaluation of the mechanical state in the rail. Finally, fatigue analysis was conducted using the multiaxial fatigue criterion, thanks to a specific analysis software. The constants a and b of the criterion are determined using fatigue limit of the rail steel determined by classical fatigue tests (rotating-bending f -1 = 450 MPa, alternate torsion t -1 = 270 MPa). The stresses are calculated at the stabilized state reached after the contact sequence mentioned previously. On Fig. 3 the loading path of three different points (gauge side, central zone, field side) of the running surface is shown. It can be seen that the fatigue crack are likely to occur in the central zone.

Snake-skin defects

The same type of approach was employed. First, the use of the previous dynamic rail-wheel software gives the conditions of contact, then we perform the repeated elastoplastic calculations to determine the stabilized mechanical state, and finally we apply the fatigue multiaxial criterion. The obtained predictions are shown on Fig. 4. We can see that the high friction coefficient increases the danger of snake-skin defect occurrence.

Application to fretting fatigue

Our objective is to propose an intrinsic methodology to predict fretting which can be applied to any kind of structure. For this purpose, the approach is based on numerical calculations in order to compute the mechanical quantities resulting from the cyclic contact loading. The obtained stabilized stress variation is then introduced in the multiaxial fatigue criterion to check if crack may appear. The first verification of the validity of the proposed method consists in the interpretation the results of a laboratory fretting test.

The following experimental set-up was used by C. Petiot in order to study this phenomenon [START_REF] Petiot | An analysis of fretting-fatigue failure combined with numerical calculations to predict crack nucleation[END_REF]: it consists of two cylindrical pads clamped by a clamping force P again the two sides of a flat uniaxial fatigue specimen tested in repeated tension under constant amplitude loading (0,S max ). The prescribed oscillatory displacement between the pads is linked to the prescribed oscillatory fatigue stress S(t) in the specimen. A special device permits one to measure simultaneously the clamping force and the tangential force T(t) due to friction and related to the accommodation displacements. The variations of T(t) are recorded for each fatigue cycle and plotted as function of the fretting fatigue stress S(t) (called hereafter fretting-fatigue loop). Due to the contact of the pads, the fatigue limit of the specimen decreases in an important manner compared to the fatigue limit without these two pads. In fact, by varying the operating parameters (P , S max ), three regimes are observed: the stick regime: fretting-fatigue loops are quite linear during the test which means that the displacement between the contacting surfaces is mainly accommodated by elastic deformation. No damage appears during the 10 7 cycles of the test; the mixed stick-slip regime: The loops start with an closed elliptical shape and stabilized linearly in the same way as in the previous regime. There is partial slip which decreases with the number of cycles. In this regime crack nucleation is observed at the edges of the contact; the gross slip regime: the loops present a trapezoidal shape corresponding to full slip; in this regime particle detachment is observed corresponding to severe wear.

The question is the following: is it possible to predict the initiation (and no initiation) of cracks in relation with the loading parameters (P , S max , R, radius of the pads, ...) and the intrinsic material fatigue properties (material Wöhler curve) determined by classical fatigue tests without contact?

In order to calculate the asymptotic stress state near the contacts, C. Petiot et al. [START_REF] Petiot | An analysis of fretting-fatigue failure combined with numerical calculations to predict crack nucleation[END_REF] used the incremental method, whereas N. Maouche et al. [START_REF] Maouche | On a new method of evaluation of the inelastic state due to moving contacts[END_REF] employed the direct cyclic method which is more convenient to use. These stress cycles are then introduced in the Dang Van multiaxial fatigue criterion in order to check the danger of fatigue. The most critical points are located at the surface at the edge of the contact. The most critical loading path (τ,p) for each regime resulting from the stress cycles in the stabilized state is plotted in the (τ,p) diagram shown on Fig. 5. In the sticking regime, the loading path (τ,p) is beneath the fatigue material line, so no damage occurs. In the stick-slip regime, the loading path intersects the material fatigue line, concluding with crack failure. The prediction of crack initiation is thus in good agreement with the experimental observations.

Prediction of fretting on a turbine disk [14]

Aero-engines incorporate a large number of mechanical joints where fretting-fatigue is an important life-controlling factor. Safety requirements of aeronautic components make prediction and avoidance of fretting-fatigue a critical reliability issue. Blade disk attachments are typical examples. They have to endure high centrifugal forces leading to elevated contact and bulk stress, combined with small relative displacements due to vibrations in operation, make an important, if not the main, design criterion. In order to explore a life method to predict and quantify an industrial fretting fatigue case, engineers of MTU aero engines company studied with the help of researchers from the École polytechnique the turbine spin-pit test performed on a turbine disk configuration [START_REF] Arrieta | Modelling attempts to predict fretting-fatigue life on turbine components[END_REF]. The disk is fitted with 47 blade 'dummies' that are connected to the disk by a fir tree-type attachment geometry. A two dimensional plane strain computational model was employed to derive the stress and the strain fields in the contact zone of the disk and the dummy. The calculations were made by incremental methods for ten loading cycles to converge to a stable shakedown state. To evaluate the effect of friction coefficient on the fretting fatigue life, several computations were performed with varying friction coefficient µ varying from 0 to 0.9. Fig. 6 shows the stress distribution on the outer contact lobe of the attachment during the 10th loading cycle for a friction coefficient of 0.3. The highest stress occurs on the trailing edge of the contact on the disk side. Fig. 7 shows the Dang Van loading paths on the critical disk locations for µ of 0, 0.3, 0.9 at one of the critical locations (inner disk lobe). From the lobe, it is clear that for no friction no crack initiation is expected, as the loading path does not cross the material line at any instant during the loading cycle. On the contrary, for higher friction coefficients (µ = 0.3, 0.9), the loading path crosses the fatigue line during the loading phase. The corresponding critical locations determined by the fatigue criterion lie on the trailing edges on the disk surface. This result matches very well with experimental observations.

Fig. 1 .

 1 Fig. 1. General methodology for designing structures.

Fig. 2 .

 2 Fig. 2. Numerical shakedown and fatigue map for line contact and full sliding.

Fig. 3 .

 3 Fig. 3. Squat modelling: loading path at gauge side, field side and central zone where crack initiates.

Fig. 4 .

 4 Fig. 4. Danger of 'snake skin' occurrence for 3 values of friction coefficients 0.1, 0.7, 1.5.

Fig. 5 .

 5 Fig. 5. Loading path in Dang Van's diagram in stick regime (no crack initiation) and mixed-stick regime (crack initiation because the loading path cross the fatigue criterion).

Fig. 6 .

 6 Fig. 6. Distribution of equivalent stress near the contact of turbine-disk configuration at the 10th loading cycle for µ = 0.3.

Fig. 7 .

 7 Fig. 7. Loading path in Dang Van diagram for critical locations on the disk surface for µ = 0 (no crack initiation), µ = 0.3, 0.9 (crack initiation).


Mariko Tatsumi

Silviu-Ioan Filip

Caroline White

Olivier Sentieys

Guy Lemieux

Mixing Low-Precision Formats in Multiply-Accumulate Units for DNN Training

The most compute-intensive stage of deep neural network (DNN) training is matrix multiplication where the multiply-accumulate (MAC) operator is key. To reduce training costs, we consider using low-precision arithmetic for MAC operations. While low-precision training has been investigated in prior work, the focus has been on reducing the number of bits in weights or activations without compromising accuracy. In contrast, the focus in this paper is on implementation details beyond weight or activation width that affect area and accuracy. In particular, we investigate the impact of fixed-versus floating-point representations, multiplier rounding, and floatingpoint exceptional value support. Results suggest that (1) lowprecision floating-point is more area-effective than fixed-point for multiplication, (2) standard IEEE-754 rules for subnormals, NaNs, and intermediate rounding serve little to no value in terms of accuracy but contribute significantly to area, (3) lowprecision MACs require an adaptive loss-scaling step during training to compensate for limited representation range, and (4) fixed-point is more area-effective for accumulation, but the cost of format conversion and downstream logic can swamp the savings. Finally, we note that future work should investigate accumulation structures beyond the MAC level to achieve further gains.

I. INTRODUCTION

With the increasing adoption of deep learning (DL), particularly with embedded/edge computing scenarios and the increasing complexity of state-of-the-art DNNs, there is great interest in improving the cost/performance ratio of DNN training. To make training more affordable, one has to inevitably use low-precision arithmetic (e.g. 8-bit and lower formats), as opposed to the current de facto standard of 32-bit IEEE-754 single-precision floating-point (FP32). Smaller data sizes can reduce the hardware needed and can potentially improve memory footprint, bandwidth, and power.

While low-precision DNN training has been investigated by prior art, few works directly examine hardware implementation complexity. Accordingly, the precise hardware resources saved by replacing FP32 with a low-precision data type is left ambiguous. However, reducing the resource demand per MAC unit, while maintaining accuracy, is crucial to deploying more computational power within a fixed hardware budget.

In order to address these ambiguities we apply different low-precision data formats to the MAC units used within a general matrix multiply (GEMM) accelerator for DNN training and investigate the impact on training accuracy and hardware resource utilization in an FPGA context.

In some previous work, the exact methodology used for implementing low-precision calculations is not described. To do so accurately, one would need to avoid standard accelerated libraries and emulate each low-precision calculation, requiring greater effort to code and leading to a significant slowdown of the training process. One quick technique, used for instance by QPyTorch [START_REF] Zhang | QPyTorch: A Low-Precision Arithmetic Simulation Framework[END_REF], applies a quantization layer after each composite computational layer like convolution or GEMM. In such a scheme, all intermediate calculations use a higher precision, and only the final results are 'downconverted' to the reduced precision. Such techniques do not properly model precision loss, leading to overly optimistic results with inflated training accuracy. In this paper, we accurately model the arithmetic of each individual operator and storage unit within a customwritten GEMM kernel that has two implementations (FPGA and GPU) which are cross-validated for bit-level accuracy.

Consider, for example, the validation accuracy of training LeNet5+MNIST shown in Fig 1, where all multiplier inputs are quantized to 5 exponent bits and 1 or 2 mantissa bits (E5M1 and E5M2, respectively). With FP32 accumulation followed by quantization (blue/orange curves), both E5Mm formats look viable, but with low-precision accumulation in the E5M2 format (green curve) accuracy is slightly degraded. Worse yet, the E5M1 format (red curve) result is catastrophic and never converges. This shows the importance of accurately modeling all stages of low-precision arithmetic.

The main contributions/observations are as follows:

acceleration to study the resource-accuracy tradeoffs of using custom data types and operators; • relaxing IEEE-754 rules (i.e., removing subnormal support and repurposing exceptional value encodings) for multiplication saves area and maintains accuracy; • individual operator area is minimized using different number representations (i.e., a floating-point multiplier and a fixed-point accumulator), but format conversion and additional downstream logic (due to the wider output) quickly overcomes any apparent savings; and • in a system like ours (where no effort is made by hardware or software to limit the number of terms that appear in an accumulation chain), the size of a fixed-point accumulator is more sensitive to the target model/dataset than a floating-point multiplier/accumulator.

II. BACKGROUND AND RELATED WORK

We first review DNN training together with a technique called loss scaling that is used in low-precision training to maintain representation range and accuracy.

A. DNN Training

DNNs are a sequence of layers which usually take the form

f (x) = g(W T x + b) (1)
where x ∈ R m is a layer input vector, f (x) ∈ R n is the layer output, W ∈ R m×n and b ∈ R n are trainable parameters (weight and bias terms), and g is a non-linear activation function applied element-wise. The application of Eq. ( 1) on all layers in the network is called inference. Equation ( 1) is at the core of most layers, convolutions included. A common approach is to cast a two-dimensional convolution operation into a matrix multiplication operation coupled with so-called im2col and col2im transforms [START_REF] Jia | Learning Semantic Image Representations at a Large Scale[END_REF].

DNN training consists of finding a set of parameters that (approximately) minimizes an empirical loss function

L(W) = 1 B B i=1 ℓ i (W) (2)
where W is the complete set of trainable parameters in the network, each ℓ i (W) is an averaged loss function on a batch of N input data samples, and B is the number of training batches. To simplify notation, we have avoided specifying the dependence of loss functions on input data. A form of stochastic gradient descent (SGD) is generally used to minimize L. At each iteration t, SGD will randomly pick a batch from the training set (i.e., a loss ℓ from {ℓ i } B i=1 ) and update the parameters according to

W (t+1) = W (t) -α t ∇ℓ(W (t) ) (3)
with α t ∈ R + the current learning rate and ∇ℓ(W (t) ) the parameter gradient. This update is applied to every parameter in each layer of the network by using the chain rule. The process is known as back-propagation, with the loss gradient with respect to the output of a layer k called the activation gradient (or delta), which we denote by δ k . If z(x) = W T x, back-propagation through a layer k requires computing the forward path using Eq. ( 1) and the backward path for the delta and weight gradient respectively:

δ k-1 = ∂z ∂x ∂g ∂z δ k and ∇ℓ(W) = ∂z ∂W ∂g ∂z δ k . (4)
Since ∂z/∂x = W and ∂z/∂W = x, both equations in Eq. ( 4) reduce to a matrix multiplication operation with ∂g ∂z δ k . In total, each iteration calculates Eq. ( 3) of batch DNN training with N > 1 batch size usually requires three GEMM operations per layer: two for computing the gradients as just mentioned, and one for computing the layer inputs, which require storing the intermediate results evaluated during the forward path.

B. Loss Scaling

One of the potential problems of using low-precision during training is the limited dynamic range. This is especially important when computing gradient values, which decrease in magnitude as training progresses. Their values might thus fall below the minimum values that low-precision data types can represent, potentially leaving most of their representation range unused. A possible solution [START_REF] Micikevicius | Mixed Precision Training[END_REF] is to scale the loss function by a large value S before the start of the backward pass so that the new (scaled) gradients do not underflow in the low-precision format. The parameter update in Eq. ( 3) is then performed with the gradient divided by S.

Adaptive Loss Scaling: One difficulty with loss scaling is picking a good scaling factor. Choosing an appropriate value is essential in certain cases, but searching for the best value can be time consuming. The approach proposed in [START_REF]Train with Mixed Precision[END_REF] dynamically adjusts the factor as follows: start with a large value of S (e.g. 1024) and if no overflow is detected within a certain number of iterations P , the value of S is doubled; if overflow is detected, the current iteration is skipped and S is halved. The impact on runtime is negligible.

C. Related Work

As DNNs have increased in size, reducing precision below FP32 has been studied for both training and inference.

The current trend in low-precision training with hardware support focuses on 16-bit floating-point formats with half-precision (FP16) [START_REF] Micikevicius | Mixed Precision Training[END_REF] and bfloat16 [START_REF] Kalamkar | A Study of BFLOAT16 for Deep Learning Training[END_REF] alternatives from NVIDIA and Google. More aggressive options using 8-bit [START_REF] Sun | Hybrid 8-bit Floating Point (HFP8) Training and Inference for Deep Neural Networks[END_REF]- [START_REF] Mellempudi | Mixed Precision Training With 8-bit Floating Point[END_REF] and even 4-bit [START_REF] Sun | Ultra-Low Precision 4-bit Training of Deep Neural Networks[END_REF] formats have also been shown to achieve FP32-level results on specific training scenarios. The NVIDIA Hopper GPU architecture offers support for two different FP8 formats [START_REF]NVIDIA H100 Tensor Core GPU Architecture[END_REF]. Since submission, other work has confirmed that training with FP8 formats coupled with scaling factors [START_REF] Micikevicius | FP8 Formats for Deep Learning[END_REF] or various exponent bias values [START_REF] Noune | 8-bit Numerical Formats for Deep Neural Networks[END_REF] can match 32-bit level baselines. Accumulation in these approaches is usually done in 16-bit or 32-bit floating-point arithmetic.

Other approaches have used fixed-point types during training [START_REF] Banner | Scalable Methods for 8-Bit Training of Neural Networks[END_REF]- [START_REF] Wu | Training and Inference with Integers in Deep Neural Networks[END_REF]. For instance, [START_REF] Banner | Scalable Methods for 8-Bit Training of Neural Networks[END_REF] uses INT8/16 for storage and INT32 for accumulation to yield FP32-comparable results. Later work [START_REF] Yang | Training High-Performance and Large-Scale Deep Neural Networks with Full 8-bit Integers[END_REF] used INT16 accumulation during GEMM operations, but prediction accuracy degraded, whereas [START_REF] Fu | FracTrain: Fractionally Squeezing Bit Savings Both Temporally and Spatially for Efficient DNN Training[END_REF] dynamically changed the storage format from INT4 to INT8 during training to achieve FP32-comparable results.

The efficient implementation of quantized DNNs is also widely studied in the FPGA community [START_REF] Luo | RNA: An Accurate Residual Network Accelerator for Quantized and Reconstructed Deep Neural Networks[END_REF]- [START_REF] Xilinx | Higher Performance Neural Networks with Small Floating Point[END_REF]. In particular, [START_REF] Wu | Low Precision Floating Point Arithmetic for High Performance FPGA-based CNN Acceleration[END_REF] and [START_REF] Xilinx | Higher Performance Neural Networks with Small Floating Point[END_REF] combined floating-point multipliers and fixed-point accumulators. They determined that this combination has the potential to achieve a good prediction accuracyresource utilization tradeoff for inference tasks; in present work, we consider using a similar combination for DNN training tasks. On the training side, [START_REF] Fox | Training Deep Neural Networks in Low-Precision with High Accuracy Using FPGAs[END_REF] successfully trained VGG16 and PreResNet-20 networks [START_REF] He | Identity Mappings in Deep Residual Networks[END_REF] on an embedded FPGA using the SWALP [START_REF] Yang | SWALP: Stochastic Weight Averaging in Low Precision Training[END_REF] algorithm by replacing the GEMM multipliers with 8-bit block floating-point operators but retaining 32-bit integer accumulators; the impact of reducing accumulator size is not investigated.

On the simulation side, QPyTorch [START_REF] Zhang | QPyTorch: A Low-Precision Arithmetic Simulation Framework[END_REF] extends the Py-Torch library by adding rounding modules before and after each computation module. This approximates low-precision arithmetic by converting high-precision results into smaller formats. It allows quick testing of ideas, but internal operations within a composite operation (e.g. GEMM) are still done in FP32, making it inaccurate because it neglects overflows and rounding that must occur in low-precision hardware.

Another relevant framework is Barista [START_REF] Vink | Caffe Barista: Brewing Caffe with FPGAs in the Training Loop[END_REF], which extends the Caffe library by adding FPGA acceleration support. One advantage of our work over Barista is that we are using training results measured on real FPGA hardware, whereas Barista is limited to a performance model.

There is also a large body of work on quantization aware training methods (e.g., the Brevitas [START_REF] Pappalardo | Xilinx/Brevitas[END_REF] library for FPGAs), but such methods generally modify the training process to compress models and accelerate inference, not training.

III. DESIGN AND IMPLEMENTATION

To explore the MAC design space, we developed Archimedes-MPO, a DNN training framework with mixedprecision operations. It supports CPUs, GPUs, and FPGAs.

The FPGA-enabled version, shown in Fig. 2, uses PetaLinux, TinyDNN, and an HLS-based GEMM accelerator [START_REF] De Fine Licht | Flexible Communication Avoiding Matrix Multiplication on FPGA with High-Level Synthesis[END_REF]. FPGA support does provide acceleration, but its main purpose is to validate training accuracy and implementation area of low-precision MAC operations within an accelerated system. As a result, it lacks many system-level optimizations expected by end applications.

GPU simulation support, which is bit-accurate with the FPGA version, allows training to scale out so we can explore multiple training scenarios concurrently.

In both GPU and FPGA versions, we implement custom precision GEMM kernels since these are the core operations in both convolution and fully connected layers. All other network components (e.g., activation functions and batch normalization), are performed in software on the host CPU using FP32. We have limited our exploration to GEMMs for expedience and because they amount to the most computation and are the most likely to exhibit precision problems during training. A more in-depth study that also customizes the arithmetic of other operations during training, notably the parameter update in Eq. ( 3), is left for future work.

A. Host Implementation

Archimedes-MPO is a hybrid software/hardware approach that builds upon the TinyDNN [START_REF] Noumi | Tiny-DNN: Header only, Dependency-free Deep Learning Framework in C++14[END_REF] framework. Written in C++14, TinyDNN is a pure software, header-only, dependency-free, templated and optimized DL library capable of multi-threaded SIMD execution on x86 and ARM.

Archimedes-MPO modifies TinyDNN in three key ways: (a) it adds more data types (i.e., fixed-point and floating-point with custom precisions) to extend each layer type with C++ templates and operator overloading; (b) it employs custom lowprecision GPU or FPGA versions for the GEMM kernel; and (c) it implements col2im and im2col transforms needed to offload convolutional layers using the GEMM kernel.

With templates and overloading, Archimedes-MPO can create networks with a different precision and data type for each layer instance. Some layers can specify different precisions for different operations, e.g. multiplier inputs, outputs, and accumulation in GEMM. However, all layer instances that use the FPGA-accelerated GEMM hardware must use the same formats because only one hardware kernel is synthesized.

We explore the use of custom data types within the GEMM kernel and compute all other layers on the host using FP32. Conversion between FP32 and the custom types is done in the accelerated GEMM kernels; this improves compute efficiency, but does not improve memory footprint or bandwidth.

By leaving network operations on the host in FP32, our trained models have FP32 weight values, and that parameter updates are also done in FP32. Reducing the precision during these updates is subject of other research, e.g., [START_REF] Zamirai | Revisiting BFloat16 Training[END_REF]. Nevertheless, our GEMM calculations do use low precision.

B. GPU Implementation

Archimedes-MPO on GPU enables fast batch-mode training on x86/GPU servers running Linux. Since it uses the same TinyDNN core as the FPGA version, it is easy to track changes between the two (e.g. if we modify training strategies). The GEMM kernel uses CUDA to emulate low-precision arithmetic in a way that is bit-accurate with the FPGA accelerator. To decrease the overhead caused by data transfer between layers, all operations (other than the image data loading at the start of each iteration) are implemented as CUDA kernels.

C. FPGA Implementation

The FPGA version of Archimedes-MPO uses a heavily customized version of an open-source GEMM hardware kernel [START_REF] De Fine Licht | Flexible Communication Avoiding Matrix Multiplication on FPGA with High-Level Synthesis[END_REF] written in C using Xilinx Vitis HLS. The Vitis framework also relies upon PetaLinux and the Vitis Embedded Platform (an FPGA logic shell).

The original GEMM-HLS kernel [START_REF] De Fine Licht | Flexible Communication Avoiding Matrix Multiplication on FPGA with High-Level Synthesis[END_REF] supports only standard-width primitive data types (i.e., FP32 and uint8). We modified it to support small custom data types, such as FP8 with different exponent-mantissa breakdowns. In addition, we modified data reading/writing modules in the GEMM kernel to convert data between FP32 and low-precision during transfer. While TinyDNN can use any precision in each layer, FP32 is used on the host side to keep things simple, reduce the overhead of emulating other precisions, and avoid unnecessary data type conversions. During conversion, values outside the representation range are cast to ±∞. Similarly, multiplier and accumulator operators propagate ±∞ when necessary.

We also modified the GEMM-HLS kernel to directly read both transposed (n × m) and non-transposed (m × n) matrix data. This does not affect the resulting clock speed, and kernel throughput remains over 97% of peak GOPS when tested with sufficiently large (2048×2048) transposed and non-transposed matrices. This allowed us to remove the redundant input data transpose during back-propagation. Finally, we also modified GEMM-HLS to remove the use of DSP blocks everywhere outside of the PEs, such as address generation logic.

The FPGA kernels are implemented by targeting a Xilinx ZCU104 development board. It contains a Zynq UltraScale+ MPSoC device (ZCU7EV) and one 2GB bank of 64-bit DDR4 memory. The GEMM kernel is synthesized into the FPGA logic using Vitis HLS version 2020.2. The software portion of our framework utilizes all four Arm Cortex-A53 cores running at 1.2GHz. Compiler auto-vectorization is enabled to use NEON instructions. The FPGA GEMM kernel uses a linear array of 16 PEs, with each PE containing 4 MAC units, for a total of 64 parallel MAC units running at 280MHz. In this system, each MAC sums an entire dot product and the number of addends depends upon the network model. 1) Multiplier: Conventionally, a floating-point product is computed by adding the exponents and multiplying the mantissas of two operands. With larger mantissa fields, the mantissa product is normally computed using the FPGA DSP blocks. Since we are working with small mantissa fields of 3 or fewer bits, we have instead used look-up tables (LUTs). We call this base implementation CFG1. Given that FPGAs typically contain 6-input LUTs, a single LUT can compute a single output bit of a 3 × 3 product. LUTs are far smaller than DSPs, which are normally intended for 16-bit or larger mantissas. For example, in a Xilinx UltraScale+ VUP13 FPGA, for every DSP block there are 140 6-input LUTs [START_REF] Xilinx | Higher Performance Neural Networks with Small Floating Point[END_REF]; these LUTs can implement 23 independent 3 × 3 integer multipliers.

We also implemented a fixed-point multiplier following standard conventions: the multiplier calculates the product of two integer values and rounds the result back to the input data type. This implementation uses DSP blocks because the required fixed-point formats are much wider.

Resource Utilization: Table I shows the resources used by a single multiplier using various formats. Floating-point formats are reported as EeMm, where e and m indicate the exponent and mantissa sizes. Fixed-point formats use Qi.f notation, where i and f are the integer and fractional part sizes. The precise formats reported were chosen because they achieved FP32-comparable results in our early experiments (not reported here). Note that Table I results include support for subnormals, NaNs, infinity and round-to-nearest.

To see the reduction in resource utilization, Table I also provides baseline values for 16-and 32-bit multipliers. The low-precision floating-point multipliers save considerable resources and notably use zero DSP blocks. The mantissa size more strongly affects area than the exponent length (e.g. E5M2 vs. E5M1 and E5M2 vs. E6M2 in Fig. 3). This matches expectations since extra mantissa bits in multipliers typically require more logic than extra exponent bits.

The fixed-point multipliers tend to use more LUTs than EeMm multipliers, plus they all need DSP blocks. Given that EeMm needs only 6 to 10 total bits while Qi.f needs 12 to 32 bits to represent values, Qi.f multipliers may also need more system support for wider buffering and interconnect than EeMm. Given that low-precision floating-point multipliers also have better range and good training results (shown later), fixed-point multiplication seems untenable.

2) Multiplier Variants: In this subsection, we enumerate several possible multiplier implementations. Their impact on DNN training is evaluated in Section IV-B2. Each variant below includes the changes of all previous variants.

CFG1 -Default: The first multiplier variant discussed in the previous subsection has ancillary logic to handle exceptional values, namely subnormal inputs and outputs, one output rounding mode, and NaN values.

CFG2 -Removal of Subnormal Outputs: The second multiplier variant removes support for subnormal output values. It saves several LUTs, yet it leads to information loss as subnormal outputs are all truncated to 0.

CFG3 -Removal of Output Rounding: The third multiplier variant produces an exact product by removing output rounding and widening the output sufficiently to avoid all precision loss. The product of two input values, each with m i + 1 mantissa bits after including the implicit bit, generates 2m i +2 output bits. Reducing this to fewer bits requires rounding and renormalizing, which incurs extra logic overhead. Removal of output rounding simultaneously reduces LUT count for multiplier and improves accuracy by avoiding all information loss. Note that output saturation to ±∞ remains, and inputs to the subsequent adder also change to 2m i + 2 bits.

CFG4 -Removal of NaNs: This variant removes NaN output values because they waste some encoding range. By treating NaNs as normalized values, we can extend the encoding range without additional logic. This requires remapping ±∞ to use a mantissa field of all ones, which is also more logical and consistent. Table II shows the expanded range.

CFG5 and CFG6 -Alternative Subnormal Inputs: The fifth and sixth multiplier variants attempt to save more LUTs by changing the way subnormal inputs are treated. Although this also affects the size of data conversion logic from FP32, the focus is on the size of MAC units which are far more numerous. Multiplication of denormalized input values requires dedicated logic, so CFG5 treats this part of the encoding range as normalized values to obviate the dedicated logic. Meanwhile, CFG6 truncates all subnormal input values to 0 so they become trivial.

Table II shows the custom encoding used for E5M2 subnormal input values. A CFG5 multiplier can still represent the upper portion of conventional subnormal numbers, but with extra precision. However, it loses all subnormal values with smaller exponent values. In contrast, CFG6 loses all subnormal values and is more likely to negatively impact training.

Resource Utilization: Fig. 3 shows the LUT counts used to implement the 6 multiplier variations across multiple data types: E6M2, E5M2, E6M1, and E5M1. For all data types, the LUT count of a single multiplier decreases from CFG1 to CFG6 as expected, and DSP count is always 0.

The first thing we notice is a considerable LUT count reduction when removing the logic for handling subnormal output values (i.e., CFG1 vs. CFG2).

Similarly, the reduction that arises from the removal of output rounding is shown when comparing CFG2 to CFG3. Interestingly, when the multiplier benefits a lot from the removal of subnormal output handling logic, for example saving more than 10 LUTs, its LUT count does not differ significantly between CFG2 and CFG3 (e.g. E6M1, and E5M1). On the other hand, when the reduction amount by the first modification is not remarkable, the multipliers tend to benefit more from removing the output rounding logic (e.g. E6M2). In all cases, multipliers having the same number of mantissa bits settle down to a similar LUT count after removing both subnormal outputs and output rounding. We believe this variation is caused by the sometimes inconsistent results of HLS synthesis optimizations.

In contrast, removing the extra shifters required by subnormal input values from CFG4 to CFG5 or CFG6 consistently decreases the LUT count by 10 to 13. Overall, both CFG5 and CFG6 are more than 50% smaller than CFG1, demonstrating a higher area sensitivity to multiplier policy than size in bits. Since CFG5 has a larger range with the same area, it is superior to CFG6.

3) Accumulator: Similar to the multiplier investigation, we studied accumulator resource use with low-precision floatingpoint and fixed-point data types. Both floating-point/fixedpoint accumulators are custom implemented following a standard approach for adders. Fixed-point accumulation includes saturation logic. Floating-point accumulation includes subnormals, swapping and shifting of operands plus extra bits (implied, guard, round, and sticky bits), hence using an m+4bit adder with m-bit mantissas. The baseline FP32 accumulator uses the Xilinx LogiCORE IP Library; it does not support subnormals and implements shifters in DSP blocks.

Resource Utilization: The Table III Accumulator columns show the resources used for different data types. For each floating-point accumulator, we also show the expected input type for the associated floating-point multiplier (CFG5). The fixed-point accumulators use considerably fewer LUTs than the floating-point ones because they do not require input/output alignment shifters or rounding. 4) Full MAC Unit: Based on the results so far, using a low-precision floating-point multiplier with a fixed-point accumulator uses the lowest LUT count for each operator. However, a MAC unit built this way also requires a floatto-fixed converter which uses costly data shifters.

The Table III Converter column shows the LUT count for converting the specified multiplier output precision to the Q8.13 fixed-point format; the Q8.13 precision was selected based upon training accuracy results in Section IV-B4.

The converters are larger than the Q8.13 accumulator itself. Using a total of 122 to 171 LUTs, they are closer to floatingpoint accumulators in size. Hence, the implementation cost of type conversion cannot be ignored.

IV. TRAINING RESULTS

To explore the accuracy impact of different low-precision GEMM configurations we trained ResNet-20 [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] and VGG16 [START_REF] Simonyan | Very Deep Convolutional Networks for Large-Scale Image Recognition[END_REF] networks on the CIFAR-10 [START_REF] Krizhevsky | Learning Multiple Layers of Features from Tiny Images[END_REF] dataset. On a slightly larger scale, we also train a ResNet-50 model using a subset of ImageNet. We use the GPU version for rapid experimentation; while the FPGA version works, GPUs are more readily available to run thousands of training experiments.

A. Training Settings

We have used SGD with a momentum value of 0.9 for all experiments. The initial learning rate and weight decay are set to 0.1/0.01 and 0.0001/0.0005 for ResNet-20 and VGG16 training, respectively, following the experimental settings of the original papers. All experiments employ a cosine annealing scheduler to decay the learning rate as the training proceeds. The batch size is set to 128, and we have trained the ResNet-20 model for 165 epochs and the VGG16 model for 200 epochs; both epoch counts are calculated based on the training iteration count reported in the original papers. ResNet-50 training uses an initial learning rate of 0.01 and runs for 100 epochs with a batch size of 64. All experiments also use an adaptive loss scaling technique (see Sec. II-B), with initial scaling factor S = 1024 and update frequency P = 200. For image preprocessing, the 8-bit pixel values are converted to FP32 by dividing by 255 and normalized using the mean and standard deviation of the dataset. When the first layer is lower precision, the input values are quantized accordingly.

B. Training Results

At first, we evaluate top-1 test accuracy after independently replacing either the multiplier or accumulator with lower precision. After that, we investigate accuracy of replacing the multiplier with one low-precision choice, and the accumulator with a different choice. The goal is to remain within 1 percentage point of FP32 results, which we call FP32-comparable.

1) Adaptive Loss Scaling: Results for using adaptive loss scaling with low-precision multiplication are shown in Fig. 4(a). Using more than 5 exponent bits were all FP32comparable (not shown), whereas E5Mm and E4Mm suffered from accuracy degradation (dashed lines) due to the limited representation range. By applying adaptive loss scaling (solid lines), gradients are shifted into the representable range and FP32-comparable results are restored for all except E4M1 and E4M2 multipliers with ResNet-20. Although not shown, low-precision accumulators showed similar sensitivity where loss scaling boosted accuracy to FP32-comparable results for E5Mm with 4/5 or more mantissa bits, and for E4Mm with 5/6 or more mantissa bits with VGG16/ResNet-20, respectively. Also not shown, non-adaptive loss scaling produced lower quality of results. We deem adapative loss scaling to be essential, so we adopt it for all subsequent experiments.

2) Multiplier Investigation: Fig. 4(b) shows test accuracy of E4Mm multiplier variants for different mantissa sizes. Removing subnormal output in CFG2 degrades accuracy across all mantissa sizes; this loss is restored with CFG3 by widening the output and removing rounding. Next, CFG4 reuses the NaN encoding space and improves accuracy slightly. On VGG16, both CFG5 and CFG6 are FP32-comparable, but on ResNet-20 CFG6 degrades unacceptably while CFG5 degrades only 2%. 1,2 Below, we will consider E4M2 with CFG5, since it is almost FP32-comparable, as well as E5M2.

3) Accumulator Investigation: Fig. 4(c) shows the test accuracy using different floating-point accumulators, while Fig. 4(d) uses fixed-point accumulators.

FP32-comparable training results are achieved with E7Mm/E6Mm/E5Mm accumulators using more than 5 mantissa bits, as well as E4Mm accumulators using more than 4/5 mantissa bits on ResNet-20/VGG16, respectively. As mentioned earlier, adaptive loss scaling was instrumental at improving results of E5Mm and E4Mm with larger m, but with smaller m problems persist, likely due to truncation error, also known as swamping [START_REF] Wang | Training Deep Neural Networks with 8-Bit Floating Point Numbers[END_REF], [START_REF] Higham | The Accuracy of Floating Point Summation[END_REF], when adding small values to a larger value. Compensated summation [START_REF] Higham | The Accuracy of Floating Point Summation[END_REF] can mitigate swamping, but widening the mantissa is easier.

With fixed-point accumulation, both ResNet-20 and VGG16 models can reduce the accumulator size to Q8.12 while maintaining FP32-comparable results. 4) MAC Investigation: Fig. 5 shows test accuracy when both multiplier and accumulator use lower precision. Observations from Sec. IV-B2 and Sec. IV-B3 lead us to start with an E4M2 (CFG5) multiplier and either E5M5 or Q8.12 accumulators. E5M5 accumulation converges but accuracy is 3.7% lower than FP32. Fixed-point accumulation fails to converge. Given the accuracy degradation in both cases, we first increased the exponent width of the multiplier to E5M2. This restores FP32-comparable accuracy, but fixed-point accumulation still fails to converge. Increasing the accumulator to Q8.13 achieves FP-32 comparable accuracy.

Table IV summarizes the system-level total LUT count and accuracy. Mixed-precision training requires around 25% fewer LUTs and no DSPs compared to the FP32 baseline, while limiting accuracy loss to within 1 percentage point.

Even after accounting for data conversion, we did not expect the fixed-point accumulator configuration to have a higher system-level LUT count. The increase was due to an increase in downstream interconnect and buffering logic from the wider output of the fixed-point accumulator (21 bits versus 12 bits for floating-point). To save downstream logic, we attempted to truncate the final accumulation value into a smaller fixed- point representation, but accuracy quickly degraded. There is an opportunity to explore other accumulation techniques, e.g. [START_REF] De Dinechin | An FPGA-Specific Approach to Floating-Point Accumulation and Sum-of-Products[END_REF], at the MAC and system level.

C. Extensibility to Larger Training Tasks

To test further, we used a subset of ImageNet called Imagewoof [START_REF] Howard | Imagenette, Imagewoof and Imagewang[END_REF] with a ResNet-50 model. Imagewoof reduces training effort compared to ImageNet by 99%, but keeps recognition accuracy challenging, by using just 10 dog classes out of the total 1,000 classes. Results in Table IV indicate the E6M5 accumulator configuration achieves FP32-comparable results, but the fixed-point accumulator requires 3 more bits (Q11.13) to converge.

From this, we suspect the size of a fixed-point accumulator is more sensitive to the model and dataset than a floatingpoint accumulator or multiplier; the lack of an exponent field requires a much wider word length. Future work on systemlevel techniques may reduce accumulator size and limit this sensitivity, e.g. chunk-based accumulation [START_REF] Wang | Training Deep Neural Networks with 8-Bit Floating Point Numbers[END_REF], or fixed-to-float conversion of accumulated results.

V. LIMITATIONS

The study is limited to only a few CNN models. We use smaller CNNs because our initial focus is on edge-based training, but larger CNNs are needed to validate the results. (We also trained MLPs on CIFAR10 and MNIST, but excluded those results for brevity.) However, the intent is to go beyond CNNs and explore other network types, such as transformers, and use-cases, such as transfer learning.

This work does not use larger datasets (e.g. full ImageNet) due to limitations with our training infrastructure. Being bitaccurate is very slow; we estimate a single ImageNet training instance to require over 60 days. Unfortunately our GPU systems impose time quota of 14 days which cannot be extended. To circumvent this limitation, we are adding checkpointing into our TinyDNN framework.

In addition, this work does not consider data formats other than floating-point and fixed-point formats. Other research has shown good potential for block floating-point [START_REF] Yang | SWALP: Stochastic Weight Averaging in Low Precision Training[END_REF], [START_REF] Das | Mixed Precision Training of Convolutional Neural Networks using Integer Operations[END_REF]- [START_REF] Drumond | Training DNNs with Hybrid Block Floating Point[END_REF], posits [START_REF] Gustafson | Beating Floating Point at its Own Game: Posit Arithmetic[END_REF]- [START_REF] Uguen | Evaluating the Hardware Cost of the Posit Number System[END_REF], logarithmic number systems [START_REF] Zhao | Low-Precision Training in Logarithmic Number System using Multiplicative Weight Update[END_REF]- [START_REF] Lee | LogNet: Energy-efficient neural networks using logarithmic computation[END_REF] and even binary formats [START_REF] Wang | Enabling Binary Neural Network Training on the Edge[END_REF]. The investigation of mixed data formats DNN training including those data types, especially for block mini-floats [START_REF] Fox | A Block Minifloat Representation for Training Deep Neural Networks[END_REF], is left for future work.

In addition, other than loss scaling, this work has not investigated significant changes to the training algorithms which may affect range and precision needs. This would likely improve results further. Although the work quantifies that FP8 saves bits and area, the resulting energy and memory footprint savings are left to future work as well.

VI. CONCLUSIONS

We have investigated the use of mixed-precision operations at the MAC level of a GEMM kernel during DNN training from a resource-accuracy tradeoff perspective. To do so, we have developed a framework, Archimedes-MPO, that allows the use of custom fixed and floating-point arithmetic at each network layer via custom GEMM kernels for GPU and FPGA. We then explored low-precision sensitivities on image classification tasks with the CIFAR-10 dataset on ResNet-20 and VGG16 CNNs, and then performed a final evaluation on ResNet-50 using a subset of ImageNet. Archimedes-MPO is availble under an open-source license.

Our results agree with prior work that 8-bit floating-point formats can be used to successfully train some CNNs, and that adaptive loss scaling is essential. However, we go beyond the bit width investigations done in previous work. For example, we provide evidence that avoiding certain IEEE-754 conventions for multiplication (e.g., output rounding, subnormals, NaN and ∞ encodings) at the MAC level can lead to both lower area and better accuracy.

In addition, when optimizing the × and + operators independently, we find that floating-point is more resource-efficient for multiplication, and its area is more sensitive to multiplier policy than size. In particular, we note that 8-bit floating-point formats use very little area because two 3-bit mantissas can be multiplied using 6-input LUTs, i.e., no DSPs are needed. In contrast, fixed-point is more efficient for accumulation, but its size is sensitive to the network model and dataset.

When combined at the MAC and system level, the area of data conversion and growth of downstream logic become a concern for fixed-point accumulation.

Despite the limited models and datasets used, we believe our overall conclusions to be sound: (1) multipliers should use a narrow floating-point format, (2) significant multiplier area can be saved, yet precision maintained, by eliminating output rounding and reducing exceptional value support, (3) for small mantissas, multiplier area is more sensitive to these policies than the precise mantissa size, and (4) fixed-point accumulators are very attractive (small size, no precision loss from rounding), but large structures are needed for format conversion (from the floating-point multiplier output) and systemlevel downstream logic (e.g. interconnect and buffering).

Although we did attempt to minimize the number of bits for low-precision training, we do not believe our values to be the 'final answer' due to previously stated limitations. To answer this, future work must investigate using a variety of network models, other training datasets, and changes to the training process, e.g., limiting accumulation chain lengths.

In our experiments, we demonstrated that the area used by the multiplier can be kept small if exceptional value support is modified and the number of mantissa bits is kept small. Still, it may be worthwhile to investigate other changes to training algorithms and number systems, such as multiplierless training where logarithmic numbers are used [START_REF] Zhao | Low-Precision Training in Logarithmic Number System using Multiplicative Weight Update[END_REF], [START_REF] Nair | Multiplierless MP-Kernel Machine For Energy-efficient Edge Devices[END_REF].

There are many directions for future work. Accumulation architecture at the MAC and system levels can be improved. For example, future work should consider breaking long accumulation chains, either in software or hardware, to limit the number of required accumulator bits. To save area, individual MAC adders and final accumulator outputs should be narrow while addressing the challenges of swamping, range and final test accuracy. Since GEMM hardware kernels can differ in how, where and when addition is done, there is a rich design space at the hardware level; this can be further augmented with software tools which break NN models into smaller chunks.

This work would also benefit from an error analysis to help fine-tune format requirements and further support our practical findings in a more principled manner. Beyond using adaptive loss scaling, assessing the resource-accuracy impact of other training computations (e.g. parameter update, activation function evaluation) is also important and should help guide the design of training algorithms and processes that are tailored for a low-precision regime.

Finally, we intend to expand upon our study by using more models (beyond CNNs) and datasets.

Fig. 1 .

 1 Fig. 1. Validation accuracy observed using LeNet5+MNIST.

Fig. 2 .

 2 Fig. 2. Archimedes-MPO Block Diagram for FPGAs.

Fig. 3 .

 3 Fig. 3. Area of floating-point multiplier variants.

Fig. 4 .

 4 Fig. 4. The graphs show the best top-1 test accuracy achieved by training with different configurations. The dashed red line always shows the best accuracy baseline achieved with FP32 training: (a) E5Mm (orange) and E4Mm (blue) multiplier accuracy before (dashed lines) and after (solid lines) adaptive loss scaling; (b) accuracy of E4Mm multiplier variants, where the blue/green/orange/cyan/magenta lines represent m=5/4/3/2/1 precisions; (c) accuracy of floatingpoint accumulators, where green/magenta/orange/blue lines represent e=7/6/5/4 exponent bits and the x-axis indicates the m mantissa bits; and (d) accuracy of fixed-point accumulators, where the x-axis indicates the number of fractional bits that follow 8 integer bits.

Fig. 5 .

 5 Fig. 5. Best CIFAR-10 test accuracy achieved with ResNet-20 training. The green/blue bars show the accuracy of using E5M5/E6M5 accumulators, and orange/red bars show the accuracy of using Q8.12/Q8.13 accumulators.

TABLE IV

 IV

	SYSTEM-LEVEL AREA AND TEST ACCURACY
				ResNet-20/	ResNet-50/
				CIFAR-10	Imagewoof
	configurations	LUTs	DSPs	Acc. (%)	Acc. (%)
	FP32	58,420	320	91.85	57.57
	E5M2(CFG5) + E6M5	42,640	0	91.04	59.79
	E5M2(CFG5) + Q8.13	43,471	0	90.95	10.92
	E5M2(CFG5) + Q11.13	44,876	0	not tested	59.38

• a new, extensible, open-source DNN training framework called Archimedes-MPO which uses FPGA or GPU

As explained in[START_REF] Tatsumi | Using Mixed Low-Precision Formats in Multiply-Accumulate (MAC) Units for DNN Training[END_REF], we adopt S = 4096 with CFG5/6 multipliers.

Although not shown, E5m is FP32-comparable for CFG5 and CFG6.