C Soize
email: c.soize@u-pem.fr

R Ghanem
email: ghanem@usc.edu

C Safta
email: csafta@sandia.gov

X Huan
email: xhuan@sandia.gov

Z P Vane
email: zvane@sandia.gov

J Oefelein
email: joefelein@sandia.gov

G Lacaze
email: glacaze@sandia.gov

H N Najm
email: hnnajm@sandia.gov

Q Tang
email: tang30@llnl.gov

X Chen
email: chen73@llnl.gov

Entropy-based closure for probabilistic learning on manifolds

Keywords: Statistical learning, Probabilistic learning, Concentration of probability, Measure concentration, Probability distribution on manifolds, Random sampling generator, MCMC generator, Diffusion maps, Smoothing parameter, Entropy principle

In a recent paper, the authors proposed a general methodology for probabilistic learning on manifolds. The method was used to generate numerical samples that are statistically consistent with an existing dataset construed as a realization from a non-Gaussian random vector. The manifold structure is learned using diffusion manifolds and the statistical sample generation is accomplished using a projected Itô stochastic differential equation. This probabilistic learning approach has been extended to polynomial chaos representation of databases on manifolds and to probabilistic nonconvex constrained optimization with a fixed budget of function evaluations. The methodology introduces an isotropic-diffusion kernel with hyperparameter ε. Currently, ε is more or less arbitrarily chosen. In this paper, we propose a selection criterion for identifying an optimal value of ε, based on a maximum entropy argument. The result is a comprehensive, closed, probabilistic model for characterizing data sets with hidden constraints. This entropy argument ensures that out of all possible models, this is the one that is the most uncertain beyond any specified constraints, which is selected. Applications are presented for several databases.

A lower case letter such as x or η is a real deterministic variable. A boldface lower case letter such as x or η is a real deterministic vector. An upper case letter such as X or H is a real random variable. A boldface upper case letter such as X or H is a real random vector. A lower case letter between brackets such as [x] or [η] is a real deterministic matrix. A boldface upper case letter between brackets such as [X] or [H] is a real random matrix.

Introduction

Objective and novelty of the paper. This work is a continuation and final installment for previous work [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF] devoted to a methodology for probabilistic learning on manifolds. Starting with a specific dataset, the algorithm constructs a statistical model of the data, together with a numerical generator to sample additional statistically consistent realizations. Specifically, the procedure involves three steps consisting of (1) identifying the non-Gaussian probability distribution of a random vector for which the solely available information consists of a database made up of independent samples (generated by experiments and/or by numerical simulations), [START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF] delineating an intrinsic manifold in the data and, (3) constructing a Markov Chain Monte Carlo (MCMC) generator of realizations on the manifold with the non-Gaussian distribution as its invariant measure. This approach introduces two hyperparameters. The first one, ε, is related to the parameterization of the isotropic-diffusion kernel related to the transition matrix of the diffusion maps used to discover the manifold. The second one, L, is a cut-off threshold for determining the dimension m of the diffusion maps basis used for constructing the MCMC generator. This probabilistic learning algorithm (that depends on ε and L) facilitates the detection and the construction of the intrinsic underlying probability distribution from which the database has been generated. The parameter L is associated with scale separation and a reasonable universal value for it has been deduced from considerations of a large number of datasets from a variety of disciplines. The novelty of the present paper is to propose a selection criterion for ε, resulting in a data-driven model that is algebraically closed.

Role played by statistical learning in computational physics and engineering sciences. Statistical and probabilistic learning methods have been extensively developed (see for instance, [START_REF] Schölkopf | Kernel principal component analysis[END_REF][START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF][START_REF] Aggarwal | Mining Text Data[END_REF][START_REF] Dalalyan | Sparse regression learning by aggregation and langevin monte-carlo[END_REF][START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF][START_REF] Balcan | Statistical active learning algorithms[END_REF][START_REF] James | An Introduction to Statistical Learning[END_REF][START_REF] Dong | Knowledge vault: A web-scale approach to probabilistic knowledge fusion[END_REF][START_REF] Ghahramani | Probabilistic machine learning and artificial intelligence[END_REF][START_REF] Taylor | Statistical learning and selective inference[END_REF]) and play an increasingly important role in computational physics and engineering science [START_REF] Ghanem | Handbook of Uncertainty Quantification[END_REF], in particular for design optimization with underlying stochastic operators and constraints using large scale computational model, and more generally in artificial intelligence for extracting information from big data. In this framework, statistical learning methods have been developed in the form of surrogate models from which approximations of the expensive functions can easily be evaluated [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Queipo | Surrogate-based analysis and optimization[END_REF][START_REF] Byrd | On the use of stochastic Hessian information in optimization methods for machine learning[END_REF][START_REF] Homem-De Mello | Monte Carlo sampling-based methods for stochastic optimization[END_REF]. Although Gaussian process models are most commonly used in this context [START_REF] Kleijnen | Constrained optimization in expensive simulation: Novel approach[END_REF], alternative approaches based on Bayesian methods have been proposed [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Wang | Bayesian optimization in a billion dimensions via random embeddings[END_REF][START_REF] Xie | Bayesian optimization via simulation with pairwise sampling and correlated pair beliefs[END_REF]. For the evaluations of expensive stochastic functions in the presence of uncertainties, current computational challenges remain significant enough to require some degree of probabilistic approximation [START_REF] Du | Sequential optimization and reliability assessment method for efficient probabilistic design[END_REF][START_REF] Byrd | On the use of stochastic Hessian information in optimization methods for machine learning[END_REF][START_REF] Eldred | Design under uncertainty employing stochastic expansion methods[END_REF][START_REF] Yao | Review of uncertainty-based multidisciplinary design optimization methods for aerospace vehicles[END_REF].

A classical example of the use of the probabilistic learning on manifolds is the following. We consider an expensive large-scale stochastic computational model that discretizes a complex system modeled by a boundary value problem. The vector-valued random response is written as Q = f(W, U) in which W is a non-Gaussian vector-valued random parameter controlling the system, U is a non-Gaussian vector-valued random parameter representing uncertainties, and f is a deterministic nonlinear mapping representing the computational model. For instance U corresponds to the spatial discretization of a non-Gaussian tensor-valued random field that is a coefficient of a partial differential operator of the boundary value problem on which the computational model is based. The only available information is a given initial dataset (training set) of length N, which is constructed as the set of N points {x j d = (q j d , w j d ), j = 1, . . . , N} in which q j d = f(w j d , U(θ j )) are N independent realizations of Q (calculated with the expensive computational model) where w j d and U(θ j ) are N independent realizations of W and U. Consequently, {x j d , j = 1, . . . , N} are N independent realizations of the non-Gaussian random vector X = (Q, W). Knowing only this initial dataset, the objective is, for instance, to construct, for any given w in its admissible set, an estimate h (N) (w) of h(w) that is defined, for instance, by h(w) = E{H(Q)|W = w} in which E is the conditional mathematical expectation given W = w and where H is a given deterministic mapping (for instance, h(w) could be the objective function of an optimization problem for which w would be the design parameter). If each evaluation q j d is expensive in CPU time, then N will be, generally, not sufficiently large for obtaining a good convergence of h (N) (w) towards h(w). The probabilistic learning on manifolds allows for generating ν sim N additional independent realizations {x 1 ar , . . . , x ν sim ar } of random vector X, without using the expensive computational model, which allows for deducing ν sim additional realizations (q ar , w ar ), = 1, . . . , ν sim } such that (q ar , w ar ) = x ar . We can then construct a better estimate of the conditional expectation that is required for computing h (N) (w). We are then considering a probabilistic machine learning for the small-data challenge (N small) in computational science.

Brief description of probabilistic learning under consideration and definition of the problem to be solved. As explained above, the authors have proposed a general probabilistic learning on manifold [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF] for generating additional realizations of an R n -valued random variable X = (X 1 , ..., X n ) for which the available information is only made up of a given set of N n points {x 1 , . . . , x N } that are N independent realizations in R n of random vector X. This method is based (1) on the use of nonparametric statistics for estimating the probability density function x → p X (x) with respect to the Lebesgue measure dx on R n of random vector X and then for deducing the probability density function

[x] → p [X] ([x]) with respect to the volume element d[x] on M n,N of the random matrix [X] = [X 1 , . . . , X N ]
with values in M n,N in which X 1 , . . . , X N are N independent copies of X, (2) on the construction of nonlinear Itô stochastic differential equation (ISDE) on M n,N , formulated for a dissipative Hamiltonian dynamical system for which p

[X] ([x]) d[x]
is the invariant measure, (3) on the use of diffusion maps for discovering the geometrical structure of the given set of points, and (4) on the construction of a reduced ISDE on M n,m obtained by projecting the original ISDE on a subspace of dimension m N spanned by a subset of the diffusion-maps basis and allowing additional realizations of X to be generated. This approach introduces two hyperparameters. The first one is the isotropic-diffusion kernel hyperparameter ε allowing the transition matrix of the diffusion-maps approach to be constructed from the given set of points. The second one is dimension m of the projection. Concerning hyperparameter ε related to the isotropic diffusion kernel, the authors did not find any proposed robust method in the open literature for estimating it on the basis of a quantitative criterion independent of the application under consideration. In addition, in the framework of the probabilistic learning on manifolds proposed in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF], dimension m is deduced from ε once a scale separation threshold, L, has been set. In this paper, we propose a selection model of ε based on the use of an entropy principle using only the given set of points. The algebraic basis constructed through the diffusion manifold approach is critical to the projected Itô equation developed in the proposed probabilistic learning method. Other manifold detection procedures, such as kernel PCA [START_REF] Schölkopf | Kernel principal component analysis[END_REF], that are endowed with such algebraic constructs could be used as substitutes for the manifold detection and characterization step. The probabilistic learning method proposed in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF] has been extended to polynomial chaos representation of databases on manifolds [START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF] and to probabilistic nonconvex constrained optimization with a fixed number of function evaluations [START_REF] Ghanem | Probabilistic nonconvex constrained optimization with fixed number of function evaluations[END_REF].

Recently, the robustness of this approach has been tested with success on different problems such as the optimal well-placement [START_REF] Ghanem | Optimal well-placement using a probabilistic learning, Data-Enabled Discovery and Applications[END_REF][START_REF] Thimmisetty | Homogeneous chaos basis adaptation for design optimization under uncertainty: Application to the oil well placement problem[END_REF], the prediction of maximum daily precipitation with data generated from large scale climate models, the Continental United States (CONUS) Regionally-Refined Model (RRM) [START_REF] Roesler | Climatology of the accelerated climate model for energy's community atmospheric model, version 5, configured with variable resolution over the continental united states[END_REF] of Energy Exascale Earth System Model (E3SM) V0 [START_REF] Terai | The atmospheric hydrologic cycle in the acme v0.3 model[END_REF], and the enhancing of the predictability of a two-dimensional (2D) computational model for a scramjet [START_REF] Soize | Enhancing model predictability for a scramjet using probabilistic learning on manifold[END_REF].

Organization of the paper. The objectives and the organization of the paper are as follows. Despite the fact that the reader can find the details of the considered probabilistic learning algorithm on manifolds in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF], and taking into account that the objective of the paper is to propose a selection model for its hyperparameters, ε and m, it seems reasonable to start by summarizing this algorithm. Consequently, Section 2 presents such a summary of the probabilistic learning on manifolds that underscores the role played by hyperparameters ε and m. Section 3 is devoted to a selection model for calculating an optimal value ε opt of hyperparameter ε using an entropy principle and a method for computing an adapted dimension m L (ε) of hyperparameter m(ε) as a function of ε and consequently, deducing the optimal value m opt of m(ε) associated with ε opt . Section 4 deals with a reanalysis of Application 2 presented in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF] (random 3D-data around a helix) using the proposed selection model. In Section 5, a validation of the proposed selection model of the hyperparameters is presented for several databases: Scramjet-d8 and ScramJet-d16 that correspond to simulations of a ScramJet performed by Sandia National Laboratories in Livermore with two large scale complex computational models [START_REF] Soize | Enhancing model predictability for a scramjet using probabilistic learning on manifold[END_REF], and Climate-LS-21-34 and Climate-LS-21-50 that correspond to the maxima of daily precipitations, which have been predicted by Lawrence Livermore National Laboratory with a climate computational model [START_REF] Bader | Accelerated climate modeling for energy (acme) project strategy and initial implementation plan[END_REF]. Section 6 presents an analysis of the results obtained.

Summary of the probabilistic learning on manifolds

The probabilistic learning on manifold proposed in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF] uses only a given set of N points {x 1 d , . . . , x N d } in R n , which are assumed to be N independent realizations of a random vector X with values in R n . The probability distribution of X is unknown and is assumed to be concentrated in a neighborhood of a subset of R n (a manifold) that is also unknown and that has to be discovered. We introduce the matrix [

x d ] = [x 1 d . . . x N d ] ∈ M n,N
that is one realization of the random matrix [X] = [X 1 , . . . , X N ] with values in M n,N in which X 1 , . . . , X N are N independent copies of X. For the dataset that will be used in this paper, vector X will consist of the parameters of a computational model as well as the quantities of interest that are obtained as output of the computational model. The objective of the probabilistic learning on manifold is to construct a probabilistic model using only the given dataset represented by [x d ]. This model is then used to generate additional independent realizations {x 1 ar , . . . , x ν sim ar } in R n of random vector X. The proposed method preserves the concentration of the additional realizations around the manifold. For the applications described subsequently in the paper, a very large number, ν sim N, of additional realizations can then be generated for the quantities of interest, which then allow for estimating the probability density functions of various quantities of interest (QoI), including extreme values statistics. The main steps of the probabilistic learning algorithm on manifold are summarized in Section 2.1. Section 2.2 deals with the important issue concerning the convergence analysis with respect to hyperparameter m for a given value of ε and the non-adaptation of such a convergence analysis to the objective of the probabilistic learning method. Section 2.3 deals with the mathematical formulation of an efficient method for calculating dimension m as a function of ε.

Probabilistic learning algorithm on manifold

The 8 steps of the probabilistic learning algorithm used in this paper are summarized next.

1) The initial given dataset is made up of unscaled data represented by matrix [

x uns d ] in M n,N . The matrix [x d ] in M n,N of the scaled given dataset (simply called the given dataset) is constructed such that [x d ] k j = ([x uns d ] k j -min j [x uns d ] k j )/(max j [x uns d ] k j -min j [x uns d ] k j
) for all k = 1, . . . , n and j = 1, . . . , N. This step is generally required to avoid numerical problems in the second step introduced below. 2) Let [c] ∈ M n be the empirical estimate of the covariance matrix of X. We consider the eigenvalue problem

[c] ϕ k = µ k ϕ k . Let ν ≤ n be the number of positive eigenvalues {µ k } ν k=1 with 0 < µ 1 ≤ µ 2 ≤ . . . ≤ µ ν and let [ϕ] be the (n × ν) matrix such [ϕ] T [ϕ] = [I ν ], whose columns are the associated orthonormal eigenvectors ϕ 1 , . . . , ϕ ν . A principal component analysis [X] = [x] + [ϕ] [µ] 1/2 [H] of random matrix [X]
is carried out in order to normalize/reduce the dataset in which [x] ∈ M n,N is the empirical estimate of the mean value of [X] and where [µ] is the positive diagonal (ν × ν) real matrix such that [µ] kk = δ kk µ k . The columns H 1 , . . . , H N of random matrix [H] with values in M ν,N are independent random vectors with values in R ν . This normalized/reduced representation allows for computing a new normalized dataset of N points {η

1 d , . . . , η N d } in R ν , represented by the matrix [η d ] = [η 1 d . . . η N d ] in M ν,N that is computed by [η d ] = [µ] -1/2 [ϕ] T ([x d ] -[x]
). If ν is chosen as the rank of matrix [c], then the constructed representation of [X] induces no error and corresponds to a pure normalization. If ν is chosen less that the rank of [c] in order to construct a reduced-order representation of [X], then ν can be selected for obtaining a given tolerance of the relative mean-square error defined by error pca [START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF] of the classical multidimensional Gaussian kernel-density estimation method [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis[END_REF][START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF] is then used to construct an estimate of the probability density function

(ν) = [x d ] -[x (ν) d ] F / [x d ] F where [x (ν) d ] = [x] + [ϕ] [µ] 1/2 [η d ]. 3) A modification
[η] → p [H] ([η]) with respect to the volume element d[η] on M ν,N of random matrix [H]. 4)
An MCMC generator for random matrix [H] is constructed using the approach proposed in [START_REF] Soize | Construction of probability distributions in high dimension using the maximum entropy principle. applications to stochastic processes, random fields and random matrices[END_REF][START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF] belonging to the class of Hamiltonian Monte Carlo methods [START_REF] Soize | Construction of probability distributions in high dimension using the maximum entropy principle. applications to stochastic processes, random fields and random matrices[END_REF][START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF][START_REF] Neal | MCMC using hamiltonian dynamics[END_REF], which is an MCMC algorithm [START_REF] Spall | Introduction to Stochastic Search and Optimization[END_REF]. The realizations of random matrix [H] could be obtained by solving a (ν × N) matrix-valued ISDE that corresponds to a stochastic nonlinear dissipative Hamiltonian dynamical system, for which

p [H] ([η]) d[η]
is the unique invariant measure (in fact, we will not need to generate realizations of [H] using this ISDE). 5) The diffusion-maps approach [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF][START_REF] Coifman | Diffusion maps, applied and computational harmonic analysis[END_REF] is then used to discover and characterize the geometrical structure of the normalized dataset [η d ] thus defining a manifold embedded in R N .

The transition matrix of a Markov chain relative to the geometrical structure of the given normalized dataset [η d ] is constructed thanks to the introduction of the isotropic diffusion kernel,

k ε (η, η ) = exp(-1 4ε η -η 2 ), defined on R ν × R ν in which ε > 0 is a hyperpa- rameter of the model. Let [K] be the symmetric matrix in M N with positive entries such that [K] i j = k ε (η i d , η j d ) for i and j in {1, . . . , N}. Let [b] be the positive-definite diagonal real matrix in M N such that [b] i j = δ i j N j =1 [K] i j and let [P] be the matrix in M N such that [P] = [b] -1 [K]. Matrix [P]
that has positive entries satisfying N j=1 [P] i j = 1 for all i = 1, . . . , N, can be viewed as the transition matrix of a Markov chain that yields the probability of transition in one step. We consider the first m largest eigenvalues of the gen-

eralized eigenvalue problem [K] [ψ] = [b] [ψ] [Λ] such that [ψ] T [b] [ψ] = [I m ], where [I m ] is the identity matrix in M m and where [ψ] T [K] [ψ] = [Λ]. The diagonal matrix [Λ] in M m is made up of the eigenvalues λ 1 , . . . , λ m that are such that 1 = λ 1 > λ 2 ≥ . . . ≥ λ m . The corresponding eigenvectors ψ 1 , . . . , ψ m in R N are such that [ψ] = [ψ 1 . . . , ψ m ] ∈ M N . The diffusion-maps basis of dimension m is then defined by the m vectors g 1 , . . . , g m in R N such that g α = λ α ψ α . We introduce the (N × m) matrix [g] = [g 2 . . . g m ]
with m = m -1, which will be used below for performing the projection of random matrix [H]. In general, m can be selected such that m N (see Section 2.2). Note that g 1 is always a vector whose components are all equal. Since the random matrix [H] is centered, this vector can be removed from matrix [g] used for performing the projection of [H]. Finally, since eigenvalues λ 2 , . . . , λ m and matrix [g] depend on ε, we will rewrite them as λ 2 (ε), . . . , λ m (ε) and

[g ε,m ]. 6) As proposed in [1], a reduced-order representation [H ε,m ] of [H] is introduced such that [H ε,m ] = [Z ε,m ] [g ε,m
] T is constructed on the manifold in which [Z ε,m ] is a random matrix with values in M ν,m for which m N and a M ν,m -valued reduced-ISDE is obtained by projecting the M ν,N -valued ISDE onto the diffusion manifold by using the reduced-order basis represented by matrix [g ε,m ] T . It should be noted that such a projection corresponds to a reduction of the dataset dimension and not to a reduction of the physical components of random vector H (this latter reduction already results from a statistical reduction introduced in step 2). Such a projection preserves the concentration of the generated realizations around the manifold.

], . . . , [x n MC ar ] of [X ε,m ] such that [x ar ] = [x] + [ϕ] [µ] 1/2 [η ar ].
The algorithm for generating these additional realizations is given in Appendix A. 8) Reshaping these n MC matrices yields the ν sim = N × n MC N additional realizations {x 1 ar , . . . , x ν sim ar } in R n of random vector X ε,m .

2.2. Remark concerning the convergence analysis with respect to hyperparameter m for a given value of ε and its non-adaptation to the objective of the probabilistic learning method For a given application for which the value of hyperparameter ε is fixed, the corresponding probability density function It should be noted that for m = N, the family of vectors {g 1 , g 2 , . . . , g N } constitutes a vector basis of R N . An estimation of m could be constructed with a convergence analysis aimed at reducing the relative mean-square error defined by error pro j ( m)

[x] → p [X ε,m ] ([x]; ε, m), with respect to d[x], of random matrix [X ε,m ] corresponds
= [x d ] -[x (ν, m) d ] F / [x d ] F in which [x (ν, m) d ] = [x] + [ϕ] [µ] 1/2 [z d ] [g] T with [z d ] = [y d ] [g] ([g] T [g]) -1 .
As this error goes to zero, m goes to N and the probability measure

p [X ε,m ] ([x]; ε, N) d[x] approaches the probability measure p [X] ([x]) d[x].
In such a limit, the concentration of data around the manifold would not be leveraged, and there would be no gain achieved using the proposed probabilistic learning algorithm [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF]. Thus, the L 2 convergence analysis described above is not consistent with the objectives of the proposed probabilistic learning methodology. As already explained, the main objective of this probabilistic learning method is to take into account the geometrical structure of the given dataset in order to enrich the construction of the probability measure p

[X] ([x]) d[x] substituting it by p [X ε,m ] ([x]; ε, m) d[x]
that preserves the concentration of the additional realizations computed with steps 7 and 8 of the algorithm presented in Section 2.1. Thus, rather than pursuing a convergence analysis by increasing m, an optimal value of m is selected for each value of ε. In several of the numerical examples presented in Section 5, probability density estimates have also been constructed for extreme values statistics of random QoIs. The foregoing discussion has the implication that the convergence of these extreme values statistics with respect to m should not be of concern, as m is selected with the criterion of sufficient scale separation. The model is completely specified by selecting numerical values for m and ε. Consequently, these statistics correspond to the extreme values statistics of the QoI that are defined relative to this selected model.

Method for calculating an adapted dimension m L (ε) as a function of ε

As explained in Section 1, for a given application, the value of hyperparameter ε is currently selected arbitrarily, following a trial and error procedure with no quantitative criteria.

A plausible framework for separating scales in a given dataset would first involve evaluating, for different values of ε, a graph of the eigenvalues α → λ α (ε) for α ≥ 0, similar to the ones displayed in Fig. 1. It is clear that different values of ε result in different decays in the eigenvalues, and a different value of an optimal m corresponding to the rank of the eigenvalue at which significant drop in the spectrum is observed. It is also clear from these graphs that the lagest eigenvalue is always equal to 1 and its eigenvector is constant, an artifact of the normalization of the Laplacian of the data (ie the matrix obtained from discretizing the diffusion kernel). In all subsequent calculations, this largest "unit" eigenvalue is therefore ignored when comparing spectral content of different models. A method for identifying ε objectively is presented in Section 3. In the rest of this section, we present a procedure for evaluating the optimal value of m for a fixed value of ε. The model proposed for calculating the adapted value m L (ε) of m(ε) as a function of ε can be written as,

m L (ε) = -1 + arg min α | α≥3 λ α (ε) λ 2 (ε) < L , (1)
and then the corresponding adapted dimension m L (ε) of dimension m is such that

m L (ε) = m L (ε) -1 . ( 2
)
Parameter L is independent of ε and is chosen for separating the existing scales in any given dataset. Following our proposed procedure for estimating ε, L will be the only remaining free parameter in our data-driven model, although its significance as a scale separation threshold is obvious form context. All the numerical applications performed by the authors for dataset coming from different fields show that L = 0.1 is, generally, a good value and this choice corresponds to one order of magnitude for a scale separation. The analyses that will be performed in the future for other datasets will allow this choice of the value of L to be confirmed or improved. For the case illustrated in Fig. 1 The criterion defined by Eq. ( 1), is based on the idea that, for a fixed value of ε, the smallest dimension m( ) has to be chosen for obtaining a scale separation and for preserving the concentration of the additional realizations generated by the probabilistic learning method that has been proposed [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF]. In this paper, we propose a selection model for calculating an optimal value ε opt of ε and then, deducing the corresponding optimal value m opt = m L (ε opt ) of dimension m. This model is based on the use of an entropy principle and is presented in Section 3.

Principles for selecting the optimal value of hyperparameter ε

The principles that are proposed for selecting the optimal value ε opt of hyperparameter ε and the corresponding optimal value m opt of hyperparameter m are the following:

1) An admissible finite set A = {ε 1 < ε 2 < . . . < ε n ε } of n ε ordered values of ε is defined.
The optimal value ε opt of ε will be searched for in A. Consequently, for a given database, set A must be chosen to contain a sufficiently large number of possible values of ε, but the selection model that is proposed does not directly depend on this choice. 2) For each value of ε in A, the adapted value m L (ε) of m(ε) is computed with Eq. ( 1) for the reasons given in Section 2.3. The corresponding adapted dimension m L (ε) of hyperparameter m is then computed using Eq. ( 2).

3) The mapping ε → m L (ε) from A into N being known, we define the subset E 0 of A, made up of all the values of ε that minimize m L (ε) on A,

E 0 = {ε ∈ A | ε = arg min ε ∈A m L (ε )} , (3)
In general, the number of elements in set E 0 is greater than 1 (but is less or equal to n ε ). This is one of the reasons why ε opt cannot simply be constructed as the value of ε that minimizes m L (ε). 4) We now construct the subset E ⊂ A of the values of ε such that

E = ε ∈ A | m min L ≤ m L (ε) ≤ m max L , (4)
in which the lower bound m min L is the minimum value of m L (ε) for ε ∈ A,

m min L = min ε∈A m L (ε) , (5)
and where m max L is the upper bound defined by

m max L = (1 + χ) m min L , (6)
with r denotes the floor integer part of r (for instance, 2.3 = 2), and where χ is a small real number such that 0 ≤ χ < 1. It can be seen that

E 0 E A . ( 7
)
The optimal value ε opt of ε will be searched into E and not in E 0 in order to construct a subset of A, the subset E, which is more robust than E 0 with respect to the small fluctuations effects that can occur with Eq. ( 1). Based on our experience of the analysis of data sets from a range of applications, for the value L = 0.1 recommended in Section 2.3, the value χ = 0.2 is a good associated value. It should be noted that χ is not a hyperparameter of the probabilistic learning approach as is parameter ε whose value depends on the database and for which a procedure is proposed for calculating it, but χ is a parameter whose value is independent of the database and for which a fixed value is proposed resulting from the analyses of several databases (may be this value of χ could be improved with future experimentations). 5) Subset E being known, the optimal value ε opt is sought as the value of ε in E, which maximizes the uncertainty of the family of random matrices {[H ε,m L (ε) ], ε ∈ E}. This is equivalent [START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Gray | Entropy and Information Theory[END_REF][START_REF] Ghanem | Handbook of Uncertainty Quantification[END_REF][START_REF] Soize | Uncertainty Quantification[END_REF] to maximizing the Shannon entropy [START_REF] Shannon | A mathematical theory of communication[END_REF] with respect to ε for the probability density function of random matrix [H ε,m L (ε) ]. This entropy S (ε) is then written as

S (ε) = - M ν,N p([η] ; ε) log(p([η] ; ε)) d[η] , (8)
in which we have introduced the following simplified notation,

p([η]; ε) = p [H ε , m L (ε) ] ([η] ; ε, m L (ε)) . (9)
Consequently, the use of the maximum entropy principle allows the optimal value ε opt of ε to be selected in E, solving the optimization problem,

ε opt = max ε ∈ E S (ε) . (10)
6) The corresponding optimal value m opt of hyperparameter m is then written as,

m opt = m L (ε opt ) . (11)
Remarks.

1) An alternative selection model would consist of replacing principles 4 and 5 above by one consisting in maximizing the entropy over A, that is to say ε opt = max ε ∈ A S (ε). However, such a principle would favor the selection of epsilon with respect to the criterion of minimizing m opt subject to scale separation.

2) The number of random variables in random matrix

[Z ε,m ] is ν × m. Consequently, for ε fixed, random matrix [H ε,m ] = [Z ε,m ] [g ε,m
] T , whose dimension (ν × N) is independent of m, depends on ν × m random variables. Consequently, the entropy S (ε) should "grow on average" with m = m L (ε). This phenomenon has been observed for all the databases analyzed (see Figs. 5, 7, 9, and 11) except for the simple database associated with a 3D helix (see Fig. 3) for which that is not true for the two small values 0.5 and 1 of ε, but what is true for ε ≥ 2. This growth of entropy with m is consistent with the observation that as m increases, less structure is delineated in the data (less localization), resulting in greater uncertainty about mechanisms underlying the fluctuations.

Computational aspects

The computational aspects concern the numerical calculation of S (ε) for ε given in A, the solution of the optimization problem, and the convergence analysis with respect to the number n MC of the number of realizations.

Numerical calculation of the entropy

The entropy defined by Eq. ( 4) is classically rewritten as,

S (ε) = -E{log(p([H ε,m L (ε) ] ; ε))} , (12)
For ε fixed in A and for m = m L (ε), the additional realizations of random matrix [H ε,m ], which are computed as explained in step 7 of Section 2.1, are denoted by [η 1 ar (ε)], . . . , [η n MC ar (ε)]. For n MC sufficiently large, the classical estimate of the right-hand-side of Eq. ( 12) is written as

S (ε) - 1 n MC n MC =1 log(p (n MC ) ([η ar (ε)] ; ε)) , (13)
in which the estimation, p (n MC ) ([η] ; ε), of p([η]; ε) is performed using the modification [START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF] of the classical multidimensional Gaussian kernel-density estimation method [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis[END_REF][START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF] and using additional realizations [η 1 ar (ε)], . . . , [η n MC ar (ε)].

Solution of the optimization problem.

The mapping ε → S (ε) is generally not convex on E. Since E is a finite subset of A, the optimization problem defined by Eq. ( 10) is solved using an exhaustive search over E.

Convergence with respect to the number of additional realizations

As the entropy is estimated using n MC samples, the optimal values ε opt and m opt (or m opt = m opt -1) depend of the number n MC of additional realizations generated with step 7 of Section 2.1. The convergence of the optimal values can be analyzed studying the simple convergence in R of the sequences of real numbers {ε opt (n MC )} n MC and { m opt (n MC )} n MC .

Analysis of a simple example: 3D-data around a helix (spiral)

For this simple example (similar to application 2 presented in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF]), the initial given dataset is constituted of N = 400 independent realizations of random vector X = (X 1 , X 2 , X 3 ) in dimension n = 3 and is represented by matrix [x d ] in M n,N . These given data points are concentrated in the neighborhood of a helix as shown in Fig. 2. The principal component analysis (step 2) yields ν = 3 for which error pca (ν) = 1.46 × 10 -16 . The additional realizations are computed (step 7) with n MC = 100 yielding ν sim = 40,000. The admissible finite set A of the values of ε is chosen to be {0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}. The graphs of functions ε → m L (ε) and ε → S (ε) defined on A are displayed in Fig. 3 (left and right). We have m min L = m max L = 4 and E 0 = E = {5, 6, 7, 8, 9, 10} ⊂ A. The entropy-based selection yields ε opt = 6 and consequently, m opt = 4, m opt = 3. The relative mean-square error is error pro j ( m opt ) = 3.41 × 10 -3 . For ε = ε opt , Fig. 4(a) shows function α → λ α for which λ 2 = 8.25 × 10 -2 , λ 4 = 8.06 × 10 -2 , and λ 5 = 6.88 × 10 -3 . As an example of additional realizations generated by the probabilistic learning algorithm, Fig. 4(b) displays 40,000 additional realizations of X = (X 1 , X 2 , X 3 ). These realizations stay concentrated in the neighborhood of the helix (the concentration is effectively preserved).

Applications to several databases

In this section, we present the application of the selection model proposed for several databases: Scramjet-d8 and ScramJet-d16 that correspond to simulations of a ScramJet performed by Sandia National Laboratories in Livermore with two large scale complex computational models [START_REF] Soize | Enhancing model predictability for a scramjet using probabilistic learning on manifold[END_REF], and Climate-LS-21-34 and Climate-LS-21-50 that correspond to the maxima of daily precipitations, which have been predicted by Lawrence Livermore National Laboratory with a climate computational model [START_REF] Bader | Accelerated climate modeling for energy (acme) project strategy and initial implementation plan[END_REF]. Section 5.1 deals with the ScramJet database, Section 5.2 is devoted to the climate database, Section 5.3 deals with the convergence analysis with respect to the number n MC of additional realizations, and finally, in Section 6, we present an analysis of the results that have been obtained.

The applications presented below are devoted to the identification of m opt and ε opt and not to other important analyses such as the convergence of the probability distribution of the additional 12 realizations constructed with probabilistic learning. Such analysis has been made in details in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF].

Scramjet database

The database corresponds to two different numerical simulations referenced as Scramjet-d8 and Scramjet-d16. For these two cases, random vector X has dimension n = 18, which is constituted of m w = 11 random parameters, w 1 , . . . , w 11 , of the computational model and of n q = 7 quantities of interest (QoI), q = (q 1 , . . . , q 7 ), related to outputs of the computational model.

Scramjet-d16

For Scramjet-d16, there are N = 172 independent realizations in the database, ν = 17, and error pca (ν) = 1.57 × 10 -7 . The additional realizations are computed with n MC = 150 yielding ν sim = 25,800.

Climate database

The climate database is made up of N = 130 realizations of the precipitation magnitudes for each hour (24 hours per day) of 149 consecutive days, which have been performed with a climate computational model. The specific model used is the E3SM, formerly known as Accelerated Climate Modeling for Energy (ACME) [START_REF] Bader | Accelerated climate modeling for energy (acme) project strategy and initial implementation plan[END_REF]. Consequently, 130 realizations are available for 3, 576 consecutive hours. Two sets of data are extracted from this climate database, denoted by Climate-LS-21-34 and Climate-LS-21-50,which correspond to the maximum of daily precipitations for 14 consecutive days (from 21st day to 34th day) and for 30 consecutive days (from 21st day to 50th day), respectively.

Climate-LS-21-34

For each one (indexed by ) of the 14 consecutive days from the 21st day to the 34th day, the maximum p of the precipitation is identified as well as the time t at which this maximum occurs. Random vector X has dimension n = 33 and is constituted of m w = 5 random parameters w 1 , . . . , w 5 of the computational model and n q = 28 quantities of interest, q = (q 1 , . . . , q 28 ), made up of t 1 , . . . , t

Climate-LS-21-50

For each one (indexed by ) of the 30 consecutive days from the 21st day to the 50th day, the maximum p of the precipitation is identified as well as the time t at which this maximum occurs. Random vector X has dimension n = 65 and is constituted of m w = 5 random parameters w 1 , . . . , w 5 of the computational model and n q = 60 quantities of interest, q = (q 1 , . . . , q 60 ), made up of t 1 , . . . , t 30 and p 1 , . . . , p 30 , and ν = 65 with error pca (ν) = 7.12 × 10 -16 . The additional realizations are computed with n MC = 150 yielding ν sim = 19,500. The convergence of the optimal values is analyzed (see Section 3.2.3) studying the simple convergence in R of the sequences of real numbers {ε opt (n MC )} n MC and { m opt (n MC )} n MC . Such a convergence analysis can also be presented as a function of the number ν sim = N × n MC of additional realizations. We present the convergence analysis for the ScramJet-d8 database for which N = 256. The results are similar for the other three databases. Fig. 13(a) displays the graphs of functions ν sim → m opt (ν sim ) and ν sim → ε opt (ν sim ). Fig. 13(b) displays the graph of function ν sim → S opt (ν sim ) in which S opt (ν sim ) = S (ε opt (ν sim ); ν sim ) (note that the entropy S depends not only on ε but also of the number of realizations, ν sim ). It can be seen that the convergence is obtained for ν sim = 38,400, that is to say for n MC = 150.

Remark. Figs. 6(b) and 8(b) (graph of (w 1 , w 3 ) → q 2 (w 1 , w 3 )) and 10(b), and 12(b) (graph of (w 2 , w 4 ) → q 2 (w 2 , w 4 )) are merely given as illustration. These graphs suggest the existence of an intrinsic data structure and the ability to generate realizations according to this structure (the generated data clouds exhibit clustering).

Analysis of the results

Figs. 4, 6, 8, 10, and 12 (left figures), which display the graphs of α → λ α (ε) for the five databases, show that value 0.1 proposed for parameter L (in the criterion defined by Eq. ( 1)) is well chosen. Figs. 3, 5, 7, 9, and 11 show that the variations of functions ε → m L (ε) (left figures) and ε → S (ε) (right figures) are relatively different from a database to another one. In particular, the entropy function ε → S (ε) is not monotonic and Fig. 5(b) is not similar, for instance, to Fig. 11(b). Figures 3,[START_REF] Aggarwal | Mining Text Data[END_REF][START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF][START_REF] James | An Introduction to Statistical Learning[END_REF], and 11 (right figures) also show that the value 0.2 that is proposed for parameter χ in Eq. ( 6) is well adapted in order to obtain a sufficiently robust algorithm for the five databases that have been analyzed). Following the second remark introduced in Section 3.1, these figures show that the entropy S (ε) effectively "grows on average" with m = m L (ε). Finally, as examples, Figs. 4, 6,8, 10, and 9 (right figures) display the clouds of the additional realizations generated by the probabilistic learning algorithm. These figures show the complexity of the discovered geometry of the subset around which the probability measure is concentrated. In addition, the additional realizations stay concentrated around the geometry of the clouds made up of the initial data points. Using the optimal values ε opt and m opt of ε and m, the MCMC generator proposed does not induce scattering.

Conclusions

In this paper, we have presented a selection model of the isotropic-diffusion kernel hyperparameter of the diffusion maps that is used in a probabilistic learning model on manifold recently proposed by the authors. The selection model has been tested on several databases and seems to be independent of the databases used. The probabilistic learning on manifold that has been proposed is algebraically closed with respect to its hyperparameters and which can be applied for some databases. However, this robustness that has been found must be still confirmed using another databases.

  in M n M n = set of all the square (n × n) real matrices M n,N = set of all the (n × N) real matrices N = set of all the integers 0, 1, . . . R = set of all the real numbers R n = Euclidean space of dimension n [x] T = transpose of the real matrix [x] [x] F = Frobenius norm (Hilbert Schmidt norm) of matrix [x]

  to a better construction than the one used for the probability density function [x] → p [X] ([xd]), with respect to d[x], of random matrix [X], which would be estimated by using only the Gaussian kernel-density estimation method from the given dataset represented by [x d ], because p [X ε,m ] ([x]; ε, m) also uses a second source of information relative to the geometrical structure of the given dataset. An interpretation of the difference between p [X ε,m ] ([x]; ε, m) and p [X] ([xd]) is necessary and is presented next.

  (a), we have m L (ε) = 20 and m L (ε) = 19 because λ 21 (ε)/λ 2 (ε) = 0.0028/0.0316 < 0.1 while λ 20 (ε)/λ 2 (ε) = 0.0065/0.0316 > 0.1. As an illustration, Fig. 1(b) displays the distributions of eigenvalues {λ α (ε)} α for ε ∈ {10, 15, 20, 40} for a given database (Scramjet-d16 database) and allows for highlighting the construction defined by Eq. (1).

Figure 1 :

 1 Figure 1: (a): choice of ε yielding a strong decreasing of the eigenvalues λ α (ε) for which the ranks are greater than 18. The graph shows that λ 2 (ε) = 0.0316, λ 18 (ε) = 0.0231, λ 20 (ε) = 0.0065, and λ 21 (ε) = 0.0028. In the present case, for the criterion defined by Eq. (1), we have m L (ε) = 20 and m L (ε) = 19. (b): illustration of the distributions of eigenvalues {λ α (ε)} α for ε = 10 (blue circle), = 15 (red diamond), = 20 (green cross), and = 40 (black star) for a given database (Scramjet-d16 database).

3 .

 3 Selection model for calculating an optimal value ε opt of hyperparameter ε For given ε and m, let [η] → p [H ε,m ] ([η] ; ε, m) be the probability density function of random matrix [H ε,m ] with respect to the volume element d[η] on M ν,N , for which the additional realizations [η 1 ar ], . . . , [η n MC ar ] are computed as explained in step 7 of Section 2.1. Note that ν and N of random matrix [H ε,m ] are independent of ε and m.

Figure 2 :Figure 3 :

 23 Figure2: Helix database: initial given dataset constituted of N = 400 data points for random vector X = (X 1 , X 2 , X 3 ) concentrated in the neighborhood of a helix.

Figure 4 :

 4 Figure 4: Helix database: For ε = ε opt , graph of function α → λ α (left figure) and ν sim = 40,000 additional realizations of X = (X 1 , X 2 , X 3 ).

Figure 5 :

 5 Figure 5: Scramjet-d8 database: graphs of functions ε → m L (ε) (left figure) and ε → S (ε) (right figure), defined on A. On the two figures, (i) E = A is the set of the values of ε such that m min L ≤ m L (ε) ≤ m max L , (ii) the values of m L (ε) surrounded by circle symbols correspond to the values of ε that minimize m L (ε) and that define the set E 0 ⊂ E = A, and (iii) the square symbol localizes the optimal value ε opt of ε.

ForFigure 6 :

 6 Figure 6: Scramjet-d8 database: For ε = ε opt , graph of function α → λ α (left figure) and ν sim = 38,400 additional realizations of a QoI, q 2 , as function of parameters w 1 and w 3 .

Figures 7 (

 7 left and right) display the graphs of functions ε → m L (ε) and ε → S (ε) defined on A = {10, 12, 15, 18, 20, 30, 40, 50, 60, 70, 80, 90, 100}. We have m min L = 18, m max L = 21, E 0 = {60, 70, 80, 90, 100} ⊂ E = {20, 30, 40, 50, 60, 70, 80, 90, 100} ⊂ A. The entropy-based selection yields ε opt = 20 and consequently, m opt = 21, m opt = 20, and error pro j ( m opt ) = 1.67 × 10 -2 . For ε = ε opt , Fig. 8(a) shows function α → λ α for which λ 2 = 2.62×10 -2 , λ 21 = 3.53×10 -3 , and λ 22 = 2.48 × 10 -3 . Similarly to the ScramJet-d8, Fig. 8(b) displays the 25,800 additional realizations of the QoI, q 2 , as function of w 1 and w 3 .

Figure 7 :Figure 8 :

 78 Figure 7: Scramjet-d16 database: graphs of functions ε → m L (ε) (left figure) and ε → S (ε) (right figure), defined on A. On the two figures, (i) E ⊂ A is the subset of the values of ε such that m min L ≤ m L (ε) ≤ m max L , (ii) the values of m L (ε) surrounded by circle symbols correspond to the values of ε that minimize m L (ε) and that define the set E 0 ⊂ E ⊂ A, and (iii) the square symbol localizes the optimal value ε opt of ε.

Figure 9 :Figure 10 :

 910 Figure 9: Climate-LS-21-34 database: graphs of functions ε → m L (ε) (left figure) and ε → S (ε) (right figure), defined on A. On the two figures, (i) E ⊂ A is the subset of the values of ε such that m min L ≤ m L (ε) ≤ m max L , (ii) the values of m L (ε) surrounded by circle symbols correspond to the values of ε that minimize m L (ε) and that define the set E 0 ⊂ E ⊂ A, and (iii) the square symbol localizes the optimal value ε opt of ε.

λ 2 =

 2 1.46 × 10 -2 , λ 39 = 1.47 × 10 -3 , and λ 40 = 1.43 × 10 -3 . Fig. 10(b) displays the cloud of the 19,500 additional realizations of QoI, q 32 , as function of parameters w 2 and w 4 .

Figures 11 (Figure 11 :Figure 12 : 5 . 3 .

 11111253 Figure 11: Climate-LS-21-50 database: graphs of functions ε → m L (ε) (left figure) and ε → S (ε) (right figure), defined on A. On the two figures, (i) E ⊂ A is the subset of the values of ε such that m min L ≤ m L (ε) ≤ m max L , (ii) the values of m L (ε) surrounded by circle symbols correspond to the values of ε that minimize m L (ε) and that define the set E 0 ⊂ E ⊂ A, and (iii) the square symbol localizes the optimal value ε opt of ε.

Figure 13 :

 13 Figure 13: ScramJet-d8 database: convergence analysis of the optimal solution as a function of the number ν sim = N×n MC of additional realizations. (a): graphs of functions ν sim → m opt (ν sim ) (black thin line) and ν sim → ε opt (ν sim ) (blue thick line). (b): graph of the entropy function ν sim → S opt (ν sim ) for the optimal value ε opt (ν sim ) of ε.

  7) The constructed reduced ISDE is then used for generating n MC additional realizations [z 1 ar ], . . . , [z n MC ar ] of random matrix [Z ε,m ], and therefore, for deducing the additional realizations [η 1 ar ], . . . , [η n MC ar ] of random matrix [H ε,m ]. Using step 2 yields n MC additional realizations [x 1 ar

Acknowledgments

Part of this research was supported by the ScramJet-UQ project funded under DARPA's EQUIPS Program. The authors thank the Combustion Research Facility of Sandia National Laboratory in Livermore for the use of their simulated data of the ScramJet and the Climate Modeling and Analysis Group of Lawrence Livermore National Laboratory for the use of their simulated climate data in order to test the proposed methodology. Part of this research was also supported as part of the Energy Exascale Earth System Model (E3SM) project, funded by the U.S. Department of Energy (DOE), Office of Science, Office of Biological and Environmental Research under the auspices of the U.S. DOE by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. This research used resources of the Argonne Leadership Computing Facility at Argonne National Laboratory, which is supported by the Office of Science of the U.S. DOE under contract DE-AC02-06CH11357.

Appendix A. Appendix A

In this appendix, we summarize the algorithm used in Step 7 for generating n MC additional realizations [z 1 ar ], . . . , [z n MC ar ] of random matrix [Z ε,m ], using the reduced ISDE for which the details can be found in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF]. Let M = n MC × M 0 be the positive integer in which M 0 is a positive integer greater or equal to 1. The reduced-order ISDE is solved on the finite interval R = [0 , M ∆r], in which ∆r is the sampling step of the continuous index parameter r. The integration scheme is the Störmer-Verlet scheme for which M + 1 sampling points {r , = 0, . . . , M} are used with r = ∆r. For = 0, . . . , M -1, the algorithm is written as

The quantities in these equations are defined as follows.

(i) For = 0, a given realization of the random matrix [Z ε,m 0 ] is known and is denoted by

and where [N ] is a random matrix with values in M ν,N for which its columns are N independent copies of a normalized Gaussian vector with values in R ν .

(ii) For = 0, . . . , M -1, we have

M independent random matrices with values in M ν,N and where, for all k = 1, . . . , ν and j = 1, . . . , N, the real-valued random variables {[∆W +1 ] k j } k j are independent, Gaussian, second-order, and centered such that E{[∆W +1 ] k j [∆W +1 ] k j } = ∆r δ kk δ j j . 19

(iii) We have β = f 0 ∆r /4 and [L + 1 2 ] is the M ν,m -valued random variable such that

)] in M ν,N is defined, for all k = 1, . . . , ν and for all = 1, . . . , N, by • If integer M 0 is chosen sufficiently large (such that ρ is much larger than the relaxation time of the dissipative Hamiltonian dynamical system), then [z 1 ar ], . . . , [z n MC ar ] can approximately be considered as independent realizations of random matrix [Z ε,m ].


Cesar Vargas

Nadège Varsier

Jean Schwoerer

Xavier Lagrange

Energy-efficient discovery process for mMTC applications

come

Energy-efficient discovery process for mMTC applications

Cesar Vargas Anamuro * † , Nadège Varsier * , Jean Schwoerer * and Xavier Lagrange † * Orange Labs, Meylan, France Email: {cesaraugusto.vargasanamuro, nadege.varsier, jean.schwoerer}@orange.com † IMT Atlantique/IRISA, Rennes, France Email: xavier.lagrange@imt-atlantique.fr

Abstract-Reducing energy consumption is fundamental to address the requirements of Massive Machine-Type Communication (mMTC) applications. In this sense, D2D relaying could play an important role in some scenarios. D2D discovery process is a key phase of D2D relaying since in this phase the MTC device (MTD) has to select an optimal relay. This process should ensure low energy consumption and low complexity while taking into account the specific constraints of the MTDs.

In this paper, we analyze the total (i.e., during the discovery phase and the data transmission phase) energy consumed by the MTD when it uses a D2D relaying mechanism. We propose a simple distributed discovery protocol, which allows reducing the MTD energy consumption. Our numerical and simulation results confirm that when using the proposed discovery mechanism, the total MTD energy consumption can be significantly reduced depending on the data packet size and the location of the MTD with respect to the base station.

Index Terms-Energy-efficient, mMTC, D2D discovery, relay selection

I. INTRODUCTION

The future Fifth Generation (5G) of cellular networks aims to support three major service groups: Enhanced Mobile BroadBand (eMBB), Ultra-Reliable Low Latency Communications (URLLC), and Massive Machine Type Communications (mMTC) [START_REF] Gpp | Study on new radio (NR) access technology physical layer aspects[END_REF] [START_REF] Gpp | Study on provision of low-cost Machine-Type Communications (MTC) User Equipments (UEs) based on LTE[END_REF]. Each of these services has different characteristics and requirements. eMBB can be seen as an extension of the 4G networks, which focuses mainly on the improvement of the data rate. URLLC groups the applications that need high level of reliability and low latency while mMTC services are characterized by a large number of Machine-type Devices (MTDs). The MTDs are usually low-complexity devices that sporadically transmit small data packets. Moreover, most of these devices are usually equipped with low-capacity batteries, which must keep the MTD operational for many years. Minimizing energy consumption due to the transmission of data has a positive impact on the battery lifetime of the MTDs.

The use of Device-to-Device (D2D) communications to reduce the MTD energy consumption has been proposed in the literature [START_REF]Proximity-Based Services (ProSe) (Release 14)[END_REF] [4] [START_REF] Rigazzi | Aggregation and trunking of M2M traffic via D2D connections[END_REF]. An MTD in poor coverage conditions sends its data to a nearby User Equipment (UE), which acts as a relay. The main advantage of such configuration is a better link performance since the relay can be closer to the MTD than the Base Station (BS). However, the MTD must search for a relay before transmitting its data. This phase, called the discovery phase, causes additional energy consumption, which is not negligible for small data transmissions.

The use of relays to improve the energy efficiency has been well studied [START_REF] Madan | Energy-Efficient Cooperative Relaying over Fading Channels with Simple Relay Selection[END_REF] [7] [START_REF] Li | Energy-Efficient Optimal Relay Selection in Cooperative Cellular Networks Based on Double Auction[END_REF]. Authors in [START_REF] Madan | Energy-Efficient Cooperative Relaying over Fading Channels with Simple Relay Selection[END_REF] analyze the total energy consumption (including the energy cost of acquiring channel information) based on a general relay selection approach. In [START_REF] Zhou | Energy-Efficient Cooperative Communication Based on Power Control and Selective Single-Relay in Wireless Sensor Networks[END_REF] the authors propose a relay selection scheme based on 802.11 RTC/CTS (Request-to-Send/Clear-to-Send) protocol, where an optimal relay is selected in order to maximize the energy efficiency. In [START_REF] Li | Energy-Efficient Optimal Relay Selection in Cooperative Cellular Networks Based on Double Auction[END_REF], the relay selection strategy is analyzed considering the double auction theory. However, in all these works, the total energy consumption is computed taking into account also the energy consumed by the relay. This is a reasonable approach in networks where all devices have the same characteristics (e.g., wireless sensor networks).

In this work, we consider two groups of devices with different characteristics, the MTDs (energy-constrained devices) and the UEs (non-energy constrained devices) as in [START_REF] Rigazzi | Aggregation and trunking of M2M traffic via D2D connections[END_REF]. Our objective is to analyze the energy consumed by an MTD when it uses a UE as a relay to transmit its data. From this analysis, we determine in which cases it is convenient for the MTD to use a D2D relay mechanism instead of a direct transmission to the BS.

A D2D relay mechanism is composed of two phases: the discovery phase, and the data transmission phase. If no relay is selected at the end of the discovery phase, the MTD transmits data to the BS, which may induce a large energy consumption. If a relay is selected, its location is very important because the closer the relay is to the MTD, the lower the MTD energy consumption will be during the data transmission phase. Therefore, the discovery phase is the key phase to minimize the total MTD energy consumption (i.e., the energy consumed in both phases).

In this paper we propose a discovery and relay selection process in a distributed manner (without the participation of the BS) inspired by the RTS/CTS protocol as in [START_REF] Zhou | Energy-Efficient Cooperative Communication Based on Power Control and Selective Single-Relay in Wireless Sensor Networks[END_REF] [9] [10]. This protocol is suitable for mMTC applications due to its simplicity and limited signal overhead. In [START_REF] He | Cooperative RTS/CTS MAC with relay selection in distributed wireless networks[END_REF] and [START_REF] Awang | Toward an analysis of energy consumption in multihop wireless sensor networks[END_REF], the neighbor UE with the best channel condition UE-BS is selected as a relay. While [START_REF] He | Cooperative RTS/CTS MAC with relay selection in distributed wireless networks[END_REF] shows simulation and analytical results, the results were obtained only through simulations in [START_REF] Awang | Toward an analysis of energy consumption in multihop wireless sensor networks[END_REF]. In [START_REF] Zhou | Energy-Efficient Cooperative Communication Based on Power Control and Selective Single-Relay in Wireless Sensor Networks[END_REF] the authors take the power control into account. Our protocol differs from [START_REF] Zhou | Energy-Efficient Cooperative Communication Based on Power Control and Selective Single-Relay in Wireless Sensor Networks[END_REF] [9] [START_REF] Awang | Toward an analysis of energy consumption in multihop wireless sensor networks[END_REF] as follows: (i) it selects a non-energy constrained device (UE) as relay; (ii) it selects a relay close to the MTD by limiting the discovery area; (iii) it solves the hidden node problem thanks to a feedback packet; (iv) and finally it considers dynamic data rate adaptation.

First, we analyze the relay discovery probability and the average number of slots used in the contention process. From these parameters, we derive the energy consumption in the discovery phase and in the data transmission phase. In order to minimize the MTD energy consumption, we derive the optimal contention window size and the optimal radius of the discovery area. Finally, we compare the energy consumption in cellular mode and the energy consumption in D2D mode taking into account the optimal values.

The remainder of this paper is organized as follows. Section II presents the system model. The relay discovery probability is derived in section III. In Section IV, we analyze the number of slots used in the contention process. We compare the energy consumption in cellular and D2D modes in Section V while in Section VI we present the numerical results. Finally, conclusions and future works are given in Section VII.

II. SYSTEM MODEL

A. Network Model

We consider a single cell in which UEs and MTDs are randomly distributed. The locations of UEs form a homogeneous Poisson point process (PPP) Φ u with density λ u in R 2 . All the devices (MTDs and UEs) are in down-link coverage. MTDs with an unfavorable link budget have up-link coverage at the cost of higher energy consumption. In order to reduce its energy consumption, the MTD may use a nearby UE as a relay. This approach consists of two phases:

1) Discovery phase: The MTD searches for a nearby UE, which can serve as a relay. 2) Data transmission phase: If the MTD finds a relay, it transmits its data to the selected relay using a D2D communication, then the selected relay transmits the MTD data to the BS using a cellular communication.

Otherwise, the MTD transmits its data directly to the BS using a cellular communication. In our analysis, we do not consider the energy consumed by the relay. We assume that the BS allocates dedicated resources for D2D links (i.e., there is no interference between D2D links and cellular links). In order to keep the analysis tractable, we make the following assumptions:

• The path loss is distance-based (i.e., without considering fading and shadowing effects). • Devices (MTDs and UEs) transmit at a fixed transmission power and have a bandwidth of B w Hz. • The channel is reciprocal (i.e., the channel from point A to B is the same as the channel from B to A).

B. Propagation model

There are two communication modes: cellular mode and D2D mode. In cellular mode, the MTD communicates directly with the BS while in D2D mode the MTD communicates with the relay (UE). We consider the path loss model described in [START_REF] Xing | The investigation of power control schemes for a device-to-device communication integrated into OFDMA cellular system[END_REF] which is based on the micro-urban channel models from ITU [START_REF] Itu-R | Guidelines for evaluation of radio interface technologies for IMT-Advanced[END_REF].

We consider only non-line-of-sight (NLoS) scenarios. Thus the path loss models for carrier frequency (f c ) 900 MHz are given in Table I, where d m,b and d m,u are the MTD-BS and the MTD-UE distances (in meters), respectively.

C. Discovery protocol description

The proposed protocol uses a distributed relay selection approach (i.e., the BS does not participate in the relay selection process). When the MTD wakes up to transmit its data, it first waits for the synchronization messages sent by the BS. Then, the procedure consists of the following steps:

1) The MTD measures the path loss value L m,b between itself and the BS when it reads the network configuration messages sent by the BS. Based on this value, it defines the contention window size W and the path loss threshold between itself and a nearby UE (L th ). Then it broadcasts a Request-for-Relay (RR) packet that carries the values of W and L th and waits for a response from a UE (contention process). The maximum waiting time is W time slots. 2) A UE that receives an RR packet measures the path loss between the MTD and itself (L m,u ) by comparing the received and the transmission power. Only the UEs having L m,u < L th participate in the contention process (these UEs are called relay candidates). Each relay candidate responds to the MTD with a Relay-Candidate (RC) packet. To send its RC packet, the relay candidate chooses randomly and uniformly one timeslot s ∈ [1, W ] in the window. The first relay candidate that transmits an RC packet without collision wins the contention process.

3) The MTD broadcasts the feedback packet as soon as it successfully receives an RC packet. The feedback packet serves to confirm the winning UE and to avoid the hidden node problem. The hidden node problem occurs when a relay candidate is not able to hear the RC packet sent by another one. The relay candidate that wins the contention is called the selected relay while the other relay candidates stop contending when they receive the feedback packet. If the MTD finds a relay, it establishes a D2D link with the relay to transmit its data. Otherwise, the MTD transmits directly to the BS using a cellular communication. The D2D relaying procedure is represented in Fig. 2. We consider that the RR, RC, and feedback packets have a duration of one time-slot. We remark that there is no collision in the transmission of the RR packet since only the MTD transmits in that time slot. In other words, all UEs around the MTD receive the RR packet. However, when two or more UEs send the RC packet in the same time slot, we assume that the MTD is not able to decode these packets (ignoring the capture effect). Therefore, it has to listen during the next time slot.

The fact that the MTD sends the feedback packet and its data as soon as an RC packet is successfully received, on one hand, allows the MTD to save energy since the time in the active state is reduced. On the hand, the relay candidates also save energy since, after sending the RC packet, if they do not receive the feedback packet they stop participating to the contention process as it means that there was a collision. Fig. 3 shows the packet exchange sequence when the MTD selects a relay during the discovery phase. In this network, W = 6 and there are four UEs that successfully received the RR packet. In slot 1, no UE responds (No RC). In slot 2, there are two UEs responding at the same time (collision). In slot 3, only UE-4 transmits an RC packet (single slot) and thus wins the contention process. UE-1 and UE-3 stop participating in the process in slot 3 since they do not hear the feedback in this slot. UE-2, which had chosen slot 6, stops participating in the process when it hears the feedback packet in slot 4.

III. RELAY DISCOVERY PROBABILITY

The relay discovery probability is the probability that the MTD finds a relay during the discovery phase.

Let r be the distance between the MTD and a UE, the received power can be derived as P rx = P m K d r -α d , thus the MTD-UE path loss is L m,u = r α d /K d , where P m is the MTD transmission power, K d and α d are respectively the path loss factor and the path loss exponent for D2D communications. According to the discovery protocol, a UE that receives an RR packet participates in the contention process if and only if L m,u ≤ L th , where L th is the path-loss threshold, which is a parameter sent by the MTD. Thus, a UE is a relay candidate when:

r ≤ (K d L th ) 1/α d . (1)
We introduce a new parameter R = (K d L th ) 1/α d called radius of the discovery area. Thus, we define the discovery area (gray disk in Fig. 1) as the region bounded by a circle of radius R centered on the MTD. All the UEs inside the discovery area are relay candidates.

The relay candidates choose a time slot s ∈ [1, W ], with a probability of 1/W (uniform random choice). Since the locations of UEs form a homogeneous PPP with density λ u , in each time slot s, the relay candidates form an independent thinning PPP with density λ u /W . Thus, the probability of receiving n RC responses in a time slot is

P(N = n) = (πR 2 λ u /W ) n n! exp(-πR 2 λ u /W ), (2)
A successful relay discovery occurs when at least one RC packet is successfully received (single slot) by the MTD. A single slot is a time slot where only one UE transmits an RC packet (i.e., a time slot without collision). Thus, the probability that a time slot s is a single slot is

P 1 = P(N = 1) = N u W exp - N u W , (3)
where N u = πR 2 λ u represents the mean number of UEs inside the discovery area.

The relay discovery probability P disc is equivalent to the probability of finding at least one single slot during the contention process:

P disc = 1-(1-P 1 ) W = 1-1 - N u W exp - N u W W . (4)
Taking the first derivative of P disc with respect to N u , it is easy to show that the relay discovery probability is maximized when W = N u . For W = N u , P 1 = e -1 and P disc = 1 -(0.63) W .

IV. NUMBER OF SLOTS USED IN THE CONTENTION

PROCESS

The maximum number of slots in the contention process is W . However, according to the discovery protocol, the contention process ends as soon as an RC packet is successfully received (i.e., at the first single slot). In case that none of the RC packets is correctly received, the number of slots used in the contention process is W .

Let S be a discrete random variable with range {1, 2, ..., W } that represents the number of slots used in the contention process. This number of slots is equivalent to finding the first single slot. From (3), the probability that a time slot is a single slot is P 1 = (N u /W ) exp(-N u /W ). Since P 1 is independent for each time slot, S can be modeled as a geometric distribution with parameter P 1 . Thus, the probability mass function (PMF) of S is

P(S = s|W ) = (1 -P 1 ) s-1 P 1 if 1 ≤ s ≤ W -1, (1 -P 1 ) W -1 if s = W . (5)
P(S = W |W ) = (1 -P 1 ) W -1 needs special attention because it is the sum of two probabilities, the probability that slot W is a single slot (1 -P 1 ) W -1 P 1 and the probability that none of the RC packets is correctly received (1 -P 1 ) W by the MTD.

Let S be the mean number of slots used in the contention process. It can be easily calculated from the PMF of S:

S = W -1 s=1 sP(S = s|W ) + W P(S = W |W ). (6)
Substituting ( 5) into ( 6), and simplifying we obtain

S = 1 -(1 -P 1 ) W P 1 , (7)
where

P 1 = (N u /W ) exp(-N u /W ).
The energy consumption in the discovery phase is directly proportional to S. Therefore, we can minimize the energy consumption in the discovery phase by minimizing S. Taking the first derivative of S with respect to N u , we can find the minimum number of time slots in the contention process. Thus, solving ∂S ∂Nu = 0, we deduce that the number of time slots in the discovery phase is minimal when W = N u . And we obtain

min S = e(1 -(1 -e -1 ) W ). (8)
Note that when W is large min S ≈ e.

V. ENERGY CONSUMPTION IN CELLULAR AND D2D

MODES An MTD can operate in one of the following states: Active (Tx or Rx), Idle, or Sleep. Each state has different and constant power requirements [START_REF]Cellular system support for ultra low complexity and low throughput internet of things[END_REF]. In our study, we analyze the energy consumed by the MTD when it transmits data (i.e., when it is in active state). Thus, the energy consumed by the MTD is computed as E = P m,T t m,T + P m,R t m,R , where P m,T and P m,R are the MTD power consumption (energy consumption per unit time) in Tx state and Rx state, respectively; t m,T and t m,R are the proportion of time spent in Tx state and Rx state, respectively. Note that P m,T is constant since we consider a fixed transmission power.

A. Energy consumption in cellular mode

In order to compute the energy consumed by the MTD when it transmits its data directly to the BS, we need to determine the transmission time (i.e., the time spent in Tx state) as in [START_REF] Vargas Anamuro | Simple modeling of energy consumption for D2D relay mechanism[END_REF]. Assuming that the MTD transmits data message of total size D bits, the transmission time is equal to D/C m,b , where C m,b is the data rate (bits per second) in the MTD-BS link.

For computing the data rate C m,b , we use a modified Shannon capacity formula proposed in [START_REF] Mogensen | LTE capacity compared to the shannon bound[END_REF]:

C m,b = B eff B w log 2 1 + θ m,b θ eff , (9)
where B w is the transmission bandwidth, B eff adjusts for the system bandwidth efficiency of LTE, θ m,b is the Signal-to-Noise-Ratio (SNR) for the MTD-BS link, and θ eff adjusts for the SNR implementation efficiency of LTE.

As we consider that the MTD transmits at fixed power P m the SNR expression for the MTD-BS link is then:

θ m,b = P m K c (d m,b ) -αc N 0 B w , ( 10
)
where N 0 is the noise power spectral density, d m,b is the MTD-BS distance, K c and α c are respectively the path loss factor and the path loss exponent for cellular communications. Thus, the MTD energy consumption due to the data transmission is

E cell = P m,T D C m,b , (11)
where

C m,b = B eff B w log 2 1 + PmKc(d m,b ) -αc N0Bwθeff
.

B. Total energy consumption in D2D mode

The total MTD energy consumption in D2D mode is the sum of the energy consumed in the discovery phase and the energy consumed in the data transmission phase.

Since data rate adaptation requires the knowledge of channel state information (CSI), we consider that the discovery phase also serves to obtain the CSI of the MTD-UE link. Thus, we consider a fixed data rate during the discovery phase and a dynamic data rate adaptation during the data transmission phase.

1) Energy consumption in the discovery phase: We consider that the RR and feedback packets have a duration of one time-slot (T s ). The MTD is in Tx state when it transmits the RR or the feedback packet, and it is in Rx state during the contention process. Thus, the average energy consumption during the discovery phase is

E D2D,disc = T s (2P m,T + SP m,R ), ( 12
)
where P m,T and P m,R are the MTD energy consumption in Tx state and Rx state, respectively; and S is the mean number of slot used in the contention process, which is defined in [START_REF] Zhou | Energy-Efficient Cooperative Communication Based on Power Control and Selective Single-Relay in Wireless Sensor Networks[END_REF].

2) Energy consumption in the data transmission phase: Once the discovery phase has been completed, the MTD will transmit its data to the selected relay (D2D communication). In case no relay has been found, due to low relay density or collisions, the MTD will transmit directly to the BS as in traditional cellular communications. Thus, the average energy consumption during the data transmission phase can be derived as

E data Tx = P disc E D2D,comm + (1 -P disc )E cell , (13)
where P disc is the relay discovery probability, which is defined in ( 4), E D2D,comm is the average energy consumption in the D2D communication, and E cell is the average energy consumption in cellular mode, which is defined in [START_REF] Xing | The investigation of power control schemes for a device-to-device communication integrated into OFDMA cellular system[END_REF]. E D2D,comm depends on the MTD-Relay distance, which is a random variable since the locations of the UEs form a PPP. Let X be a random variable that represents the distance between the MTD and the relay. Since we consider a uniform random choice of time slots in the contention process, all relay candidates have the same probability of being selected as a relay (without considering their locations). Thus, the cumulative distribution function (CDF) of X can be derived as

F X (x) = P(X ≤ x) = x 2 /R 2 , x ∈ [0, R], (14)
where R is the radius of the discovery area. Thus, the probability density function (PDF) of X is

f X (x) = 2x/R 2 , x ∈ [0, R]. (15)
To compute the energy consumption in the D2D communication (E D2D,comm ), we follow the same procedure used to derive [START_REF] Xing | The investigation of power control schemes for a device-to-device communication integrated into OFDMA cellular system[END_REF]. We consider that during the D2D communication the MTD transmits at fixed power P m . Using the modified Shannon capacity formula [START_REF] Mogensen | LTE capacity compared to the shannon bound[END_REF], we obtain

E D2D,comm = P m,T D B eff B w log 2 1 + PmK d x -α d N0Bwθeff , (16)
where x is the MTD-Relay distance, K d and α d are respectively the path loss factor and the path loss exponent for D2D communications. Thus, the average energy consumption in the D2D communication is computed as follows:

E D2D,comm = R 0 P m,T D B eff B w log 2 1 + PmK d x -α d N0Bwθeff f X (x)dx, (17) where f X (x) = 2x/R 2 is the PDF of X (15).
Taking into account that the MTD transmits to its maximum allowed transmission power (i.e., P m = 23 dBm) and the MTD-UE distance x is short, we can assume:

log 2 1 + P m K d x -α d N 0 B w θ eff ≈ log 2 P m K d x -α d N 0 B w θ eff . (18)
Then in (17) we have

E D2D,comm = P m,T D B eff B w R 2 R 0 2x log 2 PmK d x -α d N0Bwθeff dx, (19)
Integral ( 19) is simplified as follows:

E D2D,comm = Ae B E 1 (B), (20)
where

A = ln(4)P m,T D α d BeffBw , B = 2 α d ln PmK d R -α d

N0Bwθeff

, and E 1 (z) = ∞ z e -t t -1 dt is the exponential integral. Then, substituting (4) and (20) in [START_REF]Cellular system support for ultra low complexity and low throughput internet of things[END_REF] we get

E data Tx = (1-(1-P 1 ) W )Ae B E 1 (B)+(1-P 1 ) W E cell , (21)
where

P 1 = πR 2 λu W exp(-πR 2 λu W )
, and E cell is defined in [START_REF] Xing | The investigation of power control schemes for a device-to-device communication integrated into OFDMA cellular system[END_REF].

3) Total MTD energy consumption in D2D mode: The average of the total MTD energy consumption in D2D mode is

E D2D,total = E D2D,disc + E dataTx . (22)
From ( 12) and ( 21), we observe that E D2D,total depends on three variables: d m,b , W , and R.

4) Minimization of the MTD energy consumption:

We propose a discovery mechanism that aims to minimize E D2D,total . As discussed in sections III and IV, E D2D,disc is minimal when W = N u (since S is minimal). Moreover, when W = N u the maximum relay discovery probability (P disc ) is reached, and consequently E dataTx is minimized since in data transmission phase when the MTD uses a relay it consumes less energy than when it transmits directly to the BS. Therefore the optimal contention window size, minimizing the total MTD energy consumption is W * = N u and E D2D,total depends now only on two parameters (R and d m,b ).

Looking for the optimal radius of the discovery area that minimizes the total MTD energy consumption in the D2D mode, let R = ρd m,b , where 0 < ρ < 1. We then look for an optimal ρ * value that minimize E D2D,total when the distance between the MTD and the BS is d m,b . It is not possible to get this value analytically. We thus use numerical computation to find ρ * . Fig. 4 shows the energy consumption in D2D mode as a function of the distance MTD-BS and the parameter ρ. From this figure, we observe (in red) that for each distance MTD-BS there exists an optimal ρ * that minimizes the total MTD energy consumption. We plot in Fig. 5 the optimal ρ * as a function of the distance MTD-BS for λ u = {10 × 10 -6 , 100 × 10 -6 } UEs/m 2 and D = 200 bytes. The optimal radius of the discovery area is From (3), we know that N u = πR 2 λ u . Thus, the optimal contention window size is

R * = ρ * d m,b . (23)
W * = π(ρ * d m,b ) 2 λ u . (24)
From ( 1), we can derive the path loss threshold L th as follows:

L th = (ρ * d m,b ) α d K d . ( 25

VI. NUMERICAL RESULTS

The analytical and simulation results are presented to validate the accuracy of our analytical models. The simulations are performed using MATLAB and the simulation parameters are specified in Table II. We use the MTD power consumption parameters proposed in [START_REF]Recommended Multi-Service Performance Optimization Solutions for Improved E2E Performance[END_REF], and we set the bandwidth efficiency B eff and the SNR efficiency θ eff as the SISO channel parameters [START_REF] Mogensen | LTE capacity compared to the shannon bound[END_REF].

In Fig. 6 we show that the probability of finding a relay decreases as the mean number of UEs in the discovery area (N u ) decreases. It is either due to a low UE density (λ u ) or because the discovery area (πR 2 ) is very small. The relay discovery probability increases as N u increases, until reaching a maximum. However, due to the limited number of time slots (W ), collisions also increase, reducing the relay discovery probability. Fig. 7 shows the mean number of slots used in the contention process S as a function of N u . We can see that the analytical and simulation results fit well. From this figure, we observe that when N u is zero (i.e., there is no UEs inside the discovery area), the number of time slots used in the contention process is W . This is because even if there is no relay candidate, the MTD waits for an RC packet until the last time-slot. We also note that as N u increases, S decreases and reaches a minimum value at N u = W , and then S increases (i.e., a convex curve). This is because the probability of receiving an RC packet in the first slots increases as the number of relay candidates increases. But, if the number of relay candidates is greater than the number of slots available (W ), thus S increases due to the impact of the collisions. In Fig. 8, we compare the energy consumption in cellular mode and the minimum average energy consumption in D2D mode (i.e., considering the the optimal values W * and ρ * ). This figure shows that the D2D relay mechanism allows reducing the energy consumption significantly when the data packet size is relatively large compared to the discovery packet size (12 bytes), and when the MTD is far from the BS.

When the data packet size is small (e.g. 50 bytes), the energy consumption in cellular mode is less than the energy consumption in D2D mode, because the energy consumed in the discovery phase is significant in comparison to the energy consumed in the cellular mode.

The energy reduction factor is defined as the ratio between the energy consumption in D2D mode and the energy con- sumption in cellular mode. From Fig. 8, the energy reduction factor, for two different data sizes and for two different MTD-BS distances, are obtained and shown in Table III. We observe that in D2D mode, the MTD consumes 50% less energy than in cellular mode for D = 200 bytes and d m,b = 1000 meters.

VII. CONCLUSION AND FUTURE WORKS

In this paper, we analyze the total MTD energy consumption due to the transmission of data. We propose a simple distributed discovery protocol, which is efficient in terms of energy consumption. We show that using this D2D relay mechanism, the MTD energy consumption is significantly reduced when the data packet size is relatively large compared to the discovery packet size, and when the MTDs are far from the BS (unfavorable link budget). Furthermore, we derive analytical expressions for the relay discovery probability, the mean number of times slots used in the contention process, and the total MTD energy consumption. The accuracy of these analytical expressions is confirmed by simulations.

For future work, we want to extend our analysis of MTD energy consumption by taking into account the fading and shadowing effects. Moreover, we will investigate the MTD energy consumption using D2D relaying mechanisms considering the autonomous resource allocation mode for the discovery phase.

Fig. 1 .

 1 Fig. 1. Network model.

Fig. 2 .

 2 Fig. 2. D2D relaying procedure a) MTD side, and b) UE side.

Fig. 3 .

 3 Fig. 3. Packet exchange sequence in the D2D relaying mechanism when the MTD selects a relay during the discovery phase.

Fig. 4 .

 4 Fig. 4. Average total energy consumption in D2D mode as a function of the distance MTD-BS (d m,b ) and ρ, where R = ρd m,b (radius of the discovery area), considering D = 200 bytes and λu = 100 × 10 -6 UEs/m 2 .

Fig. 5 .

 5 Fig. 5. Optimal ρ * as a function of the distance MTD-BS considering D = 200 bytes and λu = {10 × 10 -6 , 100 × 10 -6 } UEs/m 2 .

  ) Then, the two parameters that the MTD has to send in the RR packet are: W * and L th . The MTD-BS distance (d m,b ) can be derived from the MTD-BS path-loss (L m,b ) measured by the MTD.

Fig. 6 .

 6 Fig. 6. Relay discovery probability as a function of the mean number of UEs inside the discovery area (Nu) for W={2, 4, 8, 16}.

Fig. 7 .

 7 Fig. 7. Mean number of time slots used in the contention process (S) versus the mean number of UEs inside the discovery area (Nu) for W = {16, 24, 32}.

Fig. 8 .

 8 Fig. 8. Comparison of the total MTD energy consumption in cellular mode and in D2D mode as a function of the distance MTD-BS (d m,b ) and the data size (D), considering λu = 100 × 10 -6 UEs/m 2 .

TABLE I PATH

 I LOSS MODELS FOR NLOS WITH DISTANCE IN METERS, fc = 900 MHZ

	Communication mode	Path loss (dB)
	Cellular	L

m,b = 21.53 + 36.7 log 10 (d m,b ) D2D Lm,u = 17.62 + 40 log 10 (dm,u)

TABLE III ENERGY

 III CONSUMPTION FOR CELLULAR AND D2D MODES

	Data size	100 bytes	200 bytes
	Distance MTD-BS (m)	500	1000	500	1000
	Energy consumption Cellular (mJ) 1.06	2.16	2.12	4.31
	Energy consumption D2D (mJ)	1.65	2.68	2.11	2.16
	Energy Reduction Factor	1.56	0.78	1.00	0.50

ACKNOWLEDGMENT

This paper is part of and thereby supported by the Horizon 2020 project ONE5G (ICT-760809) receiving funds from the European Union. The authors would like to acknowledge the contributions of their colleagues in the project, although the views expressed in this contribution are those of the authors and do not necessarily represent the project.