Lydia Abady

Jean-Stéphane Bailly

Nicolas Baghdadi

Yves Pastol

Hani Abdallah

Assessment of Quadrilateral Fitting of the Water Column Contribution in Lidar Waveforms on Bathymetry Estimates

Keywords: Lidar bathymetry, quadrilateral, recursive fitting, water column contribution, waveform processing

A new approach based on a mixture of Gaussian and quadrilateral functions was developed to process bathymetric lidar waveforms. The approach was tested on two simulated data sets obtained from the existing Water-LIDAR (Wa-LID) waveform simulator. The first simulated data set corresponds to a sensor configuration modeled after a possible future satellite bathymetric lidar sensor that was previously studied. The second simulated data set corresponds to a lidar airborne configuration modeled using the HawkEye airborne lidar parameters. In the proposed approach, the lidar waveform is fitted into a combination of three functions, two Gaussians for both the water surface and water bottom contributions and a quadrilateral function to fit the water column contribution. The results show more accurate bathymetry estimates compared with the use of a triangular function to fit the column contribution or a simple peak detection method. For the satellite configuration, the bias is improved by 16.8 and 0.8 cm compared with the peak detection method and the use of a triangular function, respectively. For the airborne configuration, the bias is improved by 10.0 and 2.4 cm compared with the peak detection method and the use of a triangular function, respectively. The proposed waveform fitting using the quadrilateral function underestimates the bathymetry by -5.0 and -6.1 cm for the simulated satellite and airborne data sets, respectively. The standard deviations of the bathymetry estimates are 6.0 and 8.2 cm, respectively. The obtained biases are inherent to overlaps between functions fitting the water surface, column, and bottom contributions.

I. INTRODUCTION

IDAR is an active remote sensing system that can be used to estimate altimetry and bathymetry. However, to date, there has not been any lidar satellite mission usable to bathymetry due to low laser power or unsuitable wavelengths to penetrate water [START_REF] Battaglia | Utility of satellite LIDAR waveform data in shallow water[END_REF]. Nevertheless, airborne lidar bathymetry (ALB) is recognized as an efficient tool that is able to swiftly collect information to evaluate the effects of natural disasters (flooding and hurricanes) and is also beneficial in hydrography and hydrological applications (see [START_REF] Abdallah | Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters[END_REF], [START_REF] Battaglia | Utility of satellite LIDAR waveform data in shallow water[END_REF], and [START_REF] Mohammadzadeh | A state of art on airborne lidar application in hydrology and oceanography: A comprehensive overview[END_REF]). The water depth is usually computed by subtracting the bottom return time from the surface return time of the laser beam and then multiplying this time difference by half the speed of light in water [START_REF] Irish | An Introduction to Coastal Zone Mapping With Airborne LIDAR: The SHOALS System[END_REF]. Typically, a part of the laser energy is lost due to Fresnel reflectance at the water surface, scattering, and absorption in the water, which will limit the strength of the bottom return and, hence, limit the maximum detectable depth (see [START_REF] Irish | An Introduction to Coastal Zone Mapping With Airborne LIDAR: The SHOALS System[END_REF] and [START_REF] Wang | Using airborne bathymetric lidar to detect bottom type variation in shallow waters[END_REF]). However, the most limiting factor of lidar bathymetry is the optical water clarity depending on scattering and absorption within water [START_REF] Irish | An Introduction to Coastal Zone Mapping With Airborne LIDAR: The SHOALS System[END_REF], which affects the water column contribution in the waveform, registering the backscattered energy along time.

A typical bathymetric lidar waveform contains a surface contribution, a column backscattered contribution, a bottom contribution, and background noise contributions. To retrieve a more accurate bathymetry from bathymetric waveforms, several postprocessing approaches can be used. Among these approaches are fitting process and mathematical approximations, heuristic methods (such as the simple peak detection method), statistical approaches, or deconvolution methods [START_REF] Allouis | Comparison of LiDAR waveform processing methods for very shallow water bathymetry using Raman, near-infrared and green signals[END_REF]. However, despite being efficient, fitting approaches do not completely capture the physics of the interaction of the laser beam with the water surface, column, and bottom. Many previous fitting approaches have been applied to represent lidar bathymetric waveforms. For example, two Gaussian functions were used to represent the surface and bottom returns. Because the column contribution shows an asymmetric shape physically corresponding to an infinite sum of successive translated Gaussian functions with an exponentially decreasing amplitude along time, i.e., water depth, it was previously approximated as a mixture of Gaussians in [START_REF] Long | What Optech's bathymetric LiDAR sees underwater[END_REF]. In Abdallah et al. [START_REF] Abdallah | Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters[END_REF], a Gaussian function, a triangular function, and a Weibull function were used to fit the surface contribution, the column contribution, and the bottom contribution, respectively. However, even with more suitable functions to fit each contribution, a bias in bathymetry estimates was also observed [START_REF] Abdallah | Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters[END_REF].

In this letter, we aim to define a new waveform fitting algorithm to retrieve bathymetry where 1) a quadrilateral function is used to better fit the water column contribution in addition to two Gaussian functions used to fit the surface and bottom contributions, and 2) it takes advantage of a nonlinear recursive optimization algorithm [START_REF] Griva | Nonlinear optimization[END_REF]. The benefit of the new fitting approach we propose is assessed on sets of both satellite and airborne simulated lidar waveforms. The results of bathymetry estimates (biases and standard deviations) computed using a quadrilateral function were compared with those obtained from a simple peak detection approach and using the triangle function fitting method initially proposed by Abdallah et al. [START_REF] Abdallah | Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters[END_REF]. Despite being basic and not robust compared with the operational methods of bathymetric lidar systems that use the correction calibration process to address biases in bathymetry (see [START_REF] Guenther | Airborne Laser Hydrology System Design and Performance Factors[END_REF] and [START_REF] Guenther | Meeting the accuracy challenge in airborne LiDAR bathymetry[END_REF]), a water depth estimate from a peak detection method was used here as the reference method to highlight the benefit of any waveform fitting method. Section II introduces the data sets, whereas Section III offers a brief presentation on the methodology used and the results obtained.

In the last section, a summary of the results is given, and a discussion on bias sources is proposed.

II. DATA SETS DESCRIPTION

Three waveforms data sets were generated in using the recently developed Wa-LID code [START_REF] Abdallah | Wa-LiD: A new LiDAR simulator for waters[END_REF]. For each of the three data sets, the waveforms are the summation of the laser pulse convolution with the impulse response functions of the surface, column, and bottom, as well as additional noises.

The simulator governs both the water properties and the used sensor's specifications. The simulator also models solar radiation and detector noise [START_REF] Abdallah | Wa-LiD: A new LiDAR simulator for waters[END_REF]. Moreover, it includes the effects of pulse stretching, where the bottom return is stretched due to bottom slopes and incidence angle, reducing its amplitude. However, the Wa-LID simulator assumes that the water column is homogeneous, with the same optical properties for the whole column.

Waveforms were generated in using probability distributions for each water parameters. These distributions came from databases collected on coastal waters around the world and from published references [START_REF] Abdallah | Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters[END_REF]. They are assumed to be representative of what we can encountered at the global scale. The probability distributions of water properties that were used to generate the databases of simulated waveforms are presented in Table I.

We used parameters that correspond to two sensors, namely, satellite and airborne (see Table II). The parameters used for satellite data set modeling were chosen by the European Aeronautic Defence and Space (EADS)-Astrium Company [START_REF] Abdallah | Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters[END_REF] as a possible future satellite lidar bathymetry system. Airborne data set parameters were modeled using HawkEye sensor parameters [START_REF] Beckmann | The Scatter of Electromagnetic Waves From Rough Surfaces[END_REF]. In our simulations, we used the usual green wavelength (532 nm) common to many ALBs [START_REF] Axelsson | Rapid topographic and bathymetric reconnaissance using airborne LiDAR[END_REF] because of the necessity for a tradeoff between the attenuation of the bottom return due to water absorption and the scattering and absorption by particles in the water. The field-of-view (FOV) angle for the HawkEye system ranges between 10 and 50 mrad. [START_REF] Abdallah | Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters[END_REF], [START_REF] Tulldahl | Analytical waveform generation from small objects in LiDAR bathymetry[END_REF], [START_REF] Guenther | Airborne Laser Hydrology System Design and Performance Factors[END_REF] The lower the FOV is, the lower is the solar radiation noise. For a very small FOV, however, the bottom contribution's power magnitude will be reduced, resulting in a lower maximum detectable depth. Thus, an FOV of 15 mrad was chosen for the airborne simulation. For the sake of simplicity, we used an onnadir rather than an off-nadir angle.

In order to compute accuracy statistics on the bathymetry estimates (biases and standard deviations), the first data set was generated. Data set 1 contains 32 000 simulated waveforms, 16 000 waveforms for each sensor (airborne and satellite). For each sensor, the 16 000 simulated waveforms come from 2000 waveforms for each one of the eight investigated water depths: 1, 2, 3, 4, 5, 6, 7, and 10 m. At a given water depth, the 2000 waveforms were generated with random water properties, according to the distributions in Table I.

In order to better investigate if the standard deviation in the bathymetry estimates comes from water column properties variations, from noise, or from the fitting process, a second data set was generated. Data set 2 contains 16 000 simulated waveforms, 8000 waveforms for each one of the two sensors. For each sensor, the 8000 simulated waveforms come from 1000 waveforms for each one of the eight investigated water depths. At a given water depth, the 1000 generated waveforms differ from data set 1 by fixing the water column parameters to slightly turbid water (Yellow substance Ay0 = 0.07 m -1 , Phytoplankton Ph = 1 mg/m 3 , and Sediment S = 2.6 mg/l) [START_REF] Smith | Optical properties of the clearest natural waters (200-800 nm)[END_REF]. The choice of slightly turbid water is to ensure better bottom peak detection at greater water depths. Computing the ratio between the obtained standard deviations of bathymetry estimates with variable (data set 1) or fixed water column properties (data set 2) allows measuring the sensitivity of standard deviation estimates to the variability of water column properties.

Finally, in order to verify that the bias of bathymetry estimates mainly comes from the water column return overlapping the surface and bottom water returns, a smaller third data set was generated. Data set 3 contains only 4000 simulated waveforms. Two thousand waveforms were generated for each sensor using a water depth of 3m. Indeed, these 4000 waveforms are a subset of the one simulated in data set 1 without the b t c ≤ ≤ column contribution. This data set assumes a totally transparent water column. Next, these waveforms were processed with a bi-Gaussian fitting process (fitting both surface and bottom returns) in order to verify that the bathymetry bias disappears when the water column contribution is removed. For all simulated waveforms (data sets 1, 2, and 3), the effects of detector noise and noise due to solar radiation were integrated.

III. METHODOLOGY

A. Peaks Detection Procedure

To perform a water depth estimate, the lidar waveform should exhibit two distinct detectable peaks, representing the surface and bottom echoes. To facilitate the identification of the two peaks, denoising of lidar signal was applied using wavelet decomposition based on Stein's heuristic [START_REF] Qiang | A noise reduction approach based on Stein's unbiased risk estimate[END_REF]. This step is essential in order to limit the effect of the solar radiation, particularly on the airborne waveforms. Next, a peaks detection procedure was used on the smoothed waveforms. It considers that a peak is any local maxima in the lidar waveform that has an amplitude much higher than the noise level (10 times for example). After identification of all peaks in the waveform, the peak with the smallest time bin is attributed to the surface position, and the peak with the largest time bin is attributed to the bottom position. Other peaks between the surface and bottom positions are commonly attributed to turbidity, marine life, or suspended particles in the water column [START_REF] Kinzel | Evaluation of an experimental LiDAR for surveying a shallow, braided, sandbedded river[END_REF].

B. Waveform Fitting Procedure

Next, a fitting procedure is applied to the waveforms. In this letter, we experimented with two fitting schemes that are both dependent on a nonlinear recursive least squares fitting algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF], and we considered the surface and bottom returns to be Gaussian functions. However, the difference lies in the presentation of the water column: In one scheme, it is fitted into a triangular function [see Fig. 1(a)]; in the newly proposed scheme, it is fitted into a quadrilateral function [see Fig. 1(b)]. Fig. 1(c) shows that both fittings are superimposed on the simulated waveform without noise, which demonstrates efficient fitting in this schematic case. However, the quadrilateral function shows a better fitting to the shape of the water column contribution than the triangular function.

To evaluate the performance of those two fitting algorithms, the bias, i.e., the mean difference between the estimated and simulated water depths, and the standard deviation of the The fitted column return f 2 (t) is given by bathymetry estimates were computed for each data set (satellite and airborne configurations).

For the proposed fitting method using a quadrilateral function

⎧ ⎪ e ⎪ 0 t ≤ a t-a a ≤ t ≤ b ⎨ ec bg b-a e)
to fit the water column contribution, the fitted lidar waveform can be expressed as

R(t) = f 1 (t) + f 2 (t) + f 3 (t). f 2 (t) = ⎪ ⎩ -+t(g- c-b ≤ ≤ g d-t c t d d-c 0 t ≥ d
The fitted surface return f 1 (t) is given by

(t-μ 1 ) 2
where a, b, c, and d are the x-axis points for the quadrilateral; and e and g are the ordinates for the third and fourth points, respectively.

f 1 (t) = a 1 - √ e 2σ 2 1
The fitted bottom return f 3 (t) is given by where a 2 , μ 2 , and σ 2 are the amplitude, the mean and the standard deviation of the Gaussian function, respectively.

σ 1 2π a 2 (t-μ 2 ) 2 -

IV. RESULTS

A. Satellite Sensor

Of the overall satellite waveforms, 39.6% were detectable, i.e., having identified peaks. Fig. 2 shows the detection rate at each depth. As expected, the detection probability increases until it reaches a water depth of 3 m as a result of the dependence of the detection rate on the merger of the surface return with the bottom's return at shallow depths. At water depths greater than 3 m, the peak detection rate decreases as the water depth increases due to attenuation and absorption.

After applying peak detection, the fitting procedure was applied to the waveforms that had a detectable bottom. Fig. 3 shows that the bias insignificantly increases as the water depth increases, where it increased from -20.3 cm at a water depth of 1 m to -25.0 cm at a depth of 10 m using the peak detection method. When using quadrilateral fitting, the bias increased from -2.0 cm at a water depth of 1 m to -8.5 cm at a depth of 10 m. For the triangular fitting, the bias increased from -2.7 cm at a water depth of 1 m to -9.4 cm at a depth of 10 m. This increase in bias is most likely due to the weaker bottom return at greater water depths. Fig. 3 also displays the standard deviation. Upon examination, we notice that, for the quadrilateral fitting, the standard deviation is slightly lower than that for the triangular fitting. Using the quadrilateral function, the standard deviation ranges between 4.8 and 8.3 cm. However, using the triangular function, the standard deviation ranges between 5.0 and 8.4 cm. The overall bias and standard deviation for all the used water depths are shown in Table III, which clearly portrays the improvement in reducing the bias by 16 cm in magnitude when using the fitting procedure based on the quadrilateral function.

Next, to scrutinize how far the overall 6-cm standard deviation obtained on data set 1 (see Table III) is caused by the waveform variability due to the water column properties variability, the standard deviation of the bathymetry estimates after applying the quadrilateral function algorithm on data set 2 was computed. The obtained overall standard deviation of bathymetry estimates is thus 1.3 cm. When comparing this 1.3 cm to the 6 cm obtained on data set 1, we can deduce that 78% (equal to 1-1.3/6) of the 6-cm standard deviation is a result of the variability of water column parameters, whereas the remaining part (22%) can be due to the variability in the fitting due to random noise. Of course, this result is dependent to the studied sensor characteristics and the water column parameters distribution we chose.

In order to verify that the bias of bathymetry estimates mainly comes from the water column return overlapping the surface and bottom returns, we computed the bias using a bi-Gaussian function fitting of the 4000 waveforms from the third data set (one Gaussian for each of the surface and bottom peaks). Next, we compared this bias to the one computed using the 4000 waveforms at a water depth of 3 m from the first data set (waveforms with surface, bottom, and column contributions). The obtained bias was near zero for data set 3 and was approximately -4.5 cm for data set 1. This result leads us to consider that the bias is inherent to the summation of the three contributions (surface, column, and bottom) that are overlapping. These overlaps cause a translation in time of the peaks, which do not correspond to central surface and bottom positions. To reduce this bias, a statistical approach introducing a calibrated offset parameter could be designed to better fit the data during actual surveys.

B. Airborne Sensor

For the airborne data set, the overall detection rate was 24.4%. The peak detection rate (see Fig. 2) has a similar distribution to that from the satellite data set, except that the detection rate starts decreasing at a water depth of 2 m, which can be explained by the differences in airborne and satellite sensor parameters (e.g., T 0 and H ). Indeed, the overlap between column and bottom contributions increases when the emitted pulse (T 0 ) increases. Moreover, the sensor altitude (H ) affects the results by acting on contribution (surface, column, and bottom) magnitude due to air column transmittance.

The bias at each water depth for the airborne data set (see Fig. 3) also slightly increases as the water depth increases, ranging from -5.7 to -9.6 cm for the quadrilateral fitting, which is better than -8.7 to -11.5 cm obtained for the triangular fitting. Additionally, it shows an obvious improvement compared with using only the peak detection method, which increased from -16.0 cm at a water depth of 1 m to -19.1 cm at a water depth of 10 m. Once again, the performance of the quadrilateral fitting is slightly better than the triangular fitting. Fig. 3 also displays the standard deviation, which is nearly identical over all water depths, ranging from 7.3 to 8.7 cm when using the quadrilateral fitting and 7.8 to 9.4 cm when using the triangular fitting. The increase in the standard deviation compared with the satellite data set is due to higher solar radiation noise in the airborne data set. The overall bias (see Table III) shows an improvement of approximately 10 cm in the bias when using the quadrilateral fitting procedure compared with simple peak detection. In order to also assess how far the 8.2-cm (see Table III) overall standard deviation in bathymetry estimates from the airborne data set is caused by the water column properties variability, a similar computation to that performed for the satellite data set was computed. The computed standard deviation of the bathymetry estimates after applying the quadrilateral function algorithm on data set 2 is 5.1 cm. Compared with the 8.2 cm obtained on data set 1 (see Table III), we can deduce that only 38% of the 8.2-cm standard deviation is a result of the variability of water column parameters.

To ensure that the bias is also coming from the column contribution overlapping the surface and bottom returns, as in the satellite data set, we used the simulated airborne data set 3 at a water depth of 3 m. Once again, the bias computed on data set 3 was negligible. This result proves that, for this data set as well, the source of the bias was due to the overlap of the column contribution on the surface and bottom returns.

V. CONCLUSION

In this letter, a new lidar fitting algorithm to retrieve bathymetry has been tested on simulated waveform data sets. Two lidar sensor configurations were used to simulate these waveforms, one using satellite sensor parameters and the other using airborne sensor parameters. The bathymetry detection rate was established depending on the number of waveforms that incorporate more than one peak. The estimation accuracy was evaluated based on the bias and standard deviation for each data set (satellite and airborne). It is highlighted that the bias depended on the column contribution being added to the surface and column returns, which caused a shift in the surface and bottom peaks that led to a bias in the water depth estimate. However, the fitting algorithm we proposed showed a clear improvement compared with the water depth estimate coming from simple peak detection, and it showed a slight improvement over fitting alternatives using a triangular function. For the simulated satellite data set using the quadrilateral function, the results showed an improvement of 16.8 cm in the bias compared with the basic peak detection method. For the simulated noisier airborne data set, the improvement in bias was lower, at approximately 10.0 cm. Subsequently, the quadrilateral function produced slightly more unbiased estimations than the triangular function in both data sets. Moreover, the fitting using a quadri-lateral function showed a slightly better standard deviation, with 6.0 cm for the satellite data set and 8.2 cm for the airborne, compared with the fitting using a triangular function, with 7.0 cm for the satellite data set and 8.5 cm for the airborne data set.

Fig. 1 .

 1 Fig. 1. (a) Triangular and (b) quadrilateral fitting of the column contribution for a simulated waveform without noise using satellite sensor parameters at a semipure water depth of 2 m (2 m = (38 -21) × 0.225 × 0.5) sampled at 1 ns. (c) Quadrilateral and triangular fittings superimposed on the simulated contributions.

Fig. 2 .

 2 Fig. 2. Peak detection rate (in percentage) versus water depth (in meters).

Fig. 3 .

 3 Fig. 3. Bias and standard deviation on bathymetry estimates (in centimeters) versus water depth (in meters) (computed on data set 1).

TABLE I PRO

 I BA BILITIES DISTRIBU TIONS OF WATER PARAMETERS O N COA S TA L WATERS[START_REF] Abdallah | Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters[END_REF]

TABLE II SENSOR

 II PARAMETER SPECIFI CATIONS

TABLE III OVERALL

 III BAT HYMETRY ACCURACY (BIASES AND STANDARD DEVIATIONS IN CENTIMETERS ) USING ALL INVESTIGATED WATER DEPTHS ) (COMPUTED ON DATA SET 1)

where a 1 , μ 1 , and σ 1 are the amplitude, the mean, and the standard deviation of the Gaussian function, respectively. f

(t) = √ e σ 2 2π

ACKNOWLEDGMENT

The authors would like to thank F. Fabre from the EADS Company for his useful advice.


Antoine Lejay
email: antoine.lejay@univ-lorraine.fr

Estimation of the bias parameter of the Skew Random Walk and application to the Skew Brownian motion

Keywords: Skew random walk, Skew Brownian Motion, maximum likelihood estimator, local asymptotic mixed normality, local time, null recurrent process

We study the asymptotic property of simple estimator of the parameter of a Skew Brownian Motion when one observes its positions on a fixed grid -or equivalently of a simple random walk with a bias at 0. This estimator, nothing more than the Maximum Likelihood Estimator, is based only on the number of passages of the random walk at 0. It is very simple to set up, is consistent and is asymptotically mixed normal. We believe that this simplified framework is helpful to understand the asymptotic behavior of the maximum likelihood of the Skew Brownian Motion observed at discrete times which is studied in a companion paper.

Published in

Statistical Inference for Stochastic Processes doi:10.1007/s11203-017-9161-9

1

A flourishing literature is devoted to the stochastic diffusion processes with discontinuous coefficients or involving local times. This covers many fields of application, ranging from geophysics [START_REF] Appuhamillage | Occupation and local times for skew Brownian motion with applications to dispersion across an interface[END_REF] to population ecology [START_REF] Lejay | Identification of Diffusion Coefficient in Nonhomogeneous Landscapes[END_REF][START_REF] Barahona | A simulation-based study on Bayesian estimators for the skew Brownian motion[END_REF][START_REF] Cantrell | Diffusion Models for Population Dynamics Incorporating Individual Behavior at Boundaries: Applications to Refuge Design[END_REF] through molecular biology [START_REF] Bossy | Probabilistic interpretation and random walk on spheres algorithms for the Poisson-Boltzmann equation in molecular dynamics[END_REF], or to study rank based statistics [START_REF] Fernholz | Two Brownian particles with rank-based characteristics and skew-elastic collisions[END_REF] as well as dealing with the tiled volatility model [START_REF] Lipton | Filling the Gap[END_REF] or directional predictability [START_REF] Alvarez | Timing in the Presence of Directional Predictability: Optimal Stopping of Skew Brownian Motion[END_REF] in mathematical finance.

In this direction of research, the Skew Brownian motion (SBM) is a process of utter importance. It plays the role of the Brownian motion for diffusions with continuous coefficients. Informally, the SBM is similar to the Brownian motion, except that when it hits 0, its behavior depends on a parameter 𝛼 ∈ (0, 1) which gives a trend to go on the positive or negative axis according to 𝛼 > 1/2 or 𝛼 < 1/2. As pointed out for example in [START_REF] Lejay | Simulating diffusion processes in discontinuous media: a numerical scheme with constant time steps[END_REF], Lamperti's type transforms allow one to reduce one-dimensional diffusion processes generated by either 𝑎△ or ∇(𝑎∇•) with 𝑎 piecewise constant to processes which behave locally as a SBM.

Formally, the SBM of parameter 𝛼 is the strong solution to the stochastic differential equation

𝑋 𝑡 = 𝐵 𝑡 + 𝜃𝐿 𝑡 , 𝜃 = 1 + 𝛼 2 ∈ (-1, 1), (1)
where 𝐵 is a Brownian motion and 𝐿 𝑡 the symmetric local time of 𝑋 at 0. Introduced through its infinitesimal generator by K. Itô and H.P. McKean (See e.g. [22, p. 220] or [23, Problem 1, p. 115]), there exist many ways to describe this process (See [START_REF] Lejay | On the constructions of the skew Brownian motion[END_REF] for a survey of some possible constructions).

The existence of a strong solution to [START_REF] Lejay | Identification of Diffusion Coefficient in Nonhomogeneous Landscapes[END_REF] has been proved first by J.M. Harrison and L.A. Shepp in [START_REF] Harrison | On skew Brownian motion[END_REF].

In the same way the Brownian motion is approximated by a simple random walk properly renormalized, the SBM is approximated by a Skew Random Walk (SRW) [START_REF] Harrison | On skew Brownian motion[END_REF]. A SRW is a simple random walk {𝜁 𝑘 } 𝑘≥0 with a bias at 0:

P[𝜁 𝑘+1 = 𝑖 |𝜁 𝑘 = 𝑗] = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1/2 if 𝑗 ̸ = 0, 𝑖 = 𝑗 ± 1, 𝛼 = (1 + 𝜃)/2 if 𝑗 = 0, 𝑖 = 1, 1 -𝛼 = (1 -𝜃)/2 if 𝑗 = 0, 𝑖 = -1, 0 otherwise. ( 2
)
Let us fix a parameter ℎ > 0 which represents the size of a grid ℎZ that discretize the state space. The successive positions {𝜉 𝑘 (ℎ)} 𝑘≥0 at which a SBM 𝑋 passes through the grid ℎZ, up to time 𝑇 , are observed and recorded (the number of recorded points is random 𝑛 ⋆ (ℎ), yet close to 𝑇 /ℎ 2 ). By doing so, the embedded Markov chain {𝜉 𝑘 (ℎ)} 𝑘≥0 is equal in distribution to the rescaled SRW {ℎ𝜁 𝑘 } 𝑘≥0 .

Our main concern is to estimate 𝛼 = (1 + 𝜃)/2 from the observations of the successive positions of a SBM on the grid ℎZ. The times at which the observations occurs are then random times. We study the asymptotic behavior of the estimators as the size ℎ of the spatial grid decreases to 0.

From now, we denote by P 𝛼 the distribution of the SBM of parameter 𝛼 ∈ (0, 1).

Since P 𝛼 [𝑋 𝑡 ≥ 0] = 𝛼 for any 𝑡 > 0, E 𝛼 [𝐴 𝑡 ] = 𝑡𝛼 where 𝐴 𝑡 is the occupation time 𝐴 𝑡 above 0 up to time 𝑡. However, 𝐴 𝑇 /𝑇 is a meaningless estimator of 𝛼 since by scaling, 𝐴 𝑇 /𝑇 law = 𝐴 1 whatever 𝑇 > 0 and 𝐴 1 follows a variant of the Arc Sine distribution [START_REF] Lamperti | An occupation time theorem for a class of stochastic processes[END_REF]Theorem 1]. Unlike most of the work in statistics of diffusion, ergodic theory could not be used here as the SBM is a null recurrent process.

The sole information useful for our estimator are the successive points (𝜉 𝑘 (ℎ), 𝜉 𝑘+1 (ℎ)) for which 𝜉 𝑘 (ℎ) = 0. The excursions of the SRW behave like the ones of a simple random walk. They carry no information about 𝛼, except from their sign.

Our estimator of 𝛼 is then simply

𝛼(ℎ) = #observed upcrossings before 𝑇 #observed passages to 0 before 𝑇 , (3)
where a passage to 0 counts to 1 when 𝜉 𝑘 (ℎ) = 0 and an upcrossing counts to 1 when 𝜉 𝑘 (ℎ) = 0 and 𝜉 𝑘+1 (ℎ) = 1. This estimator (3) is also the empirical mean number of positive excursions. It is the Maximum Likelihood Estimator (MLE) of 𝛼 for the SRW.

First, we show that 𝛼(ℎ) is consistent, i.e. converges in probability to 𝛼 as ℎ → 0 under P 𝛼 for any 𝛼 ∈ (0, 1). Second, we show that under P 𝛼 ,

ℎ -1/2 (𝛼(ℎ) -𝛼) converges in distribution to √ 𝛼(1 -𝛼) 𝑊 (𝐿 𝑇 ) 𝐿 𝑇 ,
where 𝑊 is a Brownian motion independent from 𝐵. Asymptotically, 𝛼(ℎ) behaves like a mixed normal distribution. The number 𝑛 ⋆ (ℎ) of observed positions is random, yet of order 𝑇 /ℎ 2 . For 𝑛 ⋆ (ℎ) observations, the rate of convergence is not

𝑛 ⋆ (ℎ) -1/2 (hence ℎ), but 𝑛 ⋆ (ℎ) -1/4 (hence √ ℎ).
The reason is that the number of passages to 0 is roughly of order √ 𝑛 ⋆ (ℎ). We then show that the distribution of the {𝜉 𝑘 (ℎ)} 𝑘=0,...,𝑛 ⋆ (ℎ) satisfies the Local Asymptotic Mixed Normal (LAMN) property [START_REF] Jeganathan | On the asymptotic theory of estimation when the limit of the log-likelihood ratios is mixed normal[END_REF][START_REF] Jeganathan | Some asymptotic properties of risk functions when the limit of the experiment is mixed normal[END_REF] -a generalization of the Local Asymptotic Normal (LAN) property [START_REF] Ibragimov | Has'minskiȋ. Statistical estimation[END_REF][START_REF] Lecam | On some asymptotic properties of maximum likelihood estimates and related Bayes' estimates[END_REF] -which carries some information about any asymptotically efficient estimator.

We are in a situation close to the one studied by D. Florens in [START_REF] Florens | Estimation of the diffusion coefficient from crossings[END_REF] where non-parametric estimation of the diffusion coefficient of the solution of a stochastic differential equation observed at discrete time is performed through the number of crossings of levels. The kind of limit theorem we obtain here is typical for of null recurrent "regular" diffusion processes [START_REF] Delattre | Asymptotic equivalence for a null recurrent diffusion[END_REF], although the approach shall be adapted to deal with distributions which are singular with respect to the Wiener measure [START_REF] Gall | One-Dimensional Stochastic Differential Equations Involving the Local Times of the Unknown Process[END_REF].

In [START_REF] Lejay | Is a Brownian motion skew?[END_REF], we have studied the MLE estimator for the SBM observed at frequency 𝑇 /𝑛 when 𝛼 = 1/2 (or 𝜃 = 0). By using results from J. Jacod in [START_REF] Jacod | Rates of convergence to the local time of a diffusion[END_REF], we obtained similar results, both to identify the limit and the rate of convergence. The consistency of the MLE estimator when 𝛼 ̸ = 1/2 (or 𝜃 ̸ = 0), which was left open, is now treated in the subsequent article [START_REF] Lejay | Two consistent estimators for the Skew Brownian motion[END_REF]. However, following [START_REF] Jacod | Rates of convergence to the local time of a diffusion[END_REF], these articles use estimates on the semi-groups and martingales limit theorems.

As shown in [START_REF] Lejay | Two consistent estimators for the Skew Brownian motion[END_REF], the MLE estimator 𝛼 𝑛 of 𝛼 for the SBM observed at discrete times {𝑖𝑇 /𝑛} 𝑖=0,...,𝑛 solves

𝛼 𝑛 = E 𝛼𝑛 [𝑈 |{𝑋 𝑖𝑇 /𝑛 } 𝑖=0,...,𝑛 ] E 𝛼𝑛 [𝑈 + 𝐷|{𝑋 𝑖𝑇 /𝑛 } 𝑖=0,...,𝑛 ] (4)
where 𝑈 (resp. 𝐷) is the number of intervals 𝐽 𝑖 = [𝑖𝑇 /𝑛, (𝑖 + 1)𝑇 /𝑛] for which 𝑋 𝑡 = 0 for some 𝑡 ∈ 𝐽 𝑖 and 𝑋 (𝑖+1)𝑇 /𝑛 > 0 (resp. 𝑋 (𝑖+1)𝑇 /𝑛 < 0). With the above construction, the fact that the process crosses 0 or not is considered as some hidden data (actually, (4) was derived using some ideas of the Expectation-Minimization algorithm).

In our situation, the situation is much more simpler because the crossing are observed so that the number of up-and down-crossings are exactly counted. Our MLE estimator is nothing more than a binomial one with a random number of samples. The asymptotic error, given by a CLT-like result, steams from the Berry-Esseen theorem and the P. Lévy downcrossing theorem, which relates the local time as a limit of the renormalized number of crossings of the SRW.

We then believe that our results here explain with simple computations and elementary tools the results in [START_REF] Lejay | Is a Brownian motion skew?[END_REF][START_REF] Lejay | Two consistent estimators for the Skew Brownian motion[END_REF]. Also, it could be used for other statistics involving excursions such as the ones in [START_REF] Lejay | Statistical estimation of the Oscillating Brownian Motion[END_REF].

Finally, let us mention that this approach fits well with a Bayesian framework especially using Beta random variables for prior and posterior distributions. This allows non-asymptotic point estimates and probability intervals [START_REF] Robert | The Bayesian choice. From decision-theoretic foundations to computational implementation[END_REF]. Regarding the SBM, this was addressed in [START_REF] Barahona | A simulation-based study on Bayesian estimators for the skew Brownian motion[END_REF] for observations on a fixed time grid. Statistical tests were then successfully applied to the location of South America sea lions in southern Chile.

Outline. Section 1 gives a quick account on the SBM and the SRW. The asymptotic of the estimator 𝛼(ℎ) in ( 3) is studied in Section 2. In Section 3, the LAMN property is proved.

Approximation of the Skew Brownian motion by the Skew Random Walk

The symmetric local time of a semi-martingale 𝑌 is defined as

ℒ 𝑡 (𝑌 ) = lim 𝜖→0 1 2𝜖 ∫ 𝑡 0 1 𝑌𝑠∈[-𝜖,𝜖] d⟨𝑌 ⟩ 𝑠 , 𝑡 ≥ 0. (5)
All the results we need here on the local time could be found for example1 in [START_REF] Itô | Diffusion processes and their sample paths[END_REF] or [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF].

Let 𝑋 be SBM of parameter 𝛼 ∈ (0, 1) with distribution P 𝛼 and 𝑋 0 = 0, which is the strong solution to (1) [START_REF] Harrison | On skew Brownian motion[END_REF]. We denote by 𝐿 = ℒ(𝑋) its symmetric local time at 0.

For the SBM, ⟨𝑋⟩ 𝑡 = 𝑡, 𝑡 ≥ 0 and 𝐿 = ℒ(|𝑋|) whatever 𝛼 ∈ (0, 1).

Proposition 1 (Scaling property). For any 𝑐 > 0, with

𝑋 𝑐 = 𝑐𝑋 •/𝑐 2 , (𝑋 𝑐 , ℒ(𝑋 𝑐 )) is equal in distribution to (𝑋, 𝐿) when 𝑋 0 = 0. Proof. With (5), ℒ(𝑋 𝑐 ) law = 𝑐𝐿 •/𝑐 2 , so that 𝑋 𝑐 𝑡/𝑐 2 = 𝑐𝐵 𝑡/𝑐 2 + 𝜃𝑐𝐿 𝑡/𝑐 2 and then 𝑋 𝑐 𝑡/𝑐 2 = 𝑊 𝑡 + 𝜃ℒ 𝑡 (𝑋 𝑐 ) for the Brownian motion 𝑊 = 𝑐𝐵 •/𝑐 2 .
The equality in distribution of 𝑋 𝑐 to 𝑋 follows from the uniqueness to the solution to (1).

Hypothesis 1. We consider that 𝑋 0 = 0 so that 𝜉 0 (ℎ) = 0.

This hypothesis aims at simplifying the presentation of the results. Otherwise, we have to wait the random time 𝜏 at which the SBM hits first 0. The observations' timespan is then reduced from [0, 𝑇 ] to [𝜏, 𝑇 ].

For ℎ > 0, define iteratively 𝜏 0 (ℎ) = 0 and

𝜏 𝑘+1 (ℎ) = inf{𝑡 > 𝜏 𝑘 (ℎ); |𝑋 𝑡 -𝑋 𝜏 𝑘 (ℎ) | = ℎ},
the successive passage times and positions of 𝑋 to the grid ℎZ. Also, set 𝜉 𝑘 (ℎ) = 𝑋 𝜏 𝑘 (ℎ) , the normalized passage positions. Since the SBM is a strong Markov process, (𝜉 𝑘 (ℎ)) 𝑘≥0 is a Markov chain with values in ℎZ.

The following result is an easy consequence of the results in [START_REF] Harrison | On skew Brownian motion[END_REF][START_REF] Walsh | A diffusion with discontinuous local time[END_REF] on the scale function and the speed measure of the SBM, or from the construction from K. Itô and H.P. McKean [START_REF] Itô | Diffusion processes and their sample paths[END_REF] through its excursions.

Proposition 2 (Embedding a Markov chain into a SBM). The discrete process (𝜉 𝑘 (ℎ)) 𝑘≥0 is a recurrent Markov chain with values in ℎZ, and (ℎ -1 𝜉 𝑘 (ℎ)) 𝑘≥0 has the same probability transition as the SRW (𝜁 𝑘 ) 𝑘≥0 whose transition kernel is defined by [START_REF] Alvarez | Timing in the Presence of Directional Predictability: Optimal Stopping of Skew Brownian Motion[END_REF].

In addition (𝜏 𝑘+1 (ℎ) -𝜏 𝑘 (ℎ)) 𝑘≥0 is a sequence of independent random variables, each distributed as 𝜏 1 (ℎ). Besides, this sequence is independent from the successive positions (𝜉 𝑘 ) 𝑘≥0 .

Finally, the distribution of 𝜏 1 (ℎ) does not depend on 𝛼 and is the same as the one of ℎ 2 𝜏 1 [START_REF] Lejay | Identification of Diffusion Coefficient in Nonhomogeneous Landscapes[END_REF]. The first two moments of 𝜏 1 [START_REF] Lejay | Identification of Diffusion Coefficient in Nonhomogeneous Landscapes[END_REF]

are E 𝛼 [𝜏 1 (1)] = 1 and E 𝛼 [𝜏 1 (1) 2 ] = 5/3.
Proof. The SBM is a recurrent, strong Markov process so that (𝜉 𝑘 (ℎ)) 𝑘≥0 is a recurrent Markov chain. Besides, if 𝜉 𝑘 (ℎ) ̸ = 0, then the process behaves like a Brownian motion between the times 𝜏 𝑘 (ℎ) and 𝜏 𝑘+1 (ℎ). On the other hand, P 𝛼 [𝜉 𝑘+1 (ℎ) = ℎ|𝜉 𝑘 (ℎ) = 0] is the probability that the first excursion of 𝑋 whose height is greater than 𝑘 is positive, and this probability is equal to 𝛼, an easy consequence from the result from K. Itô and H.P. McKean on the SBM [23, p. 115].

The statement about the distributions of the 𝜏 𝑘+1 (ℎ) -𝜏 𝑘 (ℎ) follows from the strong Markov property. The distribution of (𝜏 𝑘+1 (ℎ) -𝜏 𝑘 (ℎ)) 𝑘≥0 does not depend on 𝛼 as the one of 𝜏 1 (ℎ) depends only on |𝑋| which is the same law as a Reflected Brownian motion. The scaling property is a direct consequence of Proposition 1.

For the moments of 𝜏 1 (1), the Itô formula and the stopping theorem implies that for a Brownian motion 𝐵 (case 𝛼 = 1/2),

1 = E[𝐵 2
𝜏 1 (1) ] = E[𝜏 1 [START_REF] Lejay | Identification of Diffusion Coefficient in Nonhomogeneous Landscapes[END_REF]]. Again from the Itô formula, for any 𝑡 ≥ 0.

𝑡𝐵 2 𝑡 = ∫ 𝑡 0 𝐵 2 𝑠 d𝑠 + ∫ 𝑡 0 2𝑠𝐵 𝑠 d𝐵 𝑠 + ∫ 𝑡 0 𝑠 d𝑠, 𝐵 4 𝑡 = 4 ∫ 𝑡 0 𝐵 3 𝑠 d𝐵 𝑠 + 6 ∫ 𝑡 0 𝐵 2 𝑠 d𝑠.
Again from the stopping theorem,

1 = E[𝐵 4 𝜏 1 (1) ] = 6E [ ∫ 𝜏 1 (1) 0 𝐵 2 𝑠 d𝑠 ] , 1 = E[𝜏 1 (1)𝐵 2 𝜏 1 (1) ] = E [ ∫ 𝜏 1 (1) 0 𝐵 2 𝑠 d𝑠 ] + 1 2 E[𝜏 1 (1) 2 ].
Combining these two equations leads to the second moment of 𝜏 1 [START_REF] Lejay | Identification of Diffusion Coefficient in Nonhomogeneous Landscapes[END_REF].

Generalizing the Donsker theorem, J.M. Harrison and L.A. Shepp [START_REF] Harrison | On skew Brownian motion[END_REF] proved first that the SBM is the limit of a SWR (See also [START_REF] Étoré | On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients[END_REF][START_REF] Gall | One-Dimensional Stochastic Differential Equations Involving the Local Times of the Unknown Process[END_REF]). Hence, the piecewise linear interpolation of (𝜉 𝑘 (ℎ)) 𝑘≥0 along the times (𝜏 𝑘 (ℎ)) 𝑘≥0 converges in probability to 𝑋.

An excursions based estimator

The data. We consider that given a path of a SBM 𝑋, we record the positions on the spatial grid ℎZ of the embedded chain 𝜉(ℎ) = (𝜉 𝑘 (ℎ)) 𝑘=0,...,𝑛 ⋆ (ℎ) . From these observations, we construct an estimator 𝛼(ℎ) of 𝛼 ∈ (0, 1).

Indeed, using the scaling property, we may assume that we record either (𝜉 𝑘 (1)) 𝑘=0,...,𝑇 /ℎ 2 (long time horizon) or 𝜉(ℎ) for 𝜏 𝑘 ≤ 𝑇 (fixed time horizon). We state only our result for the latter.

Embedding the family SRW into a SBM allows one to study the limit of 𝛼(ℎ) as the parameter ℎ decreases to 0. For this, the underlying probability space is the one of the SBM. Of course, our estimator is also suitable when one observes only a SRW. In this case, the estimation in long time horizon has to be considered, and the results are similar in nature (the local time has to be replaced by the limit of 𝑛 -1/2 ∑ 𝑛 𝑘=0 1 𝜉 𝑘 (1)=0 , which has the same distribution as the local time. See [39, Theorem 9.11, p. 99] and also Lemma 3 by noting that 𝐿 𝑛 / √ 𝑛 law = 𝐿 1 and using Remark 6).

The estimator. For 𝑛 ≥ 0, we set

𝑁 𝑛 (ℎ) = 𝑛-1 ∑ 𝑘=0 1 𝜉 𝑘 (ℎ)=0 and 𝑁 + 𝑛 (ℎ) = 𝑛-1 ∑ 𝑘=0 1 𝜉 𝑘 (ℎ)=0, 𝜉 𝑘+1 (ℎ)=1 .
As the random walk

(𝜉 𝑘 (ℎ)) 𝑘≥0 is recurrent, 𝑁 𝑛 (ℎ) → +∞ as 𝑛 → ∞ almost surely. Finally, for 𝑇 > 0, we set 𝑛 ⋆ (ℎ) = sup{𝑘 ≥ 0 | 𝜏 𝑘 (ℎ) ≤ 𝑇 }. From the Strong Law of Large Numbers, lim 𝑛→∞ 𝜏 𝑛 (ℎ)/𝑛 = E 𝛼 [𝜏 1 (ℎ)] = ℎ 2 almost surely. As 𝜏 𝑛 ⋆ (ℎ) (ℎ) ≤ 𝑇 < 𝜏 𝑛 ⋆ (ℎ)+1 (ℎ), ℎ 2 𝑛 ⋆ (ℎ) converges almost surely to 𝑇 as ℎ → 0.
Table 1 summarizes the main variables and their scaling properties, which is justified by Proposition 1.

Our estimator of 𝛼 is

𝛼(ℎ) = 𝑁 + 𝑛 ⋆ (ℎ) (ℎ) 𝑁 𝑛 ⋆ (ℎ) (ℎ) , (6)
which is formalization of (3). ---→

𝑛→∞ 𝑌 law = 𝐿 1 𝑁 + 𝑛 (ℎ) law = 𝑁 + 𝑛 (1) Number of upcrossings in 𝑛 steps 𝑛 ⋆ (ℎ) ≈ 𝑇 ℎ -2
Number of steps of the walk before 𝑇 Table 1: Main variables and their scaling properties.

Remark 1. The likelihood is

Lik 𝛼 (𝜉(ℎ)) = 𝑛 ⋆ (ℎ)-1 ∏ 𝑘=0 𝑝 𝛼 (𝜉 𝑘 (ℎ), 𝜉 𝑘+1 (ℎ)) = 𝛼 𝑁 + 𝑛 ⋆ (ℎ) (ℎ) (1 -𝛼) 𝑁 𝑛 ⋆ (ℎ) (ℎ)-𝑁 + 𝑛 ⋆ (ℎ) (ℎ)
and the log-likelihood is

Log-Lik 𝛼 (𝜉(ℎ)) = 𝑁 + 𝑛 ⋆ (ℎ) (ℎ) log(𝛼) + (𝑁 𝑛 ⋆ (ℎ) (ℎ) -𝑁 + 𝑛 ⋆ (ℎ) (ℎ)) log(1 -𝛼).
Solving 𝜕 𝛼 Log-Lik 𝛼(ℎ) (𝜉(ℎ)) = 0 implies that 𝛼(ℎ) given by ( 6) is the maximum likelihood estimator (MLE)

𝛼(ℎ) = argmax 𝛼 Log-Lik 𝛼 (𝜉) = 𝑁 + 𝑛 ⋆ (ℎ) (ℎ) 𝑁 𝑛 ⋆ (ℎ) (ℎ)
.

Theorem 1 (Convergence of the estimator 𝛼(ℎ)). For any 𝛼 ∈ (0, 1), under P 𝛼 , 𝛼(ℎ)

a.s. --→ ℎ→0 𝛼, (7)
implying in particular that 𝛼(ℎ) is a consistent estimator of 𝛼. In addition, there exists a Brownian motion 𝑊 independent from 𝐵 such that

1 √ ℎ (𝛼(ℎ) -𝛼) dist. --→ ℎ→0 Υ 𝑇 = √ 𝛼(1 -𝛼) 𝑊 (𝐿 𝑇 ) 𝐿 𝑇 .
Moreover, for some constant 𝐶,

sup 𝑥∈R ⃒ ⃒ P[ℎ -1/2 (𝛼 ℎ -𝛼) ≤ 𝑥] -P[Υ 𝑇 ≤ 𝑥] ⃒ ⃒ ≤ 𝐶ℎ 2 9 .
Remark 2. Using the scaling property on both the Brownian motion and the local time, Υ 𝑇 law = 𝑇 -1/4 Υ 1 for any 𝑇 > 0. See [START_REF] Lejay | Is a Brownian motion skew?[END_REF] for more results about this random variable. Remark 3. Since our estimator relies on the convergence of binomial random variables toward a Gaussian one, one could make use of the large literature on confidence intervals, even for small samples, which are better than the Wald confidence interval based on the Gaussian approximation: See e.g., [START_REF] Blyth | Approximate binomial confidence limits[END_REF][START_REF] Brown | Confidence intervals for a binomial proportion and asymptotic expansions[END_REF][START_REF] Ross | Accurate confidence intervals for binomial proportion and Poisson rate estimation[END_REF]. In particular, for 𝛼 ̸ = 0, these specific approaches provide non symmetric confidence intervals around the estimator for a finite number of samples, although the limiting distribution of ℎ -1/2 (𝛼(ℎ) -𝛼) is symmetric. This was observed empirically on the numerical tests in [START_REF] Lejay | Simulating diffusion processes in discontinuous media: a numerical scheme with constant time steps[END_REF] on the MLE for the SBM.

Remark 4. While we record around 𝑛 ⋆ (ℎ) ≈ 𝑇 /ℎ 2 points of the SBM, we really use 𝑁 𝑛 ⋆ (ℎ) (ℎ) points for the estimation, as the 𝜉 𝑘 (ℎ)'s bring no information when the walk is away from 0. The number of points really useful for the estimation is roughly of order √ 𝑛 ⋆ (ℎ) ≈ √ 𝑇 /ℎ. Therefore, the rate of convergence in our CLT like result is ℎ 1/2 , that is of order 𝑛 ⋆ (ℎ) -1/4 . We are then not in the "usual" situation where the rate of convergence is 𝑛 ⋆ (ℎ) -1/2 . These are the same phenomena as the one of D. Florens [START_REF] Florens | Estimation of the diffusion coefficient from crossings[END_REF] for estimators for stochastic differential equations based on crossings, and the results of J. Jacod [START_REF] Jacod | Rates of convergence to the local time of a diffusion[END_REF].

The convergence of 𝛼(ℎ) in [START_REF] Blyth | Approximate binomial confidence limits[END_REF] depends on the next lemma and the almost sure convergence of 𝑛 ⋆ (ℎ) to infinity.

Lemma 1 (Distribution of 𝑁 + 𝑛 (ℎ) given 𝑁 𝑛 (ℎ)). Under P 𝛼 , the distribution of 𝑁 + 𝑛 (ℎ) given 𝑁 𝑛 (ℎ) is a binomial one with parameters (𝑁 𝑛 (ℎ), 𝛼), so that

𝛼(ℎ) = 𝑁 + 𝑛 (ℎ) 𝑁 𝑛 (ℎ) ---→ 𝑛→∞ 𝛼 almost surely. (8)
The convergence of ℎ -1/2 (𝛼(ℎ)-𝛼) follows from the Berry-Esseen theorem, whose application is immediate because

P[(𝑁 + 𝑛 ⋆ (ℎ) (ℎ), 𝑁 𝑛 ⋆ (ℎ) (ℎ)) = (𝑘, ℓ)] = ∑ 𝑛≥0 E[P[𝑁 + 𝑛 (ℎ) = 𝑘 | 𝑁 𝑛 (ℎ)] 1 𝑁𝑛(ℎ)=ℓ ]P[𝑛 ⋆ (ℎ) = 𝑛]
for any integers 𝑘, ℓ ≥ 0.

In a first time, we focus of the asymptotic behavior of the distribution of 𝑁 + 𝑛 (ℎ) given 𝑁 𝑛 (ℎ) as 𝑁 𝑛 (ℎ) becomes large.

Lemma 2 (Application of the Berry-Esseen inequality [START_REF] Berry | The accuracy of the Gaussian approximation to the sum of independent variates[END_REF][START_REF] Esseen | On the Liapounoff limit of error in the theory of probability[END_REF]). Let Φ be the cumulative distribution of the Gaussian distribution 𝒩 (0, 1).

With 𝜎 = √ 𝛼(1 -𝛼), sup 𝑥∈R ⃒ ⃒ ⃒ ⃒ ⃒ P 𝛼 [ 𝑁 + 𝑛 (ℎ) -𝑁 𝑛 (ℎ)𝛼 𝜎 √ 𝑁 𝑛 (ℎ) ≤ 𝑥 𝑁 𝑛 (ℎ) ] -Φ(𝑥) ⃒ ⃒ ⃒ ⃒ ⃒ ≤ 𝐶𝛼 𝜎 3 √ 𝑁 𝑛 (ℎ) , ( 9
)
where 𝐶 is a constant.

Remark 5. Many articles have been devoted to exhibit the best constant 𝐶 in ( 9) for binomial distribution (See e.g. [START_REF] Nagaev | On an estimate for the closeness of the binomial distribution to the normal distribution[END_REF][START_REF] Nagaev | On the estimation of the closeness of the binomial distribution to the normal distribution[END_REF] and the references within).

From now, the result stems from two ingredients: (i) the convergence of ℎ𝑁 𝑛 ⋆ (ℎ) toward 𝐿 𝑇 , which is granted by the P. Lévy's downcrossing theorem (See e.g. [START_REF] Chung | Downcrossings and local time. Z. Wahrscheinlichkeitstheorie und Verw[END_REF][START_REF] Itô | Diffusion processes and their sample paths[END_REF][START_REF] Khoshnevisan | Exact rates of convergence to Brownian local time[END_REF] or [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF]Exercice VI.1.19] and references within for refinements of the statement). (ii) A control on E 𝛼 [1/ √ 𝑁 𝑛 ⋆ (ℎ) (ℎ)] which follows from the convergence of 𝑁 𝑛 (ℎ)/ √ 𝑛 (more precise statements are given in [START_REF] Chung | On the zeros of ∑ 𝑛 1 ±1[END_REF] which contains some almost sure control of 𝑁 𝑛 (ℎ)/ √ 𝑛 for 𝑛 large enough).

Lemma 3 (P. Lévy downcrossing theorem). Almost surely, ℎ𝑁 𝑛 ⋆ (ℎ) (ℎ) converges to 𝐿 𝑇 as ℎ → 0.

Besides, for some constant 𝐶,

E 𝛼 [ ⃒ ⃒ ⃒ ⃒ 𝑁 𝑛 (1) √ 𝑛 - 𝐿 𝑛 √ 𝑛 ⃒ ⃒ ⃒ ⃒ 2 ] 1/2 ≤ 𝐶 𝑛 1/4 , ∀𝑛 ≥ 1. (10)
Remark 6. The cumulative distribution function of the local time 𝐿

1 is P[𝐿 1 ≤ 𝑥] = √ 2/𝜋 ∫ 𝑥 0 exp(-𝑦 2 /2) d𝑦 for 𝑥 ≥ 0.
(See e.g. [39, Theorem 9.1, p. 101]), a consequence of the reflection principle and the identity in distribution between 𝐿 1 and the supremum of a Brownian motion on [0, 1].

Remark 7. It is proved in [START_REF] Csáki | A combinatorial proof of a theorem of P. Lévy on the local time[END_REF] (see also [START_REF] Révész | Random walk in random and nonrandom environments[END_REF]Consequence 10.3,p. 110]) that for any 𝜖 > 0, 𝑁 𝑛 (1) -𝐿 𝑛 = o(𝑛 1/4+𝜖 ) (See [START_REF] Bass | Rates of convergence to Brownian local time[END_REF][START_REF] Csörgő | On best possible approximations of local time[END_REF] for more precise bounds). Thus, the rate of convergence of the L 2 -norm of 1/𝑛 1/4 in ( 10) is "nearly" optimal.

Proof. The first statement is the P. Lévy's downcrossing theorem, so that we skip its proof. It is valid for the SBM as it depends only on the absolute value of 𝑋, which is a reflected Brownian motion in distribution.

We consider that ℎ = 1. Thus, we drop any reference to ℎ.

We write sgn(𝑥) = 1 if 𝑥 > 0, sgn(𝑥) = -1 if 𝑥 < 0 and sgn(0) = 0 since we use the symmetric local time. The Itô-Tanaka formula yields that

|𝑋 𝑡 | = |𝑋 𝑠 | + ∫ 𝑡 𝑠 sgn(𝑋 𝑟 ) d𝑋 𝑟 + 𝐿 𝑡 -𝐿 𝑠 for 0 ≤ 𝑠 ≤ 𝑡. Since |𝑎| -|𝑏| ≤ |𝑎 -𝑏| for 𝑎, 𝑏 ∈ R, |𝐿 𝑛 -𝐿 𝜏𝑛 | ≤ |𝑋 𝑛 -𝑋 𝜏𝑛 | + ⃒ ⃒ ⃒ ⃒ ∫ 𝑛 𝜏𝑛 sgn(𝑋 𝑟 ) d𝑋 𝑟 ⃒ ⃒ ⃒ ⃒ for 𝑛 ≥ 0. As ∫ 𝑡 𝑠 sgn(𝑋 𝑟 ) d𝑋 𝑟 = ∫ 𝑡 𝑠 sgn(𝑋 𝑟 ) d𝐵 𝑟 since sgn(0) = 0, E 𝛼 [|𝐿 𝑛 -𝐿 𝜏𝑛 | 2 ] ≤ 2E 𝛼 [|𝑋 𝑛 -𝑋 𝜏𝑛 | 2 ] + 2E 𝛼 [|𝜏 𝑛 -𝑛|]. Since (𝑋 𝑡 -𝑋 𝑟 ) 2 = 2 ∫ 𝑡 𝑠 𝑋 𝑟 d𝑋 𝑟 + (𝑡 -𝑠) and ∫ 𝑡 𝑠 𝑋 𝑟 d𝐿 𝑟 = ∫ 𝑡 𝑠 𝑋 𝑟 1 𝑋𝑟=0 d𝐿 𝑟 = 0, 1 𝑛 E 𝛼 [|𝐿 𝑛 -𝐿 𝜏𝑛 | 2 ] ≤ 4E 𝛼 [⃒ ⃒ ⃒ 𝜏 𝑛 𝑛 -1 ⃒ ⃒ ⃒ ] ≤ 4 √ 2 √ 3𝑛 , (11)
as 𝜏 𝑛 is the sum of 𝑛 independent random variables of mean 1 and variance 2/3. On the other hand, it is easily seen that for any 𝑘 ≥ 0,

|𝜉 𝑘+1 | = |𝜉 𝑘 | + sgn(𝜉 𝑘 )(𝜉 𝑘+1 -𝜉 𝑘 ) + 1 𝜉 𝑘 =0 . (12)
Comparing ( 2) with ( 12), as 𝜉 𝑘 = 𝑋 𝜏 𝑘 ,

𝑁 𝑛 = 𝐿 𝜏𝑛 + 𝑀 𝑛 with 𝑀 𝑛 = ∫ 𝜏𝑛 0 (sgn(𝑋 𝑟 ) -sgn(𝑋 𝜋(𝑟) )) d𝑋 𝑟 ,
where 𝜋(𝑠) = 𝜏 𝑘 when 𝜏 𝑘 ≤ 𝑠 < 𝜏 𝑘+1 . It follows that

1 𝑛 E 𝛼 [𝑀 2 𝑛 ] = 1 𝑛 E 𝛼 [∫ 𝜏𝑛 0 (sgn(𝑋 𝑟 ) -sgn(𝑋 𝜋(𝑟) )) 2 d𝑟 ] ≤ 1 𝑛 E 𝛼 [∫ 𝜏𝑛 0 1 𝑋𝑟∈[-1,1] d𝑟 ] ≤ 1 𝑛 E 𝛼 [∫ 𝑛 0 1 𝑋𝑟∈[-1,1] d𝑟 ] + 1 𝑛 E 𝛼 [|𝜏 𝑛 -𝑛|]. Since |𝑋| law = |𝐵|, P 𝛼 [𝑋 𝑟 ∈ [-1, 1]] = P 𝛼 [|𝐵 𝑟 | ∈ [-1, 1]] ≤ 2/ √ 2𝜋𝑟, so that 1 𝑛 E 𝛼 [𝑀 2 𝑛 ] ≤ 1 𝑛 + 4 √ 2𝜋𝑛 . (13)
Eq. ( 10) is then proved by combining [START_REF] Chung | On the zeros of ∑ 𝑛 1 ±1[END_REF] with [START_REF] Csáki | A combinatorial proof of a theorem of P. Lévy on the local time[END_REF].

Lemma 4. For a constant 𝐶,

E 𝛼 [ 1 √ 𝑁 𝑛 ⋆ (ℎ) (ℎ) ] ≤ 𝐶ℎ 2 9 --→ ℎ→0 0. (14)
Proof. Since 𝑁 𝑛 (ℎ) ≥ 1, E 𝛼 [𝑁 𝑛 (ℎ) -1/2 ] ≤ 1 for any 𝑛 ∈ N, ℎ > 0.

From the scaling property, 𝐿 𝑛 / √ 𝑛 law = 𝐿 1 , and the law of 𝐿 1 has the density 2 exp(-𝑦 2 /2)/ √ 2𝜋 (See Remark 6). We write 𝑁 𝑛 for 𝑁 𝑛 [START_REF] Lejay | Identification of Diffusion Coefficient in Nonhomogeneous Landscapes[END_REF]. It holds that

P 𝛼 [ 𝑁 𝑛 √ 𝑛 ≤ 𝑥 ] ≤ P 𝛼 [ 𝑁 𝑛 √ 𝑛 ≤ 𝑥; ⃒ ⃒ ⃒ ⃒ 𝑁 𝑛 √ 𝑛 - 𝐿 𝑛 √ 𝑛 ⃒ ⃒ ⃒ ⃒ ≤ 𝛿 ] +P 𝛼 [ 𝑁 𝑛 √ 𝑛 ≤ 𝑥; ⃒ ⃒ ⃒ ⃒ 𝑁 𝑛 √ 𝑛 - 𝐿 𝑛 √ 𝑛 ⃒ ⃒ ⃒ ⃒ ≥ 𝛿 ] ≤ P 𝛼 [ 𝐿 𝑛 √ 𝑛 ≤ 𝑥 + 𝛿 ] + 1 𝛿 E 𝛼 [ ⃒ ⃒ ⃒ ⃒ 𝑁 𝑛 √ 𝑛 - 𝐿 𝑛 √ 𝑛 ⃒ ⃒ ⃒ ⃒ 2 ] 1/2 ≤ √ 2 𝜋 (𝑥 + 𝛿) + 𝐶 𝛿𝑛 1/4 . (15) Let 𝜑 : N → R + such that 𝜑(𝑛) ↘ 𝑛→∞ 0 and 𝑛 1/4 𝜑(𝑛) ↗ 𝑛→∞ +∞. Since 𝑁 𝑛 ≥ 1, E 𝛼 [ 1 √ 𝑁 𝑛 ] = E 𝛼 [ 1 √ 𝑁 𝑛 ; 𝑁 𝑛 < √ 𝑛𝜑(𝑛) ] + E 𝛼 [ 1 √ 𝑁 𝑛 ; 𝑁 𝑛 ≥ √ 𝑛𝜑(𝑛) ] ≤ P 𝛼 [𝑁 𝑛 < √ 𝑛𝜑(𝑛)] + 1 √ 𝑛 1/2 𝜑(𝑛) .
Choosing 𝑥 = 𝛿 = 𝜑(𝑛) in [START_REF] Delattre | Asymptotic equivalence for a null recurrent diffusion[END_REF], for some constant 𝐶,

P 𝛼 [ 𝑁 𝑛 √ 𝑛 ≤ 𝜑(𝑛) ] ≤ 𝐶𝜑(𝑛) + 𝐶 𝑛 1/4 𝜑(𝑛)
.

By balancing the expressions, 𝜑(𝑛) = 1/𝑛 1/8 is a suitable choice which ensures that

E 𝛼 [ 1 √ 𝑁 𝑛 ] ≤ 𝐶 𝑛 1/8 ---→ 𝑛→∞ 0. (16)
For 𝜅(ℎ) ↗ ℎ→0 +∞ and ℎ 2 𝜅(ℎ) ↘ ℎ→0 0, since 𝑁 𝑛 is increasing with 𝑛 and greater than 1,

E 𝛼 [ 1 √ 𝑁 𝑛 ⋆ (ℎ) (ℎ) ] ≤ E 𝛼 [ 1 √ 𝑁 𝑛 ⋆ (ℎ) (ℎ) ; 𝑛 ⋆ (ℎ) ≤ 𝜅(ℎ) ] + E 𝛼 [ 1 √ 𝑁 𝜅(ℎ) (ℎ) ] ≤ P 𝛼 [𝑛 ⋆ (ℎ) ≤ 𝜅(ℎ)] + E 𝛼 [ 1 √ 𝑁 𝜅(ℎ) (1)
] .

Using the Bienaymé-Tchebychev inequality,

P 𝛼 [𝑛 ⋆ (ℎ) ≤ 𝜅(ℎ)] ≤ P 𝛼 [𝜏 𝜅(ℎ)+1 (ℎ) ≥ 𝑇 ] ≤ 1 𝑇 E 𝛼 [𝜏 𝜅(ℎ)+1 (ℎ)] ≤ 𝜅(ℎ) + 1 𝑇 ℎ 2 .
Using ( 16), we take 𝜅(ℎ) = ℎ -16/9 to balance the rate of convergence in the two terms in the right-hand side of ( 17). This gives [START_REF] Csörgő | On best possible approximations of local time[END_REF].

Proof of Theorem 1. Since 𝑋 0 = 0, 𝑁 𝑛 (ℎ) > 0 and 𝐿 𝑇 > 0 almost surely. The convergence of 𝛼(ℎ

) = 𝑁 + 𝑛 ⋆ (ℎ) (ℎ)/𝑁 𝑛 ⋆ (ℎ) (ℎ) to 𝛼 is a direct consequence of Lemma 1 and the convergence of ℎ -2 𝑛 ⋆ (ℎ) to 𝑇 , since 𝑛 ⋆ (ℎ) is independent from {(𝑁 𝑛 (ℎ), 𝑁 + 𝑛 (ℎ))} 𝑛≥0 . From the Berry-Esseen theorem (See Lemma 2), 𝐹 𝑛 (𝑥) = P 𝛼 [ 𝑇 1/4 √ ℎ ( 𝑁 + 𝑛 (ℎ) 𝑁 𝑛 (ℎ) -𝛼 ) ≤ 𝑥 ] = E 𝛼 [ P 𝛼 [ 𝑁 + 𝑛 (ℎ) -𝛼𝑁 𝑛 (ℎ) 𝜎 √ 𝑁 𝑛 (ℎ) ≤ 𝑥 √ ℎ𝑁 𝑛 (ℎ) 𝜎𝑇 1/4 𝑁 𝑛 (ℎ) ]] = E 𝛼 [ Φ ( 𝑥 √ ℎ𝑁 𝑛 (ℎ) 𝜎𝑇 1/4 )] + 𝜖 𝑛 (ℎ) with |𝜖 𝑛 (ℎ)| ≤ E 𝛼 [ 𝐶𝛼 𝜎 3 √ 𝑁 𝑛 (ℎ)
] .

Hence,

E 𝛼 [𝐹 𝑛 ⋆ (ℎ) (𝑥)] = E 𝛼 [ Φ ( 𝑥 √ ℎ𝑁 𝑛 ⋆ (ℎ) (ℎ) 𝜎𝑇 1/4 )] + E 𝛼 [𝜖 𝑛 ⋆ (ℎ) (ℎ)].
The convergence of ℎ𝑁 𝑛 ⋆ (ℎ) (ℎ) to 𝐿 𝑇 and Lemma 4 imply the convergence in distribution of ℎ -1/2 (𝛼(ℎ) -𝛼) to 𝜎𝐺/ √ 𝐿 𝑇 for 𝐺 ∼ 𝒩 (0, 1), where 𝐺 is independent from 𝐿 𝑇 . The rate of convergence of the cumulative distribution function is given by Lemma 4.

Let us consider a Brownian motion 𝑊 independent from 𝐵. Using the scaling property of the Brownian motion,

E 𝛼 [ Φ ( 𝑥 √ 𝐿 𝑇 𝜎𝑇 1/4 )] = P 𝛼 [ 𝜎𝑇 1/4 𝑊 (1) √ 𝐿 𝑇 ≤ 𝑥 ] = P 𝛼 [ 𝜎𝑇 1/4 𝑊 (𝐿 𝑇 ) 𝐿 𝑇 ≤ 𝑥 ] .
Thus, the limit in distribution of ℎ -1/2 (𝛼(ℎ)-𝛼) could be written 𝜎𝑊 (𝐿 𝑇 )/𝐿 𝑇 .

Remark 8. The computations above are easily modified to show that for any 𝑎, 𝑏 ∈ R,

𝑎ℎ𝑁 𝑛 ⋆ (ℎ) (ℎ) + 𝑏 𝑇 1/4 √ ℎ ( 𝑁 + 𝑛 ⋆ (ℎ) (ℎ) 𝑁 𝑛 ⋆ (ℎ) (ℎ) -𝛼 ) law --→ ℎ→0 𝑎𝐿 𝑇 + 𝑏 𝜎𝑇 1/4 𝑊 (𝐿 𝑇 ) 𝐿 𝑇 .
The Cramér-Wold device proves the joint convergence in distribution of (ℎ𝑁 𝑛 ⋆ (ℎ) , ℎ -1/2 (𝛼(ℎ) -𝛼)) to (𝐿 𝑇 , 𝜎𝑊 (𝐿 𝑇 )/𝐿 𝑇 ). Since 𝛼(ℎ) converges almost surely to 𝛼 which is deterministic, the vector (ℎ𝑁 𝑛 ⋆ (ℎ) (ℎ), ℎ -1/2 (𝛼(ℎ) -𝛼), 𝛼(ℎ)) converges in distribution to (𝐿 𝑇 , 𝜎𝑊 (𝐿 𝑇 )/𝐿 𝑇 , 𝛼).

The Local Asymptotic Mixed Normal property

The Local Asymptotic Mixed Normal (LAMN) property [START_REF] Jeganathan | On the asymptotic theory of estimation when the limit of the log-likelihood ratios is mixed normal[END_REF] extends the Local Asymptotic Normal property introduced by L. Le Cam [START_REF] Ibragimov | Has'minskiȋ. Statistical estimation[END_REF][START_REF] Lecam | On some asymptotic properties of maximum likelihood estimates and related Bayes' estimates[END_REF] for limits which are mixed normal distributions. Among other properties, one may deduce results on the convergence of estimators constructed from 𝜉(ℎ) (See e.g. [START_REF] Jeganathan | Some asymptotic properties of risk functions when the limit of the experiment is mixed normal[END_REF]). for some constant 𝐶. Since ℎ𝑁 (ℎ) converges almost surely to 𝐿 𝑇 which is almost surely finite, 𝜖(ℎ, 𝛾) converges almost surely to 0 for any 𝛾 > 0, hence in probability. Thanks to Lemma 1,

𝐴(ℎ) = 1 𝑁 (ℎ) ( 𝑁 + (ℎ) 𝛼 - 𝑁 -(ℎ) 1 -𝛼 ) = 1 𝜎 2 ( 𝑁 + (ℎ) 𝑁 (ℎ) -𝛼 )
a.s.

--→ ℎ→0 0 and 𝐵(ℎ) = 1 𝑁 (ℎ)

( 𝑁 + (ℎ) 𝛼 2 + 𝑁 -(ℎ) (1 -𝛼) 2
) a.s.

--→ --→ ℎ→0 𝑊 (1) ∼ 𝒩 (0, 1).

Using again Remark 8, we obtain indeed the joint convergence of (𝑄(ℎ), 𝐺(ℎ)) in distribution. This proves the LAMN property.

𝑋 𝑡 law =

 law ℎ𝑋 𝑡/ℎ 2 SBM, same scaling of for Brownian motion𝐿 𝑡 law = ℎ𝐿 𝑡/ℎ 2Local time at 0, same scaling as for Brownian local time𝜏 𝑘 (ℎ) law = ℎ 2 𝜏 𝑘 (1) Sum of iid rv, E 𝛼 [𝜏 1 (1)] = 1, E 𝛼 [𝜏 1 (1) 2 ] = 5/3 𝜉 𝑘 (ℎ) law = ℎ𝜉 𝑘 (1) Biased random walk 𝑁 𝑛 (ℎ) law = 𝑁 𝑛 (1)Number of zeros in 𝑛 steps, 𝑁 𝑛 (ℎ)/ √ 𝑛 dist.

Proposition 3 (

 3 The LAMN property). For any 𝛾 ∈ R and ℎ > 0,log Lik 𝛼+𝛾 √ ℎ (𝜉(ℎ)) Lik 𝛼 (𝜉(ℎ)) -𝛾 √ 𝑄(ℎ)𝐺(ℎ) -𝛾 2 2 𝑄(ℎ) + 𝜖(ℎ, 𝛾) P𝛼 --→ ℎ→0 0where 𝐺(ℎ) and 𝑄(ℎ) are ℱ 𝑛 ⋆ (ℎ) -measurable,(𝐺(ℎ), 𝑄(ℎ)) 𝜎 -2 𝐿 𝑇 ) with 𝐻 ∼ 𝒩 (0, 1)and 𝐻 and 𝐿 𝑇 are independent. In other words, (P 𝛼 ∘ 𝜉(ℎ) -1 ) ℎ>0 satisfies the LAMN property.Proof. To simplify the notations, let us set𝑁 (ℎ) = 𝑁 𝑛 ⋆ (ℎ) , 𝑁 + (ℎ) = 𝑁 𝑛 ⋆ (ℎ) and 𝑁 -(ℎ) = 𝑁 (ℎ) -𝑁 + (ℎ).The logarithm of the ratio of the likelihood is, for any 𝛾 ∈ R, the Taylor development of the logarithm, 𝜖(ℎ, 𝛾) ≤ 𝐶𝛾 3 𝑁 (ℎ)ℎ 3/2

  𝑇 ) 𝜎𝐿 𝑇 .On the other hand, since ℎ𝑁 (ℎ) converges almost surely to 𝐿 𝑇 and owing to Remark 8,√ ℎ𝑁 (ℎ)𝐴(ℎ) = ℎ𝑁 (ℎ) ℎ𝑁 (ℎ)𝐴(ℎ) -𝛾 2 ℎ𝑁 (ℎ)𝐵(ℎ) + 𝜖(ℎ, 𝛾).Besides, 𝑄(ℎ) = ℎ𝑁 (ℎ)𝐵(ℎ)

The reader shall be warned that different normalizations are used in different textbooks. For example, in[40, Chapter VI], the right local time is used. For the SBM, this distinction is crucial since 𝑎 ↦ → 𝐿 𝑎 𝑡 (𝑋) is discontinuous at 0[START_REF] Walsh | A diffusion with discontinuous local time[END_REF].

Acknowledgement. This work has been developed within the framework of the Inria's Équipe Associée ANESTOC-TOSCA between France and Chile. It is associated to a joint work with E. Mordecki and S. Torres on the estimation of the parameter of the Skew Brownian motion. The author wishes to thank them for interesting discussion on this topic.

The author also wishes the referees for their careful reading and having suggested corrections, improvements and for having pointing out the attention toward Bayesian estimation.


Olivier Beaumont

Lionel Eyraud-Dubois

Thomas Lambert

Cuboid Partitioning for Parallel Matrix Multiplication on Heterogeneous Platforms

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In the case of homogeneous resources, the problem of partitioning data for Linear Algebra kernels in order to both balance the load throughout the computation and to minimize communications is well understood. 2D block-cyclic distributions, for instance, have been introduced in Scalapack [START_REF] Choi | ScaLAPACK: A Portable Linear Algebra Library for Distributed Memory Computers: Design Issues and Performance[END_REF] in order to achieve this goal. More recently, the problem has received a lot of attention in Communication Avoiding algorithms design (see [START_REF] Hoemmen | Communication-avoiding Krylov Subspace Methods[END_REF][START_REF] Anderson | Communication-avoiding QR Decomposition for GPUs[END_REF] and [START_REF] Solomonik | Communication-optimal Parallel 2.5 D Matrix Multiplication and LU factorization Algorithms[END_REF][START_REF] Ballard | Communicationoptimal Parallel Algorithm for Strassen's Matrix Multiplication[END_REF] for Matrix Multiplication specifically). In this context, the goal is to partition the set of elementary computations to be performed into a minimal number of zones, each zone being able to be processed in local memory (i.e. both input, intermediate and output data). This corresponds to maximizing the volume of computations that can be processed with a given amount of memory.

In this paper, we concentrate on Matrix Multiplication algorithms and more specifically on Matrix Multiplication algorithms that involve N 3 elementary operations of type C i,j ← C i,j + A i,k B k,j , i.e. we ignore variants such as Strassen or Coppersmith-Winograd. Note that throughout the paper, we will assume that matrices are partitioned into blocks, whose size is chosen so as to be well adapted to all types of resources (typically CPUs and GPUs). On the other hand, we consider a fully heterogeneous platform, where all nodes may have different processing capacities and we address the most general problem, where several partially aggregated copies of C can reside simultaneously in memory, such as in 2.5D algorithms [START_REF] Solomonik | Communication-optimal Parallel 2.5 D Matrix Multiplication and LU factorization Algorithms[END_REF]. In this context, the problem consists in partitioning the computational domain (the cube of N 3 points) into sub-domains allocated to the different resources. In order to balance the load between the processing units, each unit should receive a volume of computations proportional to its processing speed and the overall amount of communications, that corresponds to the overall boundary area between the zones should be minimized.

Many algorithms [START_REF] Kalinov | Heterogeneous Distribution of Computations Solving Linear Algebra Problems on Networks of Heterogeneous Computers[END_REF][START_REF] Beaumont | Partitioning a Square into Rectangles: NP-completeness and Approximation Algorithms[END_REF][START_REF] Becker | Towards Data Partitioning for Parallel Computing on Three Interconnected Clusters[END_REF][START_REF] Deflumere | Optimal Data Partitioning Shape for Matrix Multiplication on Three Fully Connected Heterogeneous Processors[END_REF][START_REF] Nagamochi | An Approximation Algorithm for Dissecting a Rectangle into Rectangles with Specified Areas[END_REF][START_REF] Fügenschuh | Exact and Approximation Algorithms for a Soft Rectangle Packing Problem[END_REF] have been proposed in the context of dense matrix multiplication based on Canon's-like algorithm, that corresponds to the 2D version of the problem, i.e. how to partition a matrix into zones of fixed area while minimizing the overall length of the boundaries. On the other hand, to the best of our knowledge, this paper is the first to consider the complexity of the 3D version of the algorithm, to prove the NP-Completeness of the underlying decision problem and to propose an approximation algorithm for it.

Related Works

The 2D version of this optimization problem has been first introduced by Lastovetsky and Kalinov in [START_REF] Kalinov | Heterogeneous Distribution of Computations Solving Linear Algebra Problems on Networks of Heterogeneous Computers[END_REF]. In [START_REF] Beaumont | Partitioning a Square into Rectangles: NP-completeness and Approximation Algorithms[END_REF], it has been proven that the problem is NP-Complete, and a first algorithm with bounded approximation ratio (1.75) has been proposed. This algorithm has been improved along two directions. On the one hand, Lastovetsky et al. have proposed to relax the assumption stating that the zones allocated to the processors should be single rectangles and have proposed optimal algorithms, but limited to 2 heterogeneous processors [START_REF] Becker | Towards Data Partitioning for Parallel Computing on Three Interconnected Clusters[END_REF] and more recently to 3 heterogeneous processors [START_REF] Deflumere | Optimal Data Partitioning Shape for Matrix Multiplication on Three Fully Connected Heterogeneous Processors[END_REF]. On the other hand, recursive partitioning algorithms have recently been proposed where at each step, the set of processors is split into two parts. Sophisticated proof techniques enabled Nagamochi and Abe [START_REF] Nagamochi | An Approximation Algorithm for Dissecting a Rectangle into Rectangles with Specified Areas[END_REF] to improve the approximation ratio down to 1.25. Recently, Fügenschuh et al.. [START_REF] Fügenschuh | Exact and Approximation Algorithms for a Soft Rectangle Packing Problem[END_REF] improved this result to 1.15, but under the assumption that if we consider processors in decreasing order of their processing speeds, there is no abrupt change in the performance between 2 successive processors. Unfortunately, such an abrupt decrease typically happens when considering nodes consisting of CPUs and GPUs, such that Fügenschuh's algorithm is limited to the case of relatively homogeneous platforms. In [START_REF] Beaumont | A New Approximation Algorithm for Matrix Partitioning in Presence of Strongly Heterogeneous Processors[END_REF], an algorithm based on the idea of non rectangular partitioning proposed by Lastovetsky and extended to any number of processors by adapting the recursive partitioning algorithm proposed by Nagamochi has been proposed. It achieves an approximation ratio of 1.15 and does not require any specific assumption on the relative speed of resources, so that it can be used in the case of nodes consisting of both regular cores and accelerators.

Besides a single heterogeneous node, this partitioning problem has been adapted to distributed hierarchical and highly heterogeneous platforms in [START_REF] Clarke | Hierarchical Partitioning Algorithm for Ccientific Computing on Highly Heterogeneous CPU + GPU Clusters[END_REF], where the partitioning is applied at two levels (intra-node and inter-node), based on sophisticated performance models. The same partitioning has also been extended to finite-difference time-domain (FDTD) methods to obtain numerical solutions of Maxwell's equations in [START_REF] Shams | On Optimization of Finite-difference Time-domain (FDTD) Computation on Heterogeneous CPU and GPU clusters[END_REF]. More dynamic settings have also been considered in [START_REF] Mohamed | DDOps: Dual-Direction Operations for Load Balancing on Non-dedicated Heterogeneous Distributed Systems[END_REF]. Recently, in order to cope with resource heterogeneity and the difficulty to build optimal schedules, the use of dynamic runtime schedulers have been proposed, such as StarPU [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF], StarSs [START_REF] Planas | Hierarchical Task-based Programming with StarSs[END_REF], or PaRSEC [START_REF] Bosilca | PaRSEC: A Programming Paradigm Exploiting Heterogeneity for Enhancing Scalability[END_REF]. In these systems, at runtime, the scheduler takes the scheduling and allocation decisions based on the set of ready tasks (tasks whose all data and control dependencies have been resolved), on the availability of the resources (estimated using expecting processing and communication times), and on the actual location of input data. The comparison between static scheduling strategies and runtime scheduling strategies has been considered in [START_REF] Beaumont | Comparison of Static and Dynamic Resource Allocation Strategies for Matrix Multiplication[END_REF], where the analysis of the behavior of static, dynamic, and hybrid strategies highlights the benefits of introducing more static knowledge and allocation decisions in runtime libraries.

Paper Outline

The paper is organized as follows. In Section 2, we formally present the partitioning problem and the notations that will be used throughout the paper. In Section 3, the complexity of the associated decision problem in the 3D case is established and a 5 6 2/3 1.51-approximation algorithm for cube-partitioning is proposed in Section 4. Conclusions and perspectives are given in Section 5.

General Context

Definition 1. Let P be a connected polyhedron included in [0, x] × [0, y] × [0, z].
We define its covering cuboid as the smallest cuboid

Cu(P ) = [x 1 , x 2 ] × [y 1 , y 2 ] × [z 1 , z 2 ]
that includes P . We define also its width w(P ) as x 2 -x 1 , its height h(P ) as y 2 -y 1 and its length l(P ) as z 2 -z 1 . Let us define Hs(P ) = h(P )l(P ) + w(P )l(P )+h(P )w(P ), ρ(P ) = max(h(P ),w(P ),l(P ))

min(h(P ),w(P ),l(P )) and ρ (P ) = max(h(P ),w(P ),l(P )) med(h(P ),w(P ),l(P )) . Finally we denote by V (P ) the volume of P .

Problem 1 (Minimizing-Surface-Cuboid-Partition (MSCuboidP)). Given a set of

n numbers {v 1 , . . . , v n } such that v k = xyz, and the cuboid Cu = [0, x] × [0, y] × [0, z], find for each v k a polyhedron P k of Cu such that V (P k ) = v k and P k = Cu minimizing Hs(P k ).
A general lower bound for this problem has been established by Ballard et al. [START_REF] Ballard | Minimizing Communication in Linear Algebra[END_REF] and comes from Loomis-Whitney inequality. It simply states that a polyhedron P of volume V (P ) minimizes the surface of its covering cuboid if and only if it is shaped as a cube. This implies the following lower bound:

Hs(P ) ≥ 3V (P ) 2 3 .
(1)

NP-Completeness

We prove in this section the NP-completeness of the decision problem associated to MSCuboidP, MSCuboidP-DEC.

Problem 2 (MSCuboidP-DEC). Given a set of n given numbers {v

1 , . . . , v n } such that v k = xyz, a cuboid Cu = [0, x] × [0, y] × [0, z] and a number K, is there a set of n polyhedra P k of Cu such that V (P k ) = v k , P k = Cu and
Hs(P k ) ≤ K ?

We start by reducing this problem to a more constrained variant named ACCuboidP, in which the goal is to partition the cuboid in cubes of specified side lengths.

Problem 3 (All-Cube-Cuboid-Partition (ACCuboidP)). Given a set of p given length {l 1 , . . . , l p } such that l 3 k = xyz, and a cuboid Cu = [0, x] × [0, y] × [0, z], is there a set of k cubes C k ∈ Cu such that V (C k ) = l 3 k and C k = Cu ? Lemma 1. ACCuboidP is NP-Complete.
It is easy to check that ACCuboidP belongs to NP. We prove NP-hardess of ACCuboidP with a method inspired from the hardness proof of the equivalent 2D problem [START_REF] Beaumont | Matrix Multiplication on Heterogeneous Platforms[END_REF], by using a reduction from 2-PART-EQUAL, a variant of the well-known 2-PART problem. The NP-completeness of 2-PART-EQUAL can be proven by a reduction from 2-PART. Indeed, adding a constant C to every element of the instance of 2-PART and then adding n (where n is the size of the instance) elements of size C to the instance itself creates an instance of 2-PART-EQUAL that has a solution if and only if the original instance of 2-PART has one. Our proof consists in two steps: from an instance of 2-PART-EQUAL, we first derive another set of numbers b i and prove that they can be partitioned into two equal sets if and only if the 2-PART-EQUAL instance has a solution. Then, we use the b i numbers to build an instance of ACCuboidP for which the existence of a packing is equivalent to partitioning the b i into two equal size sets.

i = a i + 3n × max i a i + D where D = 60M -(A mod 60M ) n = a i + M 2 + D.
In addition, let us set k = n 120 + A+nD 60M and S = 1 2 1≤i≤2n b i . One can prove that k is an integer (since n is a multiple of 120) and that S = 60k × M . In addition, let us notice that for all i, M 2 < b i (since D ≥ 0 and a i > 0) and b

i ≤ M . Indeed, D ≤ 60M n ≤ M 4 and a i + M 2 ≤ M ( 1 2 + 1 6n ) ≤ 4M 6 . Therefore b i ≤ 11M
12 ≤ M . Let us prove now that there exists a solution to our instance of 2-PART-EQUAL if and only if there exists a set

I ⊂ [1, n] such that i∈I b i = i / ∈I b i . If there exists I such that |I| = n and i∈I a i = i / ∈I a i , then i∈I b i - i / ∈I b i = i∈I a i - i / ∈I a i + (|I| -|I|)( M 2 + D) i / ∈I a i - i∈I a i = (|I| -|I|)( M 2 + D). Yet, i / ∈I a i -i∈I a i ≤ 2n × max a i = M 3 and M 2 ≤ M 2 + D. Therefore, (|I| -|I|) M 2 ≤ M 3 and (|I| -|I|) ≤ 2 3 < 1.
By symmetry, we obtain |I| = |I| and I is a solution to 2-PART-EQUAL.

Second Reduction In order to build the ACCuboidP instance that will be used in the reduction, we rely on a result from Walters [START_REF] Walters | Rectangles as Sums of Squares[END_REF] stating that it is possible to tile any cuboid with a number of cubes which is poly-logarithmic in the side lengths of the cuboid. We call the cubes in such a tiling Walters' cubes, and we denote by W S(X, Y, Z) a (poly-logarithmic size) set of cubes tiling the cuboid X × Y × Z.

Let us consider the following instance of ACCuboidP:

-A cuboid of size 11M × 15M × S (with S = 60k × M ). -20k cubes of length 6M . -24k cubes of length 5M . -30k cubes of length 4M . -20k cubes of length 3M . -∀i, a cube of length b i . -∀i, W S(M -b i , b i , b i ) and W S(M, M -b i , b i ).
with M , k and the b i 's defined from the a i 's as in the first reduction described above. One can see that the reduction is polynomial, since the sizes of the Walters' cubes sets are poly-logarithmic functions of the b i 's.

In the first part of the proof, we prove that if we can split the b i items in two equal sets, then the above set of cubes can be packed into the cuboid.

Let us first consider, for each i, the cube of length b i and the two associated Walters' cubes sets. Figure 1 Figure 3 shows how to tile a 11M × 15M rectangle with the corresponding squares, where both M × M squares represent a slice of the M × M × S cuboids presented above. This arrangement can be repeated for a total length of S, since

10k × 6M = 12k × 5M = 15k × 4M = 20k × 3M = 60k × M = S.
Hence, this provides a tiling of the whole 11M × 15M × S cuboid. For the second part of the proof, we need to prove if the cuboid can be tiled with the set of cubes, then a partition of the b i values in two equal size sets exists. We start by proving that in any valid tiling of the cuboid, the 11M ×15M rectangle can only be tiled as shown on Figure 3 (or under the same pattern but with an horizontal symmetry).

Let us note that, except the b i -cubes and the Walters' cubes, all cubes have length that are multiple of M . Therefore one can see that the resulting projections of the b i -cubes and the Walters' cubes on the 11M × 15M rectangle can be seen as several M × M squares.

× 20k + 120M 3 × k = 660M 3 × k.
On average over all the S slices of the cuboid, this represents a surface of 660M 3 ×k S = 11M 2 . We now prove that it is actually impossible to tile this 11M × 15M face with less than 11M 2 surface coming from these cubes.

Let s be the surface coming from these cubes divided by M 2 , s ∈ [0, 10]. Let p 6 , respectively p 5 and p 4 , the number or 6M -cubes, respectively 5M -cubes and 4M -cubes, used to tile the 11M × 15M -face. We have computed the possible values such that (15 × 11 -s)M 2 = (36p 6 + 25p 5 + 16p 4 )M 2 , they can be found in Table 1. We now consider each case one by one.

Case s p6 p5 p4

(1) 0 3 1 2 (2) 1 0 4 4 (3) 1 0 8 (4) 2 2 3 1 (5) 3 0 2 7 (6) 4 1 5 0 (7) 2 1 4
Case s p6 p5 p4 5 (8) 6 1 3 3 (9) 7 3 2 0 (10) 8 0 5 2 (11)

1 1 6 (12) 9 3 0 3 (13) 10 0 3 5 Table 1: Possible values of p 6 , p 5 and p 4 as a function of s.

Case (1) s = 0, p 6 = 3, p 5 = 1, p 4 = 2: 11 and 15 are odd, therefore there must be a least one fraction of an odd square in every line or column (we can see the rectangle as a grid 11 × 15). The odd cubes are the 5M -cubes, the 3M -cubes and the ones we can produces with b i -cubes and Walters' cubes. Since we have only one 5M -cube, the total length of odd squares is less than 15, so we cannot have one fraction of an odd square in every column. Therefore the tiling in this case is not possible.

Case (2) s = 1, p 6 = 0, p 5 = 4, p 4 = 4: As 5×4 > 15, we cannot have more than three 5-squares in a line and therefore, with four 5-squares, we have to be in one of the dispositions shown in Figure 4(a) and Figure 4(b). In both cases, there exist at least two columns where two 5-squares are superposed (denoted d i in the figures). However, the only way to complete these columns to a height of 11 is to use a square of size 1, and we have only one M -square since s = 1. Therefore the tiling in this case is not possible. 2), ( 4) and [START_REF] Becker | Towards Data Partitioning for Parallel Computing on Three Interconnected Clusters[END_REF].

Case (3) s = 1, p 6 = 2, p 5 = 0, p 4 = 8: By a parity argument similar to the one used in Case (1), we can prove the impossibility of the tiling.

Case (4) s = 2, p 6 = 2, p 5 = 3, p 4 = 1: It is easy to see that the 6-squares and the 5-squares can only be tiled in a way similar to Figure 4(c). Therefore there is no room to place the 4-square and the tiling is impossible in this case.

Case (5) s = 3, p 6 = 0, p 5 = 2, p 4 = 7: By a parity argument similar to the one used in Case (1), we can prove the impossibility of the tiling.

Case (6) s = 4, p 6 = 1, p 5 = 5, p 4 = 0: By reasoning on the number of 5-squares in a similar fashion than in Case (2), we can prove the impossibility of the tiling.

Case (7) s = 4, p 6 = 2, p 5 = 1, p 4 = 4: By a parity argument similar to the one used in Case (1), we can prove the impossibility of the tiling.

Case (8) s = 6, p 6 = 1, p 5 = 3, p 4 = 3: Since 5 + 5 + 6 > 15, we cannot have two 5-squares and one 6-square appearing on the same line. Then we have only two choices to tiles these squares. The first is shown on Figure 5(a). In this case there is a rectangle of size 9M × 6M that can be proved impossible to tile with six M -square and three 4M -squares. Therefore we have to be in the case of Figure 5(c), where the dashed zone does not intersect the last 5M -square (otherwise we are in the case in Figure 5(b) that is strictly harder to tile than the case of Figure 5(a)). Therefore the dashed zone, a square of size 5M × 5M has to be tiled with only fractions of 4M -squares and M -squares. In order to fill the gap of size 5M the only way to begin is the one shown in Figure 5(d).

To complete the last column, we have to fill the 7M -gap, and 7 = 5 + 2 × 1 is the only possibility (there are only two M -squares remaining). This yields the situation shown in Figure 5(e), and one can see that the tiling can not be finished with the remaining 4M -squares: the tiling is impossible in this case.

Figure 5: Tiling in case [START_REF] Beaumont | A New Approximation Algorithm for Matrix Partitioning in Presence of Strongly Heterogeneous Processors[END_REF].

Case (9) s = 7, p 6 = 3, p 5 = 2, p 4 = 0: It is impossible to have three 6Msquares. Indeed 2 × 6 > 11 and then we cannot have two 6M -squares appearing on the same column. In the same time, 3 × 6 > 15 and then we cannot have three 6M -squares appearing on the same line (see Figure 4(d) to have a better visualization). Hence there is no room to store three 6M -squares, and the tiling in this case is not possible.

Case (10) s = 8, p 6 = 0, p 5 = 5, p 4 = 2: By reasoning on the number of 5-squares in a similar fashion than in Case (2), we can prove the impossibility of the tiling.

Case (11) s = 8, p 6 = 1, p 5 = 1, p 4 = 2: By a parity argument similar to the one used in Case (1), we can prove the impossibility of the tiling.

Case (12) s = 9, p 6 = 3, p 5 = 0, p 4 = 3: By a parity argument similar to the one used in Case (1), we can prove the impossibility of the tiling.

Case (13) s = 10, p 6 = 0, p 5 = 3, p 4 = 5: We consider two sub-cases: either we have one 3M -rectangle and one M -rectangle, or we have ten M -rectangles. In the first case we notice that with the available rectangles there are only two ways to achieve a exact length of 11M : 1M + 5M + 5M and 3M + 4M + 4M . Both options create gaps that cannot be filled with the remaining rectangles, see Figure 6(a). In the second case (ten M square and no 3M -square), there are more options to achieve a length of 11M . The first one is to use two 4M -square and three M -square. This creates a 3M × 4M gap, and there are too few M -squares to fill it, see Figure 6(b), therefore we cannot tile this way. The second option is to use two 5M -squares and one M -square. In this case we have to use four more M -squares to complete and then we have again a length 11M to achieve, that, with the remaining resources, can only be done either with two 4M -squares and three M -squares, and we have proven above that it is not feasible, or with a 4M -square, a 5M -square and two M -squares, but in this case we cannot tile the resulting 2M ×3M rectangle because we have too few M -squares, see Figure 6(c). Another option is to use a 4M -square and seven M -squares, but this results in a gap of length 7M to fill and one can see that is not possible, see Figure 6(d) (replacing the 5M -square of the figure by a 4M -square yields the same problem). Another possibility is to use a 5M -square and six M -squares, for which the only reasonable continuation is the case shown on Figure 6(e) and this creates a gap of length 6M which cannot be tiled with the remaining squares. The last option is to use a 5M -square, a 4M -square and two M -squares but the tiling is still impossible, as it is shown on Figure 6(f). Therefore, in any sub-case, the tiling is impossible in this case. This proves that at every slice of the main cuboid, the surface coming from the 3M -cubes, the b i -cubes and the Walters' cubes is at least 11M 2 . Since this is exactly the average value, this implies that this surface is exactly 11M 2 at each slice. We can now observe that this requires at least two M -squares in any slice and once again this is exactly the average value of M -squares. Therefore, in any slice, there are exactly two M -squares and one 3M -square. The argument used to create Table 1 now shows that the slice necessarily includes two 6M -squares, two 5M -squares and two 4M -squares, and one can prove that the only way to tile the 11M × 15M rectangle with these squares is as shown on Figure 3.

We have built a pattern that must appear on each slice of the cuboid, and in which the b i cubes have to be included into two separate parts of the tiling. Let us denote by I the indexes of the b i cubes which appear in the leftmost M × M × S cuboid. Since by construction, b i > M 2 for all i, these cubes have to be arranged as depicted on Figure 2. This shows that

i∈I b i ≤ S and i / ∈I b i ≤ S.
Since the total sum is 2S, this implies that both sums are in fact equal to S, and thus that there exists a solution to the original 2-PART-EQUAL instance.

Note that the pattern in Figure 3 can be horizontally reversed. Furthermore, both possibles patterns can be present on the final tiling. However, even in this case, considering the b i -cubes on the left side still yields a set I such that 1) is tight if and only if there exist a partitioning where only cubes are used). Yet, thanks to Lemma 1, ACCuboidP is NP-complete. Therefore MSCuboidP-DEC is NPcomplete.

i∈I b i = i / ∈I b i = S.

Approximation Algorithm

In this section, we present 3D-NRRP, an approximation algorithm for the case where the cuboid to partition is cubic, what corresponds to the multiplication of square matrices. It is inspired by the NRRP algorithm proposed in [START_REF] Beaumont | A New Approximation Algorithm for Matrix Partitioning in Presence of Strongly Heterogeneous Processors[END_REF] and by Nagamochi et al. in [START_REF] Nagamochi | An Approximation Algorithm for Dissecting a Rectangle into Rectangles with Specified Areas[END_REF].

Presentation and correctness of 3D-NRRP

3D-NRRP (see Algorithm 1) is based on the divide and conquer principle, and its analysis relies on the following invariant: at each step, the aspect ratio of the cuboid to be partitioned is smaller than 3. In what follows, we define the aspect ratio ρ of a cuboid as the ratio of the largest length by the smallest length. We also define the second aspect ratio ρ as the ratio of the largest length to the median length.

At each step of the algorithm, the current cuboid (whose aspect ratio is smaller than 3) is split into two parts, and the same routine is recursively applied to each part. To ensure that the resulting parts have an aspect smaller than 3, the splitting is performed according to three modes. The first mode is the general case, in which the cuboid is partitioned in two disjoint cuboids by cutting along the largest length (Lines 7 to 18 in Algorithm 1, and Figure 7(a)). This is possible if there exists an index k such that k-1 i=1 v i ≥ v 3ρ2 . Indeed, Lemmas 2 and 3 show that this condition is sufficient to prove the invariant for both parts. More specifically, Lemma 2 states that in that case, both resulting cuboids have a total volume greater than one third of the overall volume, and Lemma 3 states that the aspect ratio of both cuboids is smaller than 3, under the assumption that the previous one had also a ratio less than 3, what ensures the correctness of the algorithm.

2 √ 3

 23

Problem 4 ( 2 -

 42 PART-EQUAL). Given a set of 2n integers {a 1 , . . . , a 2n }, is there I ⊆ [1, n] such that |I| = n and i∈I a i = i / ∈I a i First Reduction Let us now consider an instance of 2-PART-EQUAL, {a 1 , . . . , a 2n } and let us denote 2A = 1≤i≤2n a i and M = 6n × max i a i . Let us suppose, without loss of generality, that n is a multiple of 120 greater than 240 and let us define a new set {b 1 , . . . , b 2n } as ∀i, b

Figure 1 :Figure 2 W

 12 Figure 1: Tiling of a b i × M × M cuboid.

Figure 3 :

 3 Figure 3: Tiling of a 11M × 15M rectangle.

Figure 4 :

 4 Figure 4: Tiling in cases (2), (4) and (9).

Figure 6 :

 6 Figure 6: Tiling in case (13).

Theorem 1 . 2 / 3 k

 123 MSCuboidP-DEC is NP-complete.Proof. There is a reduction from ACCuboidP to MSCuboidP-DEC. Indeed,ACCuboidP ({l 1 , . . . , l p }, [0, x]×[0, y]×[0, z]) is true if and only if M SCuboidP -DEC({l 3 1 , . . . , l 3 p }, [0, x] × [0, y] × [0, z],3v ) is true (the bound in (

vi and α = v /v;

In the second mode, v n is significantly larger than the other v i 's. Splitting in two cuboids would result in the smallest cuboid having an aspect ratio larger than 3. Therefore, whenever possible (Lines 23 to 24 of Algorithm 1, and Figure 7(b)), 3D-NRRP shapes the smallest part as a cube included in the covering cuboid of the other part, which is made of one element only, namely v n . Lemma 4 states that we can put a such cube in the current cuboid under the condition of Line 23 of Algorithm 1. When the condition is not met, the edge length of this cube would be greater than the smallest dimension of the cuboid; thus in that case we create a cuboid like in Figure 7(c) (Line 26-32 of the Algorithm 1) by setting one dimension to the smallest dimension of the current cuboid and making the remaining two dimensions equal. Lemma 5 states that in this case, such a partitioning is possible (the edges of this small cuboid are smaller than the ones of the containing cuboid) and its aspect ratio is below 3.

Figure 7: The three splitting modes of 3D-NRRP. In all cases, the gray polyhedra is attributed to {v 1 , . . . , v k-1 }, the white one to {v k , . . . , v n }.

Lemma 2. Let {v 1 , . . . , v n } be a set of positive values sorted in non-decreasing order, ρ ≥ 1, and v = i v i . Let us assume that there exists an index k such that

, and let us consider the smallest such integer. Then

which is a contradiction with the fact that the v i 's are sorted in non-decreasing order.

Lemma 3. Let Cu be a cuboid of dimension w × h × l, with volume V = hwl, aspect ratio ρ, and second aspect ratio ρ . Let us assume that Cu 1 and Cu 2 are obtained by cutting Cu along the largest length, with V (Cu 1 ) and V (Cu 2 ) not smaller than V 3ρ . Then, ρ(Cu 1 ) ≤ max(3, ρ) and ρ(Cu 2 ) ≤ max(3, ρ). Proof. Let us suppose that w ≤ h ≤ l without loss of generality. In this case, w = w(Cu 1 ) = w(Cu 2 ), h = h(Cu 1 ) = h(Cu 2 ), ρ = l/w and ρ = l/h. Let us denote l 1 = l(Cu 1 ) and l 2 = l(Cu 2 ) and let us consider cuboid Cu 1 . Then, there are 3 cases to consider:

. Therefore

Thus, in all cases, ρ(Cu 1 ) ≤ max (3, ρ). By symmetry, the same proof applies to Cu 2 . Lemma 4. Let Cu be a cuboid of dimension w × h × l, with volume V = hwl, aspect ratio ρ, and second aspect ratio ρ . For any α > 0 such that αρ 2 ≤ ρ , then 3 √ αV ≤ min(w, h, l).

Proof. Without loss of generality, we can suppose w ≤ h ≤ l and therefore l = ρw and l = ρ h and finally h = ρ ρ w. Thus:

Lemma 5. Let Cu be a cuboid of dimension w × h × l, with volume V = hwl, aspect ratio ρ, and second aspect ratio ρ . Let α be such that ρ ρ 2 < α < 1 3ρ . Let us denote min d = min(w, h, l). Then α

Proof. Without loss of generality we can suppose min d = w ≤ h ≤ l and therefore l = ρw and l = ρ h and finally h = ρ ρ w. Thus:

w , and the previous result yields ρ(Cu ) ≤ h w = ρ ρ ≤ ρ.

Approximation Ratio

This section is devoted to the proof of Theorem 2, which states that the 3D-NRRP achieves a 5 6 2/3 approximation ratio ( 5

1.51).

Theorem 2. 3D-NRRP is a 5 6 2/3 -approximation when the given cuboid is cubic. The sketch of the proof is as follows. First, we prove that if {P 1 , . . . , P n } is an output of 3D-NRRP, then any output polyhedron P i satisfies Hs(Pi) 3V (Pi)

.

Remind that Equation (1) states that Hs(P * i ) ≥ 3V (P * i )

2 3 for any (optimal) solution {P * 1 , . . . , P * n }. By summing up these inequalities for all i, we get that i Hs(P i ) ≤ 5 The rest of the section is devoted to the proof that any polyhedron returned by 3D-NRRP satisfies the property stated above. One can see that there are only two situations in which 3D-NRRP returns a singleton: 2 and Line 33. In the first case, the returned zone is a cuboid with aspect ratio less than 3. In this case Lemma 7 provides the desired result, since 5

. We first start with a technical result in Lemma 6.

Proof. One can show that ∂f ∂x (x, y) = 2x-1-y 9y 1/3 x 4/3 . Since x ≥ 0 and y ≥ 0, ∂f ∂x (x, y) ≥ 0 if and only if x ≤ 1+y 2 . Hence, x → f (x, y) is decreasing on [1, 1+y 2 ] and increasing in [ 1+y 2 , y], which implies f (x, y) ≤ max(f (1, y), f (y, y)). f (1, y) = y+2 √ y : the previous statements show that f (x, y) ≤ g(y). One can prove that g (y) = 2(y-1) 9y 4/3 . Therefore g is increasing on [1, ∞] and, in the considered case, g(y) ≤ g(3) = 5 3 3 √ 3 , which proves the desired result.

Lemma 7. If P is a cuboid with ρ(P ) ≤ 3, then Hs(P ) 3V (P ) 2/3 ≤ 5 3 3 √ 3 . Proof. We denote ρ(P ) = ρ and V (P ) = V . We suppose that w = w(P ) ≤ h = h(P ) ≤ l(P ) = l without loss of generality. We denote x = h/w. In this case l = ρw and V = whl = ρxw 3 . Therefore

where f is as defined in Lemma 6, what ends the proof.

In the other case, 3D-NRRP returns a cuboid minus a cube (Line 22 of Algorithm 1, as described on Figure 7(b)). The bound on the volume of the cube is such that the conditions of Lemma 9 are fulfilled. As before, this relies on the technical Lemma 8.

Proof. One can prove that ∂f ∂x (x, y) = 1+y 3 x 2 +( 7 3 +y)yx-2y 2 9(yx-x 2 3 ) 5/3

. Then ∂f ∂x has the same sign than P y (x) = 1+y 3 x 2 +( 7 3 +y)yx-2y 2 . Since P y is a second order polynomial with P y (0) = -2y 2 < 0 and lim x→+∞ p y (x) = +∞, P y is either positive on [1, y],

negative on [1, y], or negative then positive. Therefore x → f (x, y) on [1, y] is either increasing, decreasing or decreasing and then increasing. In any case, f (x, y) ≤ max(f (1, y), f (y, y)).

Let g(y) = f (1, y) = 1+2y 3(y-1/3) 2/3 . One can show that g (y) = 2(y-2)

9(y-1/3) 5/3 , hence g is decreasing on [START_REF] Anderson | Communication-avoiding QR Decomposition for GPUs[END_REF][START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF] and increasing on [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF][START_REF] Ballard | Communicationoptimal Parallel Algorithm for Strassen's Matrix Multiplication[END_REF]. Therefore g(y) ≤ max(g(1), g(3)) = max(( 32 ) 2/3 , 7

12y 2 , hence h is increasing on [START_REF] Anderson | Communication-avoiding QR Decomposition for GPUs[END_REF][START_REF] Ballard | Communicationoptimal Parallel Algorithm for Strassen's Matrix Multiplication[END_REF]. Therefore h(y) ≤ h(3) = 5 6 2/3 . Putting it all together, we get f (x, y) ≤ max(( 32 ) 2/3 , 5 6 2/3 ) = 5 6 2/3 .

Lemma 9. If V (P ) ≥ (1 -1 3ρ (P ) )V (Cu(P )) and ρ(P ) ≤ 3, then Hs(P ) 3V (P ) 2/3 ≤ 5 6 2/3 . Proof. Let us denote ρ = ρ(P ), ρ = ρ (P ), V = V (P ), and let us denote by w, h, l for the dimensions of Cu(P ). Let us suppose that w ≤ h ≤ l without loss of generality. Then, l = ρw and l = ρ h, and let us also denote x = h/w = ρ/ρ . With such notations, we get V (Cu(P )) = ρxw 3 and Hs(P ) = (ρ + x(1 + ρ))w 2 . Thus, the condition on V (P ) can be written as

where f is as defined in Lemma 8, what provides the conclusion.

Conclusion

We introduce a model of the partitioning problem associated to the 2.5D matrix multiplication algorithm on heterogeneous resources, a problem of crucial importance in high performance computing. This corresponds to partitioning a cuboid into several polyhedra, each representing the volume of computations attributed to a resource. We provide two theoretical results: a proof of the NP-completeness for this problem, and an approximation result for 3D-NRRP, which generalizes the results obtained in the 2D case. This is the first known approximation result for this problem, and it provides a strong guarantee ( 5

1.51). In addition, the computational time of the algorithm is extremely low, O(n log n), where n is the number of processors, what makes it perfectly suitable for practical use. This work opens several interesting perspectives for dimensions higher than 3, what corresponds to tensor products. It would also be interesting to combine this algorithm, whose general goal is to minimize the overall volume of communications while enforcing perfect load balancing between heterogeneous resources, with algorithms that explicitly take into account memory constraints at each node.