Cédric Rommel
email: cedric.rommel@safety-line.fr

Joseph Frédéric Bonnans

Baptiste Gregorutti
email: baptiste.gregorutti@safety-line.fr

Pierre Martinon
email: pierre.martinon@inria.fr

Block sparse linear models for learning structured dynamical systems in aeronautics

come

INTRODUCTION 1.Dynamical system identification problem

We consider in this paper the task of learning the state function g defining the behavior of a dynamical system: ẋ = g(t, x, u),

where x(t) ∈ X ⊂ R dx and u(t) ∈ U ⊂ R du denote respectively the state and control variables, and t denotes the time. We are mainly focused here on the dynamics of an aircraft during its climb phase. Aircraft dynamics identification is essential in several engineering application, including for the optimization of flight trajectories, which allows a more efficient use of an airline fleet regarding, for example, CO 2 emissions and fuel consumption. In this context, common state and control variables are x = (h, V, γ, m) and u = (α, N 1 ) (see table 1 for notations) and a standard model for the state equation is the following, where T, D, L, I sp are supposed to 1 depend on a subset of (x, u):

                       ḣ = V sin γ, V =
T (x, u) cos α -D(x, u) -mg sin γ m , γ = T (x, u) sin α + L(x, u) -mg cos γ mV , ṁ = -T (x, u) I sp (x, u) .

(2)

(3) (4)

(5) Despite this specific use case, which is central in our work, we believe that the techniques and framework presented hereafter could be used in many other engineering and scientific contexts where sets of differential equations need to be identified based on data.

In the application case of trajectory optimization, the model used for prediction is required to be fast to evaluate and differentiable. More particularly for aircraft trajectories, scalability is also a requirement, as lots of data are available for airlines which record thousands of variables every second of every commercial flight. We also restrict our scope to models which are easily interpretable by domain experts, which includes allowing to compute T, D, L and I sp . This means that learning the function g from ( 1) summarizes here into using equations (2)-( 5) to infer the hidden functions T, D, L and I sp , based uniquely on observations of u, x and ẋ.

We propose in this paper a method addressing all the requirements cited above, based on a structured version of multi-task sparse linear regression.

Related work

Dynamical system identification According to the literature [Jategaonkar, 2006, Maine and Iliff, 1986, Klein and Morelli, 2006], two widely used approaches for learning aircraft state equations are the Output-Error Method (OEM) and Filter-Error Method (FEM). Both are said to be dynamic identification techniques, in the sense that they evaluate the goodness of a certain candidate solution ĝ by integrating it given some fixed controls and comparing the states obtained to the measured states (or some other variables measured). The difference between them lies in the uncertainty assumed for the state equation integration, a Kalman filter being used in the latter to estimate the states instead of a deterministic integration scheme. Recent advances include using neural networks for the state estimation part [START_REF] Peyada | Aircraft parameter estimation using a new filtering technique based upon a neural network and Gauss-Newton method[END_REF].

The older Equation-Error Method (EEM) [START_REF] Maine | Identification of Dynamic Systems: Theory and Formulation[END_REF] differs from the former by trying to fit the state derivative information from the dynamic system instead of matching the states themselves. This makes this class of method more scalable and avoids several pratical issues related to the nested estimation problem of integrating the states of a possibly nonlinear dynamical system. Drawbacks of this method include needing access to state derivatives data (which are usually obtained by numerical differentiation of noisy time series) [Morelli, 2006] and lower accuracy than OEM and FEM in terms of predicted state variables [START_REF] Peyada | Aerodynamic characterization of hansa-3 aircraft using equation error, maximum likelihood and filter error methods[END_REF]. Furthermore, in the presence of noisy and highly correlated regressors, as it is our case, it has been demonstrated in [Morelli, 2006] that standard EEM produces biased estimations of the model parameters. Frequency domain variations of EEM have since been proposed to mitigate this phenomenon.

Multi-task learning Since the introduction of the multi-task learning concept, benefits have been evidenced by the statistics and machine learning communities, including in particular when a common representation is assumed to be shared by all tasks [Caruana, 1997, Evgeniou et al., 2005]. It has been addressed briefly by the aircraft dynamics identification community for EEM approaches, with the different tasks corresponding to different manoeuvers [Morelli, 2006] or to different dynamic equations [START_REF] Rommel | Aircraft dynamics identification for optimal control[END_REF].

Feature selection based on L 1 regularizer The motivation for sparse linear models in aircraft dynamics identification problems has been presented in [START_REF] Rommel | Aircraft dynamics identification for optimal control[END_REF]. Since the introduction of the Lasso estimator by [Tibshirani, 1994], many variations have been proposed to perform L 1 -norm based feature selection in several settings. For example, the extensively studied Group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] performs selection among groups of variables by making use of a L 1 -regularization over the L 2 -norms of subsets of features. Similarly, multi-task feature selection techniques, such as [START_REF] Argyriou | Convex multi-task feature learning[END_REF], Obozinski et al., 2006], make use of a mixed L 2,1 -norm to find common representations for all tasks. As for some applications sparsity within the groups is desirable, additionally to the selection of groups of variables, the Sparse-group Lasso was proposed by [START_REF] Friedman | A note on the Group Lasso and a Sparse group Lasso[END_REF].

Moreover, in some contexts, as for example when regressors are strongly correlated, the Lasso selection is known to be inconsistent. This translates in practice by the selection of irrelevant variables, which can be evidenced by obtaining different sparsity patterns when solving the Lasso several times with small changes in the training set. The Adaptive Lasso, proposed in [Zou, 2006], was proved to carry consistent feature selection in this context by making use of data-dependent weights added to the L 1 penalty. Nonetheless, this last approach needs in practice the tuning of at least an additional hyperparameter, when compared to the Lasso. Another consistent feature selection algorithm is the Bolasso [Bach, 2008], which is based on the idea of intersecting the sparsity patterns returned by the Lasso run on bootstrap replications of the training set.

Core contributions

Since OEM and FEM classes of approaches do not scale very well, we propose a variation of EEM delivering a model with all the desired properties cited in section 1.1. This is achieved by casting the system (2)-( 5) into a multitask regression problem, which is the subject of section 2. Such formulation ensures that the structure of the differential system is conserved and benefits from the tight coupling existing between equations. This idea is similar to [START_REF] Rommel | Aircraft dynamics identification for optimal control[END_REF], with the difference that the multi-task regression model proposed herein is linear. Despite this difference, we also assume a sparse polynomial structure for the functions T, D, L and I sp . This not only complies with information found in the flight mechanics literature, but also ensures interpretable and differentiable models, which are fast to evaluate.

Sparsity within the groups of variables of the models of T, D, L and I sp was ensured by the derivation of a variation of the Lasso [Tibshirani, 1994] in a structured multi-task setting. It relates in this sense to the Sparse-group Lasso [START_REF] Friedman | A note on the Group Lasso and a Sparse group Lasso[END_REF]. However, here we assume that the group sparsity pattern is known for each task and we only promote sparsity within the groups. It also differs from Multi-task Lasso [START_REF] Argyriou | Convex multi-task feature learning[END_REF], Obozinski et al., 2006], as we do not assume that all tasks share exactly the same features. Given that polynomial models present highly correlated features in practice, consistent feature selection is obtained by applying the Bolasso algorithm [Bach, 2008]. Section 3 contains our results on this subject.

Finally, we propose in section 4 an additional regularization, which injects more information into the problem. This allows to circumvent the issues related to the high correlation between the groups of variables, leading to better identifiability of the hidden elements T, D, L and I sp .

All these measures in conjunction allow to obtain predictions with sufficient accuracy for trajectory optimization purposes, which was verified empirically. The results of our experiments based on real data can be found in section 5.

STRUCTURED MULTI-TASK MODEL

Hidden functions parametric models

In order to ensure interpretability and scalability, we chose to work only with parametric models of the hidden functions T, D, L and I sp .

As we wanted to use the available expert knowledge of flight mechanics, new features were designed for each of these four elements based on the usual dependencies found in the literature [Roux, 2005, Hull, 2007]. Coincidentally, each of these four functions depend on different sets of 3 features, which can be obtained from the initial inputs of our problem (x, u) through nonlinear mappings of the form

ϕ : X × U → R 3 , (6)
where ∈ {T, D, L, I sp }. More precisely, using the notation from table 1, we have:

ϕ T (x, u) = N 1, ρ(h), M (h, V ) , ϕ D (x, u) = q(h, V ), α, M (h, V ) , ϕ L (x, u) = q(h, V ), α, M (h, V ) , ϕ Isp (x, u) = SAT (h), h, M (h, V ) . ( 7
)
The reader is referred to the supplementary material for further details on the expressions of these common flight mechanics dependencies. Furthermore, flight mechanics models found in the literature were mostly polynomials on the features obtained through (6). More precisely, many models found consisted in the product between one feature from the triplets listed in (7) and a linear combination of some monomials of the two remaining features. This led us to define the additional feature mapping:

R 3 → R r Φ d : (z 1 , z 2 , z 3 ) → z 1 z k 2 z j-k 3 ; j = 0, . . . , d k = 0, . . . , j , (8)
with d ∈ N * and r = d + 2 2 . By picking d > 1, this last transformation allows to project our initial features into a higher dimensional space, potentially increasing the explanatory power of our model. Hence, for some d T , d D , d L , d Isp > 1, we assume the following linear structures:

T (x, u, θ T ) = X T • θ T , D(x, u, θ D ) = X D • θ D , L(x, u, θ L ) = X L • θ L , I sp (x, u, θ Isp ) = X Isp • θ Isp , (9)
where

X = Φ d (ϕ (x, u)) ∈ R r , (10)
for ∈ {T, D, L, I sp }, and θ T , θ D , θ L , θ Isp denote vectors of parameters. Note that, as Φ d and ϕ are differentiable for all d ∈ N, it is also the case for models (9) as required.

Multi-task regression framework

In order to cast our identification problem into a regression problem, we first rewrite equations (3)-( 5) as follows, to isolate the functions to be learned in the right-hand side of the system:

     m V + mg sin γ = T (x, u) cos α -D(x, u), mV γ + mg cos γ = T (x, u) sin α + L(x, u), 0 = T (x, u) + ṁI sp (x, u). ( 11
) (12) (13)
This led us to define the following set of regression models, after injecting T, D, L and I sp linear structures previously derived in (9):

     Y 1 = X T 1 • θ T -X D • θ D + ε 1 , Y 2 = X T 2 • θ T + X L • θ L + ε 2 , 0 = X T • θ T + X Ispm • θ Isp + ε 3 , (14) (15)
(16)

where ε 1 , ε 2 , ε 3 , are random errors of mean 0 and

Y 1 = m V + mg sin γ, Y 2 = mV γ + mg cos γ, X T 1 = X T cos α, X T 2 = X T sin α, X Ispm = ṁX Isp . (17)
As commonly experienced in polynomial regression, the groups of interaction features X T 1 , X T 2 , X T , X D , X L and X Ispm are highly correlated in practice. This causes numerical difficulties in solving regression problems ( 14)-( 16), the parameter vectors θ T , θ D , θ L and θ Isp being able to compensate for each other.

In addition, we want to avoid the trivial solution consisting in θT = θIspm = 0 for problem (16) treated alone. Finally, the approach consisting in solving the three problems separately (e.g. through Ordinary Least Squares) would lead to three different estimators of T , inducing dynamics which are inconsistent with the known physics.

Because of the facts cited above, we chose to solve the three problems jointly in a multi-task regression framework:

Y = Xθ + ε, (18)
where

Y = (Y 1 , Y 2 , 0) ∈ R 3 , ε = (ε 1 , ε 2 , ε 3 ) ∈ R 3 , θ = (θ T , θ D , θ L , θ Isp ) ∈ R p
and X is a linear operator represented by the following sparse random matrix:

X =    X T 1 -X D 0 0 X T 2 0 X L 0 X T 0 0 X Ispm    ∈ R 3×p . ( 19
)
By doing so, we enforced all three components of the final estimator to share the same T function, which complies with the dynamics structure (2)-( 5). We also hoped to leverage the resulting coupling between tasks to increase the predictive accuracy, as observed in many other multi-task learning applications [Caruana, 1997, Evgeniou et al., 2005]. Furthermore, such framework helps to identify parameters θ T , θ D , θ L , θ Isp in a high correlations setting, as explained with more details in subsection 4.

BLOCK SPARSE ESTIMATORS

Block sparse Lasso

We consider in this section that we have access to N random observations of x, u, ẋ, allowing to build through the feature mappings (10), ( 17) and ( 19) a training set {(X i , Y i )} N i=1 , where the X i 's are sampled from the distribution of the random matrix (19). We assume this sample to be i.i.d and that all these random variables are centered and scaled to 1-standard deviation, except for the third component of the targets Y i 3 , which is assumed to be sampled from δ 0 a Dirac distribution in 0. In this context, the empirical risk minimization for the squared loss writes:

min θ N i=1 Y i -X i θ 2 2 , (20)
which is equivalent to

min θ T ,θ D ,θ L ,θ Isp Y 1 -X T 1 θ T + X D θ D 2 2 + Y 2 -X T 2 θ T -X L θ L 2 2 + X T θ T + X Ispm θ Isp 2 2 , (21)
where X , ∈ {T 1, T 2, T, D, L, I spm }, (resp. Y 1 , Y 2 ) denote matrices whose N rows correspond to the observations of X (resp. observations of Y 1 , Y 2 ).

Considering that the models from the literature used as inspiration for building our linear structures (9) were often composed of only a few monomials, feature selection seemed necessary here to try to avoid overfitting with an excessively complex model. Thus, the addition of an L 1 -penalization to problem (20) was considered:

min θ N i=1 Y i -X i θ 2 2 + λ θ 1 , λ > 0. (22)
Problem ( 22) can be seen as a structured multi-task version of the Lasso [Tibshirani, 1994], called herein the Block sparse Lasso, where the residuals are in R 3 and the features vector has been replaced by the sparse features matrix (19). Despite this differences, we can solve ( 22) with all the available algorithms for solving the Lasso, by noticing that the three terms forming the criterion from (21) can be written as a single-task least-squares criterion

Y -Xθ 2 2 , (23)
where the targets observations have been stacked Y = (Y 1 1 , . . . , Y N 1 , Y 1 2 , . . . Y N 2 , 0, . . . , 0) and the features observations have been used to build the block-sparse design matrix which gave the method its name:

X =    X T 1 -X D 0 0 X T 2 0 X L 0 X T 0 0 X Ispm    ∈ R 3N ×p . (24)
Notice that, although X has a structure similar to the random matrix X from (19), it is not random and it has 3N rows instead of 3.

Bootstrap implementation

As it is the case for the Lasso (see e.g. [Zou, 2006]), it happens in practice that problem (22) leads to unstable solutions in a strong correlations setting. Hence, we propose to use resampling in order to ensure consistent feature selection, as presented in [Bach, 2008]. We make the following standard assumptions on the joint distribution P X,Y of (X, Y ): (A1) The cumulant generating functions E exp(s X 2 2 ) and E exp(s Y 2 2 ) are finite for some s > 0. (A2) The joint matrix of second order moments

Q = E XX ∈ R p×p is in- vertible. (A3) E [Y |X] = X • θ and V [Y |X] = σ 2 a.s. for some θ ∈ R p and σ ∈ R * + .
In our setting, assumption (A1) is equivalent to requiring finite cumulant generating functions for Y 1 , Y 2 , X , ∈ {T 1, T 2, T, D, L, I spm }. This is the case when ε 1 , ε 2 , X , have compact support, which is verified in our case due to the physical nature of our data. Assumption (A2) summarizes here to each non-zero element of X having invertible second order moments, which has to be verified for our model to be identifiable. Assumption (A3) is of course equivalent to assuming the regression model ( 18). We denote by J * = {j, θ j = 0} the sparsity pattern of θ.

The feature selection procedure proposed for our problem is described in algorithm 1.

Algorithm 1 Block sparse Bolasso

Require:

training data T = {(X i , Y i )} N i=1 , number of bootstrap replicates m, L 1 penalization parameter λ, for k = 1 to m do Generate bootstrap sample T k , Compute Block sparse Lasso estimate θk from T k , Compute support J k = {j, θk j = 0}, end for Compute J = m k=1 J k , Compute θJ from T J = {(X i J , Y i )} N i=1 using (20).
In the previous algorithm, X i J denotes the matrix formed by the columns of X i indexed by J.

As the Block sparse Lasso is equivalent to the Lasso with observations stacked in a certain way, the Bolasso consistency results summarized in Theorem 3.1 and proven in [Bach, 2008], also apply to algorithm 1: Theorem 3.1 (Bolasso consistency) Assume (A1-3) and λ = λ 0 N -1 2 , λ 0 > 0. Then the probability that the Bolasso does not exactly select the correct model, i.e., for all m > 0, P [J = J * ] has the following upper bound:

P [J = J * ] ≤ mA 1 e -A2N + A 3 log N N 1/2 + A 4 log m m , ( 25
)
where

A 1 , A 2 , A 3 , A 4 > 0.
The main drawback of this general approach seems to be that it is based on solving a potentially large number of Block sparse Lasso problems, which makes it computationally costly. This flaw may be mitigated by the use of existing efficient algorithms to solve the standard Lasso problem, such as the modified Least Angle Regression algorithm [Efron et al., 2004], whose computational complexity is equivalent to a matrix inversion.

REGULARIZATION FOR ENHANCED IDEN-TIFIABILITY

With the parameters θ obtained by solving any of the optimization problems ( 20) and ( 22), one can use ( 9) to have access to predictions of the hidden functions T , D, L and Îsp . In practice, least-squares solutions produced by ( 20) lead to predicted thrust and specific impulse T , Îsp really small compared to the known order of magnitude of these physical quantities. Moreover, when using sparsity promoting formulations such as the Block sparse Lasso ( 22) or Bolasso, T and Îsp are systematically equal to 0, all their parameters being rejected in the feature selection. This seems to occur because we are regressing a function constantly equal to 0 in the last task ( 16), whose trivial solution is indeed setting all parameters of T and I sp models to 0. This seems to be allowed by our model because the group of features of T are highly correlated to the features of D and L. Hence, it is possible to find solutions to problems ( 20) and ( 22) where targets Y 1 and Y 2 are completely explained by D and L, their parameters having compensated the absence of T .

To avoid this behavior, we introduce an additional regularization to our model. For this, we assume the availability of some prior estimator for any one of the hidden functions T, D, L, I sp . This can be reduced to a constant estimator, equal to the known order of magnitude of given quantity. Without loss of generality, we suppose in what follows that Ĩsp is such a prior estimator of I sp and we define the regularized Block sparse Lasso problem:

min θ N i=1 Y i -X i θ 2 2 + λ 2 Ĩi sp -X i Isp • θ Isp 2 2 +λ 1 θ 1 , λ 1 , λ 2 > 0, (26)
where, for i = 1, . . . , N , we denote Ĩi sp = Ĩsp (x i , u i ) the prediction of the prior estimator for some sample observations of states and controls.

Up to a scaling factor λ 2 , the additional L 2 regularization term in ( 26) can be interpreted as if we had added one more task to system ( 14)-( 16):

Ĩsp = X Isp • θ Isp + ε 4 , (27)
where ε 4 is a random error with 0 mean. Hence, algorithm 1 can be run with the regularized Block sparse Lasso (26) instead of ( 22) by replacing

X i , Y i by Xi =      (X i T 1 ) -(X i D ) 0 0 (X i T 2 ) 0 (X i L ) 0 (X i T ) 0 0 (X i Ispm ) 0 0 0 λ 2 (X i Isp )      , (28)
and

Ỹ i = (Y i 1 , Y i 2 , 0, λ 2 Ĩi sp ).
Empirical results presented in section 5 show that predictions of the hidden functions with the correct orders of magnitude can be obtained with this technique.

EXPERIMENTS

Data set description and preprocessing

QAR data description We present in this section numerical experiments carried using real flight data, extracted from Quick Access Recorder devices (QAR) of 25 medium haul aircraft of the same type. Although the original data set contains raw measurements of thousands of different variables, only the following were used herein: (h, M, C, N 1 , SAT, β). The reader is referred to table 1 for the notations. Only data concerning the climb phase of the flights were kept, i.e. data corresponding to altitudes between FL50 = 5 000 f t and the top of climb (cruise altitude), specific to each flight. This was justified by the fact that the estimated dynamics would be used for the optimization of the climb profiles of these aircraft. The obtained data set contains 8 261 619 observations, made of 10 471 different flights sampled at 1 measurement per second.

Derivation of state and control observations

The raw data was smoothed using univariate smoothing splines and then used to derive all the state and control variables x, u presented in section 1.1. This was done by the means of standard flight mechanics formulas, which are listed in the supplementary material. Observations of state derivatives ẋ were then computed analytically, based on the previous formulas and on the splines-derivatives of the smoothed measured variables. The features used to build the linear models of T, D, L and I sp were obtained using the feature mapping ( 10

Experiments design

Feature selection In order to assess the quality of the feature selection performed by the Block sparse Bolasso presented in section 3.2, algorithm 1 was run seperately on the datasets of each aircraft. The sparsity patterns obtained were then compared. The regularization parameter λ 2 introduced in section 4 was set to 200, which is justified in section 5.3. The prior model Ĩsp used here to define such regularization was found in [Roux, 2005]. The L 1 penalty parameter λ 1 of the regularized Block sparse Lasso (26) was set using 30-fold cross-validation based on the squared-loss. This was performed separately for each dataset, on 33%-validation sets, and done a single time, prior to the m replications of the Bolasso. For all experiments, the number of bootstrap replications was set to m = 128. Solving the mutiple Block sparse Lasso problems was done using the Least Angle Regression algorithm, implemented in Python's scikit-learn.linear model library [Pedregosa et al., 2011].

Quality of the state equation estimations

The quality of the estimators obtained using the Block sparse Bolasso algorithm was assessed using a subset of q = 424 flights, corresponding to a single aircraft and comprising 334 531 observations. Our model was trained on q -1 flights, leaving out one flight (chosen randomly) for testing. The L 1 parameter was selected using cross validation, prior to the execution of algorithm 1, as previously explained for the feature selection assessment. The calibration of λ 2 was done by solving (26) with varying parameter values. We kept the parameter leading to predicted T , D, L and Îsp agreeing with the physical context.

Let N test denotes the number of observations from the test flight. As a first quality criterion, state derivatives of the test flight are predicted using our trained model and are then compared to the observed states derivatives of given flight. Considering that the goal of our approach is to deliver models that are suitable for trajectory optimization, another quality criterion was considered, based on comparing the observed state and control variables of the test flight {x i m , u i m } Ntest i=1 to the solution of the following trajectory optimization problem:

min u∈U t f 0 u(t) -u m (t) 2 u + x -x m (t) 2 x dt s.t. ẋ = ĝ(x, u), (29)
where t f denotes the time horizon of the test flight and • u , • x denote scaling norms. Note that the aircraft dynamics in this optimal control problem (29) were defined with the trained state equation model. Solutions were obtained using the software BOCOP [START_REF] Bonnans | Bocop -A collection of examples[END_REF].

Results

The results of the feature selection performed on the 25 data sets are on figure 1. Each row of these matrices corresponds to a different aircraft and each column to a different feature. The cells colors encode the frequency of selection of given feature for given aircraft across all the Block sparse Lasso executions. It can be observed that most dark columns are quite homogeneous, which indicates that similar features were selected for the majority of aircraft. This seems to validate our approach for this kind of data, as one would expect that airplanes of the same type should have physical models for the thrust, drag, lift and specific impulse with similar structures. As an example, a common sparse model for the thrust force T would be here made of the features:

X T = N 1 (1, ρM 2 , M 4 , ρ 2 M 2 , ρ 3 M, ρ 4 ).
Figure 2 illustrates the effect of the regularization presented in section 4 on the test flight predictions obtained with our estimator. Left frames show that it affects mainly T and Îsp , which mirror each other because of equation ( 16). The drag D seems to vary slightly when the value of λ 2 is changed, while L is unchanged. It is well-known that during climb, the thrust force of commercial aircraft decreases. However, it may be observed that for a small value of λ 2 = 20, T and Îsp are flat, which seems to indicate that more regularization is needed. A high value of λ 2 = 2000 brought our predicted thrust to increase quickly at the end of the climb, which does not seem correct either. We seem to have constrained excessively our model of Îsp . Hence, λ 2 = 200 seems like the best choice here. Right frames of figure 2 show that the regularization parameter λ 2 has little effect on the final predictions of the state function estimator ĝ. These predictions seem to be relatively good approximations of the observed states derivatives for our purposes. This can be seen for instance in figure 3, which shows the obtained solution for problem (29), illustrating that it is possible to resimulate the test flight using our estimator.

CONCLUSIONS

We developed in this paper a new method for aircraft dynamics identification based on a multi-task learning framework. As it has been designed for trajectory optimization purposes, it delivers interpretable models, which agree with the flight mechanics common knowledge and can be quickly evaluated and differentiated by an optimization solver. Our estimator is able to obtain consistent and structured sparsity, thanks to an original regression formulation and the use of a previously suggested bootstrap strategy: the Bolasso. We also showed that such estimator could be trained using available, well-known and efficient algorithms. The numerical experiments that we carried using real data from 25 different aircraft strongly suggest that our approach performs reliable feature selection. Furthermore, although this was not our main design criterion, satisfactory accuracy was evidenced by our results, at the cost of less than 10 minutes training using a laptop and a data set of 334 531 observations. The main drawback of our approach seems to be the need of a prior estimator of one of the hidden functions.

An interesting question raised by this approach is whether the trajectories optimized using these types of estimators lie close to the observed flights used to compute the model. This could give an indication on the reliability of such trajectories, which is currently under study. Moreover, considering the many issues caused by the strong correlations among features of T, D, L and I sp , the use of different models (e.g. using orthogonal polynomials) should also be investigated.

  ), where d T = 4 and d D = d L = d Isp = 3.

Figure 1 :

 1 Figure1: Feature selection results for the thrust, drag, lift and specific impulse models. The columns correspond to possible features for T /N 1 , D/q, L/q and I sp /SAT .

Figure 2 :

 2 Figure 2: (Left) Thrust, Drag, Lift and Specific impulse predicted for the test flight with different regularization parameters λ 2 . (Right) Respective states derivatives predictions.

Figure 3 :

 3 Figure 3: Resimulated trajectory obtained by solving problem (29) with BO-COP.

Table 1 :

 1 Variables nomenclature

	Notation Meaning
	h	Aircraft altitude
	V	Aircraft true airspeed (TAS)
	γ	Path angle
	m	Aircraft mass
	α	Angle of attack (AOA)
	N 1	Engines turbofan speed
	T	Total thrust force
	D, L	Drag and lift forces
	I sp	Specific impulse
	ρ	Air density
	M	Aircraft Mach number
	SAT	Static air temperature
	β	Pitch angle
	C	Fuel flow


Cédric Rommel
email: cedric.rommel@inria.fr

Frédéric Bonnans
email: frederic.bonnans@inria.fr

Baptiste Gregorutti
email: baptiste.gregorutti@safety-line.fr

Pierre Martinon
email: pierre.martinon@inria.fr

Quantifying the Closeness to a Set of Random Curves via the Mean Marginal Likelihood

Keywords: density estimation, functional data analysis, trajectory discrimination

In this paper, we tackle the problem of quantifying the closeness of a newly observed curve to a given sample of random functions, supposed to have been sampled from the same distribution. We define a probabilistic criterion for such a purpose, based on the marginal density functions of an underlying random process. For practical applications, a class of estimators based on the aggregation of multivariate density estimators is introduced and proved to be consistent. We illustrate the effectiveness of our estimators, as well as the practical usefulness of the proposed criterion, by applying our method to a dataset of real aircraft trajectories.

Introduction

Functional Data Analysis (FDA) has received an increasing amount of attention in the last years (e.g. [START_REF] Ramsay | Applied functional data analysis: methods and case studies[END_REF], this kind of data being present in many fields of application, such as speech recognition [START_REF] Ferraty | Curves discrimination: a nonparametric functional approach[END_REF], radar waveforms classification [START_REF] Dabo-Niang | On the using of modal curves for radar waveforms classification[END_REF] and aircraft trajectories classification [START_REF] Nicol | Functional principal component analysis of aircraft trajectories[END_REF][START_REF] Gregorutti | Grouped variable importance with random forests and application to multiple functional data analysis[END_REF]. In this paper we are interested in the general problem of quantifying how close some newly observed random curve is to a set of random functions, being mainly motivated by the practical task of assessing optimized aircraft trajectories. To our knowledge, this problem has not been studied in the literature.

We choose to adopt a probabilistic point of view, interpreting the original problem as the estimation of the likelihood of observing the new curve, given the sample of previously observed functions. This problem is hence related to the estimation of the probability density of a random variable valued on a function space.

Density estimation has been a longstanding problem in statistics and machine learning. Many parametric and nonparametric techniques have been proposed ever since to address it in a finite-dimensional setting. For functional data, density estimation has been studied for example by [START_REF] Dabo-Niang | Kernel density estimator in an infinite-dimensional space with a rate of convergence in the case of diffusion process[END_REF], who proposed an extension of the well-known kernel density estimator. Similarly, Prakasa Rao (2010a) developed a delta-sequence method for functional density estimation.

As the distribution of a functional random variable is hard to grasp and does not present good topological properties because it is defined on sets of a space which is "too large" (see e.g. [START_REF] Jacod | Lecture notes on "Mouvement brownien et calcul stochastique[END_REF]; [START_REF] Bosq | Nonparametric statistics for stochastic processes: estimation and prediction[END_REF]), alternatives were proposed for casting this problem into a finite-dimensional setting. For example, Prakasa Rao (2010b) proposed to project the random curves on basis functions, while [START_REF] Hall | Estimating and depicting the structure of a distribution of random functions[END_REF] suggested to study the structure of the distribution of a functional random variable by estimating its modes and density ascent lines. We propose another finite-dimensional approach, by estimating an aggregation of the marginal densities, which reduces to a finite sequence of multivariate density estimation problems.

Core Contribution. After introducing our method in section 2.1, an empirical version of it is presented for practical applications (section 2.2). The obtained statistic is built using marginal density estimators which are shown to be consistent in section 2.3. In section 3, we propose an implementation of our approach using the self-consistent kernel density estimator from [START_REF] Bernacchia | Self-consistent method for density estimation[END_REF] and extending it to the functional data context. We illustrate the effectiveness of our method, as well as its usefulness as an exploratory analysis tool for functional data, on a dataset of real aircraft trajectories and compare it to more standard approaches (section 4).

Mean Marginal Likelihood Estimator

Mean Marginal Likelihood

Let (Ω, F, P) be a probability space and T = [0; t f ] be an interval of R. We denote by E a compact subset of R d , d ∈ N * , endowed with the Borel σ-field B. Let Z = (Z t ) t∈T be a random variable valued in C(T, E), the set of continuous functions from T to E, and suppose that a training set of m observations of Z is available: T = {z 1 , . . . , z m } ⊂ C(T, E). We denote by µ t the marginal distribution of Z t for any t ∈ T, and we assume that it has a density f t relative to the Lebesgue measure on R d . We assume that (t, z) ∈ T × E → f t (z) is continuous. Let y ∈ C(T, E) be some arbitrary new curve that we would like to assess.

Given t ∈ T, we can interpret the quantity f t (y(t)) as the likelihood of observing Z t = y(t). By summarizing in some way the infinite collection {f t (y(t)) : t ∈ T}, which we call the marginal likelihoods of y hereafter, we hope to build a global and simple likelihood indicator. The first idea for aggregating these quantities is to average them with respect to time:

1 t f t f 0 f t (y(t))dt. ( 1
)
The main problem with this criterion is that it mixes elements from densities which may have very different shapes. Indeed, density values of likely observations at two times t 1 , t 2 ∈ T may have completely different orders of magnitude. For this reason, we propose to use some continuous scaling map ψ : L 1 (E, R + ) × E → [0; 1] prior to averaging:

MML(Z, y) = 1 t f t f 0 ψ [f t , y(t)] dt, (2)
and we call the obtained quantity the mean marginal likelihood of y given Z.

A natural choice for ψ is simply the function that normalizes f t (y(t)) over E:

ψ[f t , y(t)] = f t (y(t)) max z∈E f t (z) . (3)
However, we can also consider more meaningful scaling maps, such as the confidence level at y(t):

Definition 1. Let X be a continuous random variable on E of density function f ∈ C(E) and let a ∈ R + . We call the confidence level of f at a the probability that X lies in a region where its density is lower or equal to a:

ψ[f, a] = E f (x)1 {f (x)≤f (a)} dx = P (f (X) ≤ f (a)) . ( 4
)
In this case, ψ[f t , y(t)] corresponds to the probability of Z t falling outside the smallest confidence region containing y(t), as illustrated in figure 1. A numerical comparison of these two scalings can be found in section 4, while a class of estimators of the mean marginal likelihood is presented in the next section.

Empirical Version

Usually in FDA, one only has access to discrete observations of the random functions under study. We assume to be in this context: for 1 ≤ r ≤ m, each path z r of the training set T is assumed to be observed at n ∈ N * discrete times {t r 1 < t r 2 < • • • < t r n } ⊂ T, drawn independently from some random variable T , supposed independent of Z. Hence, we denote by T D the set of all discrete observations:

T D = {(t r j , z r j )} 1≤j≤n 1≤r≤m ⊂ T × E, (5)
where z r j = z(t r j ). Likewise, we assume that the new curve y is observed at ñ ∈ N discrete times { tj } ñ j=1 ⊂ T and we denote these observations by

Y = {( tj , y j )} ñ j=1 ⊂ T × E, (6)
where y j = y( tj ).

Had we enough observations of y, we could approximate the integral in (2) by a Riemann sum:

1 t f ñ j=1 ψ[f j , y j ]∆ tj , (7)
where f j := f tj , ∆ tj := tjtj-1 and t0 = 0. Yet, the marginal densities {f j } ñ j=1 are unknown and need to be estimated for practical use.

Our approach is based on the idea of partitioning T into q m intervals, or bins, of same length b m = t f /q m . For 1 ≤ ≤ q m , let τ := τ 0 + b m , where τ 0 = inf T = 0. We denote by B := [τ -1 ; τ ) the th bin for = 1, . . . , q m -1 and B qm := [τ qm-1 ; τ qm ]. Similarly, T = {(t r j , z r j ) : t r j ∈ B } ⊂ T D denotes the set of observations whose sampling time fall into B . For some m and 1 ≤ j ≤ ñ, let be such that tj ∈ B . For b m small enough, we estimate f j by building a density estimator with the partial data contained in T . This is done by applying a common statistic Θ :

S → L 1 (E, R + ) to T , where S{(z k ) N k=1 ∈ E N : N ∈ N * } denotes the set of finite sequences valued on E ⊂ R d : fj := Θ[T ]
. Hence, we consider a single density estimator per bin, averaging along the times in it. We denote this estimated quantities { fj } ñ j=1 and by summing them we can build the following plug-in estimator, called the Empirical Mean Marginal Likelihood hereafter:

EMML m (Z, y) := 1 t f ñ j=1 ψ[ fj , y j ]∆ tj . (8)
In the following subsection 2.3, sufficient conditions are given for the consistent estimation of the marginal densities f t , while section 3 presents a possible class of kernel density estimators to compute { fj } ñ j=1 .

Consistency of the marginal density estimations

General case. In this section we state that by using some well chosen statistic to build density estimators in the bins described in section 2.2 we obtain pointwise consistent estimations of the marginal densities of Z. We describe the main ideas of the proof here, while the technical details can be found in the supplementary material. Our consistency result is summarized in theorem 1 and relies on the following 4 assumptions:

Assumption 1. The random variable T is absolutely continuous and ν ∈ L ∞ (E, R + ), its density relative to the Lebesgue measure, satisfies:

ν + := ess sup t∈T ν(t) < ∞, ν -:= ess inf t∈T ν(t) > 0. ( 9
)
Assumption 2. The function defined by

(t, z) ∈ T × E → f t (z) (10)
is continuous on both variables and Lipschitz in time with constant L > 0: for any z ∈ E and

t 1 , t 2 ∈ T |f t1 (z) -f t2 (z)| ≤ L|t 1 -t 2 |. ( 11
) Assumption 3. The homogeneous partition {B m } qm =1 of T = [0; t f ],
where the bins have size b m := t f /q m , is such that

lim m→∞ b m = 0, ( 12
)
lim m→∞ mb m = ∞. ( 13
) Let S = {(z k ) N k=1 ∈ E N : N ∈ N * }
be the set of finite sequences with values in the compact set E ⊂ R d . We also need to assume that the statistic Θ : S → L 1 (E, R + ) used to build the density estimators leads to uniformly consistent density estimations in a standard i.i.d setting, which is summarized in the following assumption: Assumption 4. Let G be an arbitrary family of probability density functions on E.Given a density ρ ∈ G, let S N ρ be an i.i.d sample of size N valued in S. The estimator obtained by applying Θ to S N ρ , denoted by

ρN := Θ[S N ρ ] ∈ L 1 (E, R + ), (14)
is a (pointwise) consistent density estimator, uniformly in ρ:

For all z ∈ E, ε > 0, α 1 > 0, there is N ε,α1 > 0 such that, for any ρ ∈ G, N ≥ N ε,α1 ⇒ P ρN (z) -ρ(z) < ε > 1 -α 1 . ( 15
)
For m ∈ N * , let m : T → N * be the function mapping any point t ∈ T = [0; t f ] to the index of the bin containing it:

m (t) := t b m . ( 16
)
We denote by f m m (t) the estimator obtained by applying Θ to the subset of data points T m m (t) whose sampling times fall in the bin containing t.

Theorem 1. Under assumptions 1 to 4, for any z ∈ E and t ∈ T, f m m (t) (z) consistently approximates the marginal density f t (z) as the number of curves m grows:

∀ε > 0, lim m→∞ P | f m m (t) (z) -f t (z)| < ε = 1. ( 17
)
Remark 1. Note that, unlike assumption 4, the convergence in theorem 1 is written in terms of m. This is a big difference since the number of observation points used is supposed to be controlled in the general setting of assumption 4, while it is a random variable in theorem 1 (number of observations falling in the bin T m m (t) ).

Before explaining the proof of such a theorem, notice that the observations falling into a certain bin for a given number of curves m follow some distribution whose density function can be explicitly derived. Indeed, for V ⊂ E and B ⊂ T two compact sets, we have

P (Z T ∈ V|T ∈ B) = P ({Z T ∈ V} {T ∈ B}) P (T ∈ B) = V B f t (z)ν(t)dt dz B ν(t)dt . ( 18
)
This shows that (Z

T |Z T ∈ T m m (t) ) = (Z T |T ∈ B m m (t) ) follows a distribution of density f m m (t) (z) := τ m m (t) τ m m (t)-1 f t (z)ν(t)dt τ m m (t) τ m m (t)-1 ν(t)dt . ( 19
)
The proof of theorem 1 relies on the fact that f m m (t) converges pointwise to the marginal density f t as m tends to infinity. It is indeed quite straightforward to show this by using the assumptions that f t is Lipschitz in time (11) and that the bin sizes b m tend to 0 (12). From there, the idea is to try to apply the consistency result from assumption 4 to show that f m m (t) converges pointwise in probability to f m m (t) . However, two main difficulties arise here:

1. f m m (t) is trained using the observations from T m m (t) and the number of elements contained in this subset, denoted by N m m (t) , is random; 2. we need to train f on i.i.d observations whose number tend to infinity in order to apply (15). The first difficulty can be tackled by conditioning on N m m (t) . For the second one, we use the fact that, as the bin size tend to 0 and as the number n of observations per curve is fixed with respect to m, than each training subset has, with high probability, at most one observation per curve asymptotically. Hence, because the curves are independent observations of Z, the observations contained in T m m (t) for m large enough will be independently drawn from f m m (t) with probability 1. Furthermore, we can show that if the bin size does not decrease too fast, as required by (13), than N m m (t) diverges to +∞ in probability. The detailed proof of theorem 1 can be found in the supplementary material.

Example of a kernel estimator with deterministic kernel. In this paragraph we state a stronger consistency result for this particular setting.

f m m (t) (z) = 1 σN m m (t) z k ∈T m m (t) K z k -z σ = 1 N m m (t) z k ∈T m m (t) K σ (z k -z). ( 20
)
where K : R d → R is symmetric kernel summing to 1 and σ > 0 is the bandwidth, chosen to be scalar here for simplicity. More details on this type of estimators are given in section 3.

We denote the second moments of the random variables of density K σ and K 2 σ respectively by

σ 2 Kσ = w 2 K σ (w)dw = σ 2 w 2 K(w)dw = σ 2 σ 2 K , (21)
σ 2 K 2 σ = w 2 K σ (w) 2 dw = σ w 2 K(w) 2 dw = σσ 2 K 2 , (22)
and we denote the kernel risk by

R(K σ ) = K σ (w) 2 dw = 1 σ K(w) 2 dw = 1 σ R(K). (23)
For this particular setting, we state in theorem 2 that, under certain conditions, f m m (t) (z) approximates f t (z) consistently in expected squared-error, which is stronger than the convergence in probability stated in theorem 1:

Assumption 5. The function (t, z) ∈ T × E → f t (z) is C 4 (E) in z and C 1 (T) in t ; the Lipschitz constant of the function t → d 2 f t dz 2 (z) := f t (z) (24)
is denoted by L > 0: for any z ∈ E and t 1 , t 2 ∈ T,

|f t1 (z) -f t2 (z)| ≤ L |t 1 -t 2 |. ( 25
)
Theorem 2. Under assumptions 1, 3 and 5, if f m m (t) is a kernel estimator of the form (20) where the kernel K and the bandwidth σ := σ m are deterministic (i.e. do not depend on the data), such that σ

K < ∞, σ K 2 < ∞, R(K) < ∞ and if lim m→∞ σ m = 0, lim m→∞ mb m σ m = +∞, ( 26
)
then lim m→∞ E ( f m m (t) (z) -f t (z)) 2 = 0. ( 27
)
As an example, according to (13) from assumption 3, theorem 2 applies to the case of a Gaussian kernel and a bandwidth σ m = 1/ √ mb m . Unfortunately, it does not apply to the marginal density estimator presented in the next section, whose kernel is random.

Possible Choice of Density Estimator: the Self-Consistent Estimator

In the previous section we presented a general estimator of some discrepancy used to quantify how close a certain curve is to a set of other curves, called the Mean Marginal Likelihood. As explained, such plug-in estimator is based on the aggregation of other consistent density estimators trained on uniform bins. One may wonder what local density estimator to use in this situation.

As most statistical learning problems, density estimation can be tackled in a parametric or a nonparametric setting. In the first case, a specific class of density functions has to be fixed a priori, a finite set of unknown parameters needing to be tuned using information contained in the data. Such approaches, as for example Maximum Likelihood estimation, are known to be fast to train and evaluate, they have the best learning rate attainable and are usually very scalable and accurate if the model assumptions are correct. However, the nonparametric density learning techniques make little to no assumptions on the shape of the density to be estimated. As explained in section 2.1, the marginal densities at different times may greatly vary in shape, which is why we preferred to consider nonparametric estimators in this article and more precisely the popular kernel density estimators.

For d-dimensional data, kernel density estimators (KDE) have the following general form when trained on N i.i.d observations {x k } N k=1 of a random variable X of density f :

f N (x) = 1 N det H N k=1 K H -1 (x -x k ) . (28)
In ( 28), the function K : R d → R + is called a smoothing kernel and H ∈ GL d (R) is the kernel's bandwidth matrix. In general, the main drawback of kernel density estimators (28) lies on the subjective choice of the kernel K and bandwidth H. However, it is well-known folklore in the density estimation literature (see e.g. Wasserman, 2004, chapter 20.3) that KDE's accuracy is not really sensitive to the choice of K and depends mainly on the bandwidth H used. Several rules and heuristics have been suggested since then to choose such a parameter, but they are usually based on quite strict assumptions, such as Silverman's rule of thumb for the estimation of a 1-dimensional Gaussian density [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF]. Another possibility is to select H through cross-validation [START_REF] Stone | An asymptotically optimal window selection rule for kernel density estimates[END_REF] but this approach is computationally intensive, specially if d is larger than 1. For these reasons, we decided to consider a similar method proposed by [START_REF] Bernacchia | Self-consistent method for density estimation[END_REF], called the self-consistent density estimator. It consists indeed in a KDE whose kernel incorporates the bandwidth and is learned directly from the data, hence not requiring any parameter tuning. Its derivation is based on the use of a fixed-point equation to approximate the optimal kernel estimator in the sense of the Mean Integrated Squares-Error (MISE). The obtained estimator takes the form of the Fourier transform of the following characteristic function estimator:

Φsc (s) := N ∆(s) 2(N -1) 1 + 1 - (∆ min N ) 2 |∆(s)| 2 1 A N (s),
where s is the Fourier variable, N is the number of training observations {x k } N k=1 , ∆ is the empirical characteristic function

∆(s) := 1 N N k=1 e ix k •s , ∆ min N := 4(N -1) N 2 , (29)
i = √ -1 and 1 A N denotes the indicator function over an arbitrary subset A N ⊂ S N of the frequencies in

S N := s : |∆(s)| 2 ≥ (∆ min N ) 2 . ( 30
)
Bernacchia and Pigolotti (2011) proved for 1D data that, under mild assumptions on A N , the self-consistent estimator converges almost-surely to the true density f as the number of observations N grows and is hence (strongly) consistent. This result is summarized in the following theorem:

Theorem 3 ( [START_REF] Bernacchia | Self-consistent method for density estimation[END_REF].

Let [-t * ; t * ] = A N ⊂ S N be an interval of frequencies in R. Assuming that f is L 2 (E, R + ) and its Fourier transform Φ = F[f ] is L 1 (R, R + ), if the bounds of A N are such that lim N →∞ t * = ∞, lim N →∞ t * √ N = 0, ( 31
)
then the density estimator defined by

f N sc (x) := F -1 [ Φsc ](x), ∀x ∈ E (32)
converges almost surely to f as N tends to infinity: 2016)are shown to be several times faster to compute than a regular KDE. This is achieved thanks to the smart use of the Non Uniform Fast Fourier Transform [START_REF] Greengard | Accelerating the Nonuniform Fast Fourier Transform[END_REF] to compute the empirical characteristic function ∆. This property is particularly important in our case because of the potentially large number of trainings needed to compute the EMML, which is equal to the number of bins of T's partition.

P lim N →∞ f N sc (x) = f (x) = 1, ∀x ∈ E. ( 33
For all these reasons, all EMML numerical experiments presented in section 4 make use of this density estimator. More details concerning its derivation and implementation can be found in the supplementary material.

Application to the Assessment of Optimized Aircraft Trajectories

Experiments Motivation

In this section we illustrate our approach on real data recorded from m = 424 flights of the same medium haul aircraft T = {z 1 , . . . , z m }, which corresponds to 334 531 observation points. These trajectories are used to estimate the differential system describing the aircraft dynamics. Then, by numerically solving an optimal control problem defined using the estimated aircraft dynamics, a new trajectory y is obtained, supposed to minimize the overall fuel consumption for some future flight (see e.g. [START_REF] Rommel | Aircraft dynamics identification for optimal control[END_REF]. Note that:

1. The dynamics model is not guaranteed to be valid outside of the region occupied by the data used to estimate it. Hence, it is natural to want the optimized trajectory to avoid going too far from its validity region. 2. Furthermore, it is desirable for the proposed trajectory not to be too unusual compared to standard climb profiles for better acceptance by the pilots and Air Traffic Control. The two previous points motivate the need for an indicator of closeness between the optimized trajectory and the set of recorded flights.

Experiments Design

Training set. The training data used for our experiments were extracted from the Quick Access Recorder (QAR) of the same aircraft, whose sampling rate is of one measurement per second. We only used the recordings of 5 variables, which are the altitude h, the true airspeed V , the path angle γ, the angle of attack α and the throttling position N 1 . These variables are not all directly accessible and some were computed from other measurements using standard formulas from flight mechanics (see e.g. [START_REF] Rommel | Aircraft dynamics identification for optimal control[END_REF]. Only the portion corresponding to the climb phase of these signals was kept for our experiments, i.e. data corresponding to altitudes between 1524 m = 5000 ft and the top of climb (cruise altitude, specific to each flight). The training set of curves obtained by the described procedure is displayed on figure 2a.

Test set. In order to evaluate the estimated mean marginal likelihood on relevant examples, the following test flights were considered:

1. 50 real flights, extracted from the training set before training; 2. 50 simulated trajectories which were optimized as described in section 4.1, with constraints keeping the resulting speed V and N 1 between reasonable operational bounds; 3. and another 50 simulated trajectories optimized without the operational constraints. We evaluated the likelihood of these 150 test trajectories using the EMML and the competing methods described in the following section in order to assess and compare their discriminative power and computation efficiency.

Alternate Approaches Based on Standard Methods

The problem of quantifying how close a newly observed random curve is with respect to a set of observations from the same stochastic process has not been treated by the statistical learning literature to our knowledge. However, this problem is related to other more standard approaches from Functional Data Analysis and Conditional Density Estimation, which could be adapted quite straightforwardly for this purpose. For this reason, we discuss in the following paragraphs the characteristics of two of these other existing methods, before comparing them to our approach in the numerical results of section 4.5.

Functional Principal Components Analysis

Functional Principal Components Analysis (FPCA) is a standard tool in FDA capable of building a small number of descriptors which summarize the structure of a set of random functions. As explained for example in [START_REF] Nicol | Functional principal component analysis of aircraft trajectories[END_REF], this dimension reduction method can be used to project the train set of infinite-dimensional random trajectories into a finite (low) dimensional space.

Following the same reasoning used to derive the MML, our idea here consists in estimating the density function of these low dimensional representations of the training set. Then, after projecting the new trajectory y into the same descriptors, we can evaluate the density estimate at it and obtain an approximation of its likelihood.

Least-Squares Conditional Density Estimation

From a completely different point of view, we could forget for a moment that we are considering a random process Z and look at (T, Z T ) as a pair of standard random variables valued on the finite dimensional space T × E. In this case, we could see the marginal densities f t as the conditional probability density functions

f Z T |T (t, z) = f (T,Z T ) (t, z) f T (t) , (t, z) ∈ T × E. ( 34
)
We could hence estimate (34) at the observed points of the new trajectory y and use them to compute the EMML indicator (8). It is however well-known in density ratio estimation that approximating f (T,Z T ) (t, z) and f T (t) separately before building the ratio in ( 34) is not a good idea because it magnifies the errors. For this reason, [START_REF] Sugiyama | Conditional density estimation via least-squares density ratio estimation[END_REF] proposed to use a linear model for this purpose

f Z T |T (t, z) = θ φ(t, z), (35)
where θ = (θ 1 , . . . , θ p ) is a vector of scalar parameters and φ(t, z) = (φ 1 (t, z), . . . , φ p (t, z)) is a family of nonnegative basis functions. The parameters θ are then chosen so as to minimize a L 2 -penalized least-squares criterion, which is shown to have a closed-form solution. This method was coined Least-Squares Conditional Density Estimation (LS-CDE) by the authors, and is also known as Unconstrained Least-Squares Importance Fitting (uLSIF) in the density ratio estimation literature [START_REF] Kanamori | A least-squares approach to direct importance estimation[END_REF]. The extensive numerical results presented in [START_REF] Sugiyama | Conditional density estimation via least-squares density ratio estimation[END_REF] indicate that this approach have state-of-the-art accuracy in conditional density estimation.

Algorithms Settings

For all the methods tested, the altitude h played the role of "time". This is a natural assumption made when optimizing the climb profile of a civil airliner, since the altitude is an increasing function of the time and every other variable depends on it. This allowed us to reduce the dimension of our problem from 5 to 4. MML with Self-Consistent Kernel Estimator settings. The python library fastKDE (O' Brien et al., 2014;[START_REF] O'brien | A fast and objective multidimensional kernel density estimation method: fastkde[END_REF] was used to compute the marginal densities from the bins data. It contains the implementation of the Self-Consistent kernel estimator described in section 3. The precision of the density estimations were set to single. The confidence levels were approximated by numerical integration using the trapezoidal rule over a fine grid of approximately 300 points per bin.

Concerning bin sizes, we chose to use an uneven partition in our experiments. The reason for this are the climb-steps visible in the trajectories between 3000 and 4000 m, which correspond to phases during which the aircraft decreases considerably its ascent speed, leading to slowly increasing altitudes. Such behaviors translate into rapidly increasing speeds V with respect to the altitude, as well as into plummeting values of γ and N 1 (see figure 2a). This brought us to consider tighter bins around these climb-step altitudes:

• between 1524 and 3000m and between 4000 and 12000m, we partitioned the altitudes homogeneously into bins of size b

(1) m = 21m 1/ √ m (which satisfies assumption 3); • between 3000 and 4000m, we used a bin size twice smaller b

(2) m = 10m b (1)
m /2; FPCA settings. Concerning the Functional Principal Components Analysis method, all training and testing flights were resampled on an equispaced grid of altitudes, using a step size of 5m. The trajectories were then centered and decomposed into a basis of 128 cubic B-splines. The SVD decomposition was carried using the PCA class from scikit-learn python library [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. We kept 4 components for each variable (V, γ, α and N 1 ), which was enough to explain more than 90%, 65%, 60% and 75% of their respective variance. A Gaussian mixture model was used to estimate the density of the training trajectory scores obtained by the projection into the principal functions. The model was trained using a standard EM algorithm, implemented in scikit-learn as well. The number of components was selected between 1 and 5 using the Bayesian information criterion (BIC).

LS-CDE settings. For the Least-Squares Conditional Density Estimation, the python package densratio [START_REF] Makiyama | densratio, A Python Package for Density Ratio Estimation[END_REF], implementing the uLSIF method from [START_REF] Kanamori | A least-squares approach to direct importance estimation[END_REF] was used and adapted. The basis functions chosen were p = 100 Gaussian kernels with the same variance σ and different centers. These centers were randomly drawn from the training data points using a uniform distribution, as suggested in [START_REF] Sugiyama | Conditional density estimation via least-squares density ratio estimation[END_REF]. The variance σ, as well as the L 2 penalty weight λ needed for minimizing the least-squares criterion were selected by cross-validation.

Results and Comments

Figures 2b and2c show heatmaps encoding the estimated marginal likelihoods using the normalized density (3) and the confidence level (4). We notice that both figures are similar and seem to catch the shape of the plot from figure 2a, including the multi-modalities visible for N 1 below h = 4000m for example.

Table 1a contains the estimated Mean Marginal Likelihood scores averaged over each test flight category. The training was carried on each dimension separately and the average total training time was of 5 seconds on a laptop (2.30 GHz, 7.7 GB). First of all, we notice for both types of scaling functions that the test flight categories are nicely separated by three really distinct ranges of scores. Furthermore, the higher differences between the two types of optimized flights can be seen for the variables V and N 1 , which makes sense since those are the variables which are left free for Opt2 flights and constrained for Opt1 flights. As expected from figures 2b and 2c, the performances of both confidence level and normalized density based MML are comparable in terms of discrimination power and seem adequate for the task of assessing optimized aircraft climb trajectories.

Table 1b contains the estimated Mean Marginal Likelihood scores in a 2-dimensional setting, where the pairs (V, γ) and (α, N 1 ) have been treated together. The average training time needed here was 16 times larger than in the 1D case, i.e. 1 minute 20 seconds. The scores observed are globally really low and the test flight categories are not well separated. Moreover, we expected to obtain large scores for the real flights, since we used the marginal densities of their category to build the criterion, but this is not the case here. We conclude that the MML criteria based on the self-consistent estimator does not work so well in higher dimension and we suspect this to be related to the curse of dimensionality. Indeed, it is well-known (see e.g. Wasserman, 2004, chapter 21.3) that as the dimension grows, the amount of data needed to attain a given accuracy with kernel density estimators skyrockets, which may explain this poor performance.

From a practical point of view, a reference value or threshold is needed if one wanted to use our method to determine automatically whether a given optimized flight should be accepted or not. In such a context, a quite straightforward solution would be to use a leave-one-out cross-validation approach: compute the MML score of each real flight leaving it out of the training data and then averaging over the obtained scores. These reference values have been computed for our dataset and are summarized in table 2. We note that the values obtained are very close to the average scores of the real flights showed in table 1a. 12 0.09 ± 0.05 0.11 ± 0.05 0.03 ± 0.03 0.05 ± 0.03 0.06 ± 0.03 0.01 ± 0.01 (α, N1) 0.03 ± 0.02 0.01 ± 3e-3 0.01 ± 2e-3 0.02 ± 0.01 4e-3 ± 2e-3 3e-3 ± 1e-3 Mean 0.06 ± 0.03 0.06 ± 0.02 0.02 ± 0.01 0.03 ± 0.02 0.03 ± 0.02 0.01 ± 0.01 Table 3a contains the scores obtained using the Functional PCA based method presented in section 4.3.1. The training time needed here was of 20 seconds in average. As for the MML in 2D, the real flights' scores are surprisingly low and the two types of simulated trajectories are not well discriminated by the criterion. This might be caused by the fact that this method encodes each training trajectory by a single point in the 4-dimensional space spanned by the principal functions. The training set obtained is hence of m = 424 points, which might be too small too attain sufficient accuracy from the Gaussian mixture density estimator in such a high dimension. The principal functions used and scatter plots of the projected trajectories can be found in the online version.

Concerning the LS-CDE approach, because the algorithm needs large Gram matrices (of size O(nm 2 )) to be stored, we encountered several memory problems when trying to run it on our dataset of 334 531 observation points. For this reason, our results were obtained by applying it to 100 uniform batches. These batches were obtained by partitioning the data according to the altitude. Although the three categories are well-separated by this method, as shown on table 3b, the total time needed to train the estimators on every batch was approximately 14 hours.

We didn't test both alternate methods in the 2D setting since the problems observed in 1D (curse of dimensionality for FPCA and memory/time for LS-CDE) would be aggravated. In conclusion, our numerical results indicate that the MML criterion has better discriminative power than FPCA and LS-CDE for the task of assessing curves relatively to a set of "good" examples. Furthermore, the training time and memory needed for using LS-CDE in datasets of this size seems crippling. Concerning the FPCA method, it does not seem to be applicable to datasets with so few curves and present a higher training time than MML.

Conclusions

In this paper we proposed a new approach for a problem which seemed unaddressed by the statistical learning community: quantifying the closeness from a curve to a set of random functions. We introduced a class of probabilistic criteria for this context called the Mean Marginal Likelihood (MML), and analyzed two possible scaling functions used to build them. We also derived a class of estimators of our criteria, which make use of local density estimators proved to consistently approximate the marginal densities of a random process. For practical applications, we suggested a particular flexible density estimator believed to have the right properties needed in this setting, called the self-consistent kernel estimator.

Numerical experiments using real aircraft data were carried to compare the MML with other well-established approaches from Functional Data Analysis and Conditional Density Estimation. The results show that, although the MML does not take into account the temporal structure of the data as other standard functional data analysis methods, it is a good candidate for the type of applications suggested. This seems to be especially the case if the number of training trajectories is too small for using FPCA or if the total number of observation points is too large for conditional density estimation. Moreover, the training times obtained for MML are by far the shortest among the compared methods, which confirms the relevance of the self-consistent kernel estimator. Furthermore, the ease to visualize, localize and interpret discrepancy zones allowed by MML make it a good exploratory analysis tool for functional data (see e.g. figure 2). We also note that our method does not perform as well in a multidimensional setting, but that the training time should not be an obstacle. In future work we intend to test the MML with parametric density estimators, which should be less affected by the curse of dimensionality.

Appendices

A. Marginal densities consistency proof

In this section we prove theorem 1 from section 2.3. It relies on 5 lemmas stated and proved hereafter.

Lemma 1. Let t ∈ T and z ∈ E. Under assumptions 2 and 3,

f m m (t) (z) converges to f t (z): lim m→∞ |f t (z) -f m m (t) (z)| = 0. ( 36
)
Proof.

|f t (z) -f m m (t) (z)| = f t (z) - B m m (t) f s (z)ν(s)ds B m m (t) ν(s)ds , = B m m (t) (f t (z) -f s (z))ν(s)ds B m m (t) ν(s)ds , ≤ B m m (t) |f t (z) -f s (z)|ν(s)ds B m m (t)
ν(s)ds .

(37)

According to assumption 2,

|f t (z) -f s (z)| ≤ L|t -s| ≤ L|τ m m (t) -τ m m (t)-1 | = Lb m . ( 38
)
Hence,

|f t (z) -f m m (t) (z)| ≤ Lb m B m m (t) ν(s)ds B m m (t) ν(s)ds = Lb m . ( 39
)
Since b m → 0 by assumption 3, the conclusion follows.

Lemma 2. Let m ∈ N * and t ∈ T. Under assumption 3, the probability that T m m (t) (the subset of training points whose sampling time fall in the bin containing t) contains at most one observation point per curve is asymptotically equal to 1, meaning that for large enough m, the observations in T m m (t) will be independent with high probability:

lim m→∞ P(N m r, m (t) ≤ 1) = 1, r = 1, . . . , m, (40)
where N m r, m (t) denotes the number of observations of z r in T m m (t) .

Proof. Let 1 ≤ r ≤ m and

T m r, m (t) := {(t r k , z r k ) : t r k ∈ B m m (t) ; 1 ≤ k ≤ n} (41)
be the set of observations of the r th curve with times lying in the m (t) th bin B m m (t) . Let N m r, m (t) be the number of elements in T m m (t) . The random variable N m r, m (t) follows a binomial law B(n, P m m (t) ), where

P m m (t) = B m m (t) ν(t)dt (42)
is the probability of a new observation of the r th curve falling in T m r, m (t) . The probability of N m r, m (t) being at most equal to 1 writes

P(N m r, m (t) ≤ 1) = 1 k=0 n k (P m m (t) ) k (1 -P m m (t) ) n-k , = (1 -P m m (t) ) n + nP m m (t) (1 -P m m (t) ) n-1 , = (1 -P m m (t) ) n-1 (1 + (n -1)P m m (t) ). ( 43
) As P m m (t) = B m m (t) ν(t)dt ≤ b m ν + , (44)
we have

P(N m r, m (t) ≤ 1) ≥ (1 -P m m (t) ) n-1 ≥ (1 -b m ν + ) n-1 . (45)
Since b m → 0 according to assumption 3, we obtain

lim m→∞ P(N m r, m (t) ≤ 1) = 1, r = 1, . . . , m. (46)
Lemma 3. Under assumption 3 and for any t ∈ T, the number N m m (t) of observations falling in T m m (t) diverges in probability to +∞:

∀M > 0, lim m→∞ P N m m (t) > M = 1. ( 47
)
Proof. Let M > 0. As in the proof of lemma 2, we have N m m (t) ∼ B(mn, P m m (t) ) and hence,

P N m m (t) ≤ M = M k=1 nm k (P m m (t) ) k (1 -P m m (t) ) nm-k . (48)
As the sum has is finite when we fix M , the limit of expression (48) when m tends to infinity is the sum of the limits of its terms.

nm k (P m m (t) ) k (1 -P m m (t) ) nm-k ≤ nm k (ν + b m ) k (1 -ν -b m ) nm-k , ∼ m→∞ (ν + nmb m ) k (1 -ν -b m ) nm , = (nmb m ) k exp [nm log (1 -ν -b m )] , ∼ m→∞ (nmb m ) k exp [-ν -nmb m + o(mb m )] . (49)
Using ( 13) from assumption 3, we obtain that (49) tends to 0 which proves that

lim m→∞ P N m m (t) > M = 1 -lim m→∞ P N m m (t) ≤ M = 1.
In what follows, for any M > 0 we denote C M the following event:

C M := {N m m (t) > M } m r=1 {N m r, m (t) ≤ 1}. ( 50
)
Lemma 4. For any t ∈ T, if assumption 3 holds,

∀M > 0, lim m→∞ P (C M ) = 1. ( 51
)
Proof. Let M > 0. We have by definition of the conditional probability

P (C M ) = P N m m (t) N m r, m (t) ≤ 1; r = 1, . . . , m × P N m r, m (t) ≤ 1; r = 1, . . . , m . (52)
As the variables N m r, m (t) are independent

P N m r, m (t) ≤ 1; r = 1, . . . , m = m r=1 P N m r, m (t) ≤ 1 (53)
which tends to 1 when m grows according to lemma 2. Furthermore, as

N m m (t) = m r=1 N m r, m (t) ∼ B(nm, P m m (t)
), we know that the random variable (N m m (t) |N m r, m (t) ≤ 1; r = 1, . . . , m) is also a binomial random variable where at least m(n-1) Bernouilli have failed, leaving m trials. Hence, we know that (N m m (t) |N m r, m (t) ≤ 1; r = 1, . . . , m) ∼ B(m As in the proof of lemma 3,

P N m m (t) ≤ M N m r, m (t) ≤ 1; r = 1, . . . , m = M k=1 m k (P m m (t) ) k (1 -P m m (t) ) m-k , (54)
whose limit is the sum of the limits of the sum terms. As in the proof of lemma 3, we have m k

(P m m (t) ) k (1 -P m m (t) ) nm-k (mb m ) k exp [-ν -mb m + o(mb m )] , (55)
and hence, by using (13) from assumption 3,

lim m→∞ P N m m (t) ≤ M N m r, m (t) ≤ 1; r = 1, . . . , m = 0. (56)
Combining ( 52), ( 53) and ( 56) we obtain (51).

Lemma 5. For any z ∈ E and t ∈ T, if assumptions 3 and 4 hold for some function Θ used to compute f m m (t) , then

∀ε > 0, lim m→∞ P | f m m (t) (z) -f m m (t) (z)| < ε = 1. ( 57
)
Proof. The randomness of f m m (t) (z) comes from both the sample of observations drawn from the bin's density f m m (t) and the random number of observations falling in the bin N m m (t) . Hence, the idea here is to separate these two sources of randomness by conditioning on one of them. As we would like to use the result from assumption 4, it makes sense to condition on N m m (t) here. Indeed, for any m > 0, f m

m (t) (z) = Θ[T m m (t) ](z), where T m m (t) is a sample of N m m (t)
observations drawn from f m m (t) . Hence, according to assumption 4, for any ε > 0, for any α 1 > 0, there is some N ε,α1,m such that

P | f m m (t) (z) -f m m (t) (z)| < ε C Nε,α 1 ,m > 1 -α 1 . (58)
As seen in lemma 2, the conditioning on N m r, m (t) ≤ 1 comes from the fact that when T m m (t) contains at most one observation per curve, those observations are independent. Furthermore, for ε, α 1 fixed, let M := sup m {N ε,α1,m } (which exists according to assumption 4). In this case, when assumption 3 holds, we have from lemma 4 that for any α 2 > 0,there is some m ε,α1,α2 such that,

m ≥ m ε,α1,α2 ⇒ P (C M ) > 1 -α 2 . (59)
Hence, as

N m m (t) > M ⇒ N m m (t) > N ε,α1,m , we have m ≥ m ε,α1,α2 ⇒ P C Nε,α 1 ,m > 1 -α 2 . (60)
Moreover, we can also write the following useful inequality:

P | f m m (t) (z) -f m m (t) (z)| < ε > P {| f m m (t) (z) -f m m (t) (z)| < ε} C Nε,α 1 ,m = P | f m m (t) (z) -f m m (t) (z)| < ε C Nε,α 1 ,m × P C Nε,α 1 ,m . (61)
By combining (58), ( 60) and ( 61) we get that for any ε, α 1 , α 2 > 0, there is

m ε,α1,α2 such that m ≥ m ε,α1,α2 ⇒ P | f m m (t) (z) -f m m (t) (z)| < ε > 1 -α 3 , (62)
with

α 3 = α 1 + α 2 -α 1 α 2 .
Using these lemmas we can finally prove theorem 1:

Proof of Theorem 1. Let (t, z) ∈ T × E and let ε > 0. According to lemma 1, under assumptions 2 and 3, there is

m ε ∈ N * such that for any m ≥ m ε , |f m m (t) (z) -f t (z)| < ε (63) As | f m m (t) (z) -f t (z)| ≤ | f m m (t) (z) -f m m (t) (z)| + |f m m (t) (z) -f t (z)|, ≤ | f m m (t) (z) -f m m (t) (z)| + ε, (64)
we have

{| f m m (t) (z) -f t (z)| > 2ε} ⊂ {| f m m (t) (z) -f m m (t) (z)| > ε}, and P | f m m (t) (z) -f t (z)| > 2ε ≤ P | f m m (t) (z) -f m m (t) (z)| > ε . (65)
According to lemma 5 we have that

lim m→∞ P | f m m (t) (z) -f m m (t) (z)| > ε = 0, ( 66
)
which allow us to obtain ( 17) from ( 65).

In the following section we prove an even stronger convergence (in the L 2 sense) for standard kernel density estimators.

B. L 2 consistency for kernel density estimators

We assume here that f m m (t) is a standard kernel density estimator of the form (20).

Remark 2. Such an estimator is only defined for N m m (t) > 0. This is why, in the remaining section, we consider the conditioned random variables (N m m (t) |N m m (t) > 0) and ( f m m (t) |N m m (t) > 0), which are still denoted N m m (t) and f m m (t) for simplicity. Furthermore,

P •|N m m (t) > 0 = P • {N m m (t) > 0} P N m m (t) > 0 , (67)
and, according to lemma 3, lim m→∞ P N m m (t) > 0 = 1. In the following, we use the following notations for the conditional expectation and variance:

Ẽ

[•] := E •|N m m (t) > 0 , Ṽar [•] := Var •|N m m (t) > 0 . (68)
In the remaining of this section we derive the conditions under which f m m (t) converges in expected squared-error to f t . We recall that a sufficient condition for this is having its bias and variance tending to 0. The proof was greatly inspired by the derivations presented in [START_REF] Scott | Multivariate density estimation: theory, practice, and visualization[END_REF] for the standard multivariate case.

Similarly to lemma 1, we can prove the following convergence result: Lemma 6. Under assumption 5, for any

(t, z) ∈ T × E, d 2 f m m (t) dz 2 (z) := f m m (t) (z) = τ m τ m -1 f t (z)ν(t)dt τ m τ m -1 ν(t)dt , (69)
and

lim m→∞ f m m (t) (z) -f t (z) = 0. ( 70
)
Proof. Similar argument to lemma 1.

Furthermore, the following bounds of the estimator's bias and variance hold under the same assumption:

Lemma 7. If assumption 5 holds, Ẽ f m m (t) (z) -f m m (t) (z) ≤ 1 2 f m m (t) ∞ Ẽ σ 2 Kσ , (71)
Ṽar f m m (t) (z) ≤ f m m (t) (z) Ẽ R(K σ ) N m m (t) + 1 2 f m m (t) ∞ Ẽ σ 2 K 2 σ N m m (t) + 1 4 f m m (t) 2 ∞ Ẽ (σ 2 Kσ ) 2 . ( 72
)
Proof. From the law of total expectation (e.g. Wasserman, 2004, p.55)

Ẽ f m m (t) (z) = Ẽ Ẽ f m m (t) (z)|N m m (t) . ( 73
)
For large enough m, f m m (t) writes as an empirical average of independent identically distributed random variables Z m,t of density f m m (t) and hence

Ẽ f m m (t) (z)|N m m (t) = Ẽ K σ (z -Z m,t )|N m m (t) = Ẽ [K σ (z -Z m,t )] , = K σ (z -x)f m m (t) (x)dx = K σ (w)f m m (t) (z -w)dw. (74)
The second order integral Taylor expansion of f m m (t) around z gives

f m m (t) (z -w) = f m m (t) (z) -f m m (t) (z)w + 1 0 (1 -x) f m m (t) (z -xw)w 2 dx, (75) which leads to Ẽ f m m (t) (z)|N m m (t) = f m m (t) (z) K σ (w)dw -f m m (t) (z) wK σ (w)dw + 1 0 (1 -x) f m m (t) (z -xw)w 2 K σ (w)dx dw. (76)
As K σ is a symmetric probability density function, we have

K σ (w)dw = 1, wK σ (w)dw = 0, ( 77
) which leads to Ẽ f m m (t) (z)|N m m (t) = f m m (t) (z) + 1 0 (1 -x) f m m (t) (z -xw)w 2 K σ (w)dx dw. ( 78
)
By combining ( 73) and ( 78) we obtain the following bias expression

Ẽ f m m (t) (z) -f m m (t) (z) = Ẽ 1 0 (1 -x) f m m (t) (z -xw)w 2 K σ (w)dx dw , (79)
proving bound (71).

Concerning the variance, we apply the law of total variance (e.g. Wasserman, 2004, p.55):

Ṽar f m m (t) (z) = Ẽ Ṽar f m m (t) (z)|N m m (t) + Ṽar Ẽ f m m (t) (z)|N m m (t) . ( 80
)
As in (74), we may express the conditional variance of f m m (t) (z) using the conditional variance of the kernel, which brings the first term in (80) to

Ẽ Ṽar f m m (t) (z)|N m m (t) = Ẽ 1 N m m (t) Ṽar K σ (z -Z m,t )|N m m (t)
.

(81)

Furthermore, the kernel variance can be developped as follows

Ṽar K σ (z -Z m,t )|N m m (t) = Ẽ K σ (z -Z m,t ) 2 |N m m (t) -Ẽ K σ (z -Z m,t )|N m m (t) 2 , (82)
which is hence smaller than the first term, i.e. the kernels second moment. Using integral Taylor expansion of f m m (t) around z (75) truncated to the first order, such a quantity can be written as follows:

Ẽ K σ (z -Z m,t ) 2 |N m m (t) = K σ (z -λ) 2 f m m (t) (λ)dλ = K σ (w) 2 f m m (t) (z -w)dw = f m m (t) (z)R(K σ ) -f m m (t) (z) wK σ (w) 2 dw =0 + 1 0 (1 -x)w 2 f m m (t) (z -xw)K σ (w) 2 dx dw, (83)
where we used the fact that K 2 σ is an even function. By using ( 81)-( 82) and ( 83) we get

Ẽ Ṽar f m m (t) (z)|N m m (t) ≤ f m m (t) (z) Ẽ R(K σ ) N m m (t) + 1 2 f m m (t) ∞ Ẽ σ 2 K 2 σ N m m (t)
.

We still need to bound the second term in (80). We have

Ṽar Ẽ f m m (t) (z)|N m m (t) = E Ẽ f m m (t) (z)|N m m (t) -Ẽ Ẽ f m m (t) (z)|N m m (t) Ẽ f m m (t) (z) 2 . ( 85
)
By denoting

G m = 1 0 (1 -x) f m m (t) (z -xw)w 2 K σ (w)dx dw (86)
and plugging ( 78)-( 79) in ( 85), we obtain

Ṽar Ẽ f m m (t) (z)|N m m (t) = Ẽ G m -Ẽ [G m ] 2 ≤ Ẽ G 2 m . (87) Furthermore, Ẽ G 2 m = Ẽ (σ 2 K σ 2 m ) 2 . ( 91
)
Proof. Obtained directly by using ( 21)-( 23) together with lemma 7.

We know from lemmas 1 and 6 that, under assumptions 3 and 5:

lim m→∞ f m m (t) (z) = f t (z) < ∞, ( 92
) lim m→∞ f m m (t) (z) = f t (z) < ∞. (93)
Hence, by looking at expressions (90) and (91), it seems clear that the convergence of the bias and the variance to zero will strongly depend on the asymptotics of Ẽ 1/N m m (t) . This motivates the following lemma:

Lemma 9. If assumption 1 to 3 hold, as m → ∞ Ẽ 1 N m m (t) = O 1 nmP m m (t) = O 1 mb m . (94)
Proof. We recall that N m m (t) ∼ B(nm, P m m (t) ), with

P m m (t) = B m m (t)
ν(t)dt. According to the main theorem proved in [START_REF] Cribari-Neto | A note on inverse moments of binomial variates[END_REF],

S 1 := E 1 1 + N m m (t) = nm k=0 nm k (P m m (t) ) k (1 -P m m (t) ) nm-k 1 + k = O 1 nmP m m (t) . ( 95
)
As noted in remark 2, f m m (t) is not defined for N m m (t) = 0. This motivates the conditioning of the variable 1

N m m (t)
by the event N m m (t) > 0 and we have

S 0 := E 1 N m m (t) N m m (t) > 0 = nm k=1 nm k (P m m (t) ) k (1 -P m m (t) ) nm-k k . (96)
By computing the difference between the S 1 and S 0 we obtain the following bound:

S 0 -S 1 = - nm 0 (1 -P m m (t) ) nm + nm k=1 nm k (P m m (t) ) k (1 -P m m (t) ) nm-k k(k + 1) ≤ nm k=1 nm k (P m m (t) ) k (1 -P m m (t) ) nm-k k + 1 ≤ S 1 . (97)
We conclude that S 0 ≤ 2S 1 and hence

E 1 N m m (t) N m m (t) > 0 = O 1 nmP m m (t) . (98)
Finally, the last equality in (94) comes from the fact that P m m (t) ≥ ν -b m .

In conclusion, we can now prove theorem 2 stating that the kernel marginal density estimator will be consistent in expected squared-error (which implies convergence in probability stated in theorem (1)).

Proof of Theorem 2. By ( 90) and ( 26), the bias of f m m (t) (z) converges to 0:

lim m→∞ Ẽ f m m (t) (z) -f m m (t) (z) = 0. ( 99
)
Similarly, as σ m converges to 0, the two last terms in (91) shrink. Concerning the first term, we can conclude from lemma 9 and from condition (26) that

lim m→∞ Ẽ 1 N m m (t) 1 σ m = 0, (100)
which means that the variance of f m m (t) (z) also converges to 0:

lim m→∞ Ṽar f m m (t) (z) = 0. ( 101
)
From the bias-variance decomposition of the expected squared-error between f m m (t) (z) and f m m (t) (z) we have

Ẽ ( f m m (t) (z) -f m m (t) (z)) 2 = Ẽ f m m (t) (z) -Ẽ f m m (t) (z) + Ẽ f m m (t) (z) -f m m (t) (z) 2 , = Ẽ f m m (t) (z) -Ẽ f m m (t) (z) 2 + Ẽ Ẽ f m m (t) (z) -f m m (t) (z) 2 .
Furthermore, by Jensen inequality we have that,

1 2 Ẽ ( f m m (t) (z) -f t (z)) 2 ≤ Ẽ ( f m m (t) (z) -f m m (t) (z)) 2 + (f m m (t) (z) -f t (z)) 2 . (102)
Finally, by using lemma 1 in conjunction with (99), ( 101) and ( 102), we obtain that both terms in (102) tend to 0, leading to result (27).

C. Derivation of the Self-Consistent Density Estimator

In this section we present in more details the self-consistent estimator proposed by [START_REF] Bernacchia | Self-consistent method for density estimation[END_REF]

C.0.1. Optimal Kernel Density Estimator

Let S N = {z k } N k=1 ⊂ E be a sample of N observations drawn from a common density function f . We suppose that f ∈ L 2 (E, R + ) and that N > 0 is deterministic. We consider a kernel estimator of f

f (z) := 1 N N k=1 K(z -z k ), ∀z ∈ E, (103)
where

K : R d → R is a smoothing kernel in L 2 (R d , R).
One can interpret (103) as a kernel density estimator with an implicit bandwidth H ∈ GL d (R + ) hidden inside the expression of the kernel function K

(z) = K(H -1 z)/ det H, ∀z ∈ R d .
Denote by E the expectation relative to the random sample S N , defined for any deterministic function ϕ :

E N → R by E ψ(S N ) := E • • • E ϕ(z 1 , . . . , z N )f (z 1 ) . . . f (z N )dz 1 . . . dz N . (104)
In this context, a common quality measure of density estimators is the Mean Integrated Squared Error:

M ISE := E E ( f (z) -f (z)) 2 dz . (105)
As we will show, it becomes relatively easy to minimize such a criterion with regard to the choice of the kernel K once we've shifted it to the Fourier domain. Hence, for any function v ∈ L 2 (R d , R), we define its Fourier transform hereafter with the following convention

F[v](s) := R d v(z)e iz•s dz, ∀s ∈ R d , (106)
where i = √ -1, its inverse being defined by

F -1 [v](z) := 1 2π R d v(s)e -iz•s ds, ∀z ∈ R d . (107)
As f, f ∈ L 2 , Plancherel's theorem gives that

M ISE = 1 2π E R d | Φ(s) -Φ(s)| 2 ds , (108) where
The function ∆ is commonly called the empirical characteristic function (ECF). Plugging (110) in the MISE expression (108) and expanding the square gives:

M ISE = 1 2π R d |κ| 2 E |∆| 2 + |Φ| 2 -κ(E [∆] Φ * + E [∆ * ] Φ) ds, (113)
the arguments s being omitted for lighter notation and c * denoting the complex conjugate of any c ∈ C. Furthermore, as shown in [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF], section 1.3, lemma 1.2) for example, the ECF is an unbiased estimator of the characteristic function

E [∆(s)] = 1 N N k=1 E e is•z k = 1 N N k=1 +∞ -∞ e is•z k f (z k )dz k = Φ(s), (114)
and its second moment is

E |∆(s)| 2 = E [∆(s)∆(s) * ] = E [∆(s)∆(-s)] (115) = E   1 N 2 j,k:j =k e is•(zj -z k )   + 1 N (116) = 1 N 2 j,k:j =k E e is•zj e -is•z k + 1 N (117) = 1 N 2 j,k:j =k E e is•zj E e -is•z k + 1 N (118) = 1 N 2 j,k:j =k Φ(s)Φ(-s) + 1 N (119) = N -1 N Φ(s)Φ(-s) + 1 N (120) = N -1 N |Φ(s)| 2 + 1 N . ( 121
)
Passing from line (117) to ( 118) is based on the assumption that the random variables {z k } N k=1 are independent. Replacing ( 114) and ( 121) in (113), we obtain

M ISE = 1 2π R d |κ| 2 N (1 -|Φ| 2 ) + |Φ| 2 (1 -κ) 2 ds. (122)
As initially shown in [START_REF] Watson | On the estimation of the probability density[END_REF], expression (122) can be minimized with respect to the transformed kernel κ. Indeed, for any s ∈ R d , we get the following first order optimality condition by differentiating the quadratic integrand in (122) relative to κ(s): One practical problem with the self-consistent estimator is that f N sc = F -1 [ ΦN sc ] is not lower-bounded by zero. This can be corrected by translating f N sc downwards until the positive part integrates to one and then setting the negative part to 0. Indeed, it was proven by [START_REF] Glad | Correction of density estimators that are not densities[END_REF] that such a transformation induces no cost in terms of MISE accuracy.

1 N κ(1 -|Φ| 2 ) -|Φ| 2 (1 -κ) = 0, (123)
Another practical drawback with estimator ΦN sc is that the direct computation of the empirical characteristic function ∆ can be expensive: O(N • M ) exponential evaluations, where M is the number of frequency points s ∈ R d . Noting from definition (112) that the expression of ∆ is equivalent to some Discrete Fourier Transform

∆(s) = 1 N N k=1 a k e is•z k , ( 135
)
where the Fourier coefficients a k are all equal to 1, the idea of using the Fast Fourier Transform algorithm (FFT) proposed by [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF] seems natural. However, the latter only applies to the case of uniformly spaced data, which is not the case of {z k } N k=1 . For this reason, [START_REF] O'brien | Reducing the Computational Cost of the ECF Using a nuFFT: A Fast and Objective Probability Density Estimation Method[END_REF] proposed to use an implementation of Nonuniform Fast Fourier Transform (NUFFT) developped by [START_REF] Greengard | Accelerating the Nonuniform Fast Fourier Transform[END_REF].

It consists in interpolating the original data points {z k } N k=1 on a new equispaced grid {z ). The analysis conducted in [START_REF] Greengard | Accelerating the Nonuniform Fast Fourier Transform[END_REF] indicates that simple precision can be achieved in (137) by normalizing the data {z k } N k=1 to the range [0, 2π] and setting 2N c = 12, Ñ = 2N and σ 2 = 24/N 2 . Hence, for a desired precision, this step of the algorithm introduces no additional hyperparameter to be tuned (see [START_REF] Greengard | Accelerating the Nonuniform Fast Fourier Transform[END_REF] for the double-precision settings).

Figure 1 :

 1 Figure 1: Illustration of the confidence level in the case of a univariate bimodal distribution

)

  It has been demonstrated through extensive numerical experiments in Bernacchia and Pigolotti (2011); O'Brien et al. (2014); O'Brien et al. (2016) that the self-consistent estimator achieves state-of-the-art MISE accuracy for many types of underlying densities. Furthermore, modern implementations of this estimator proposed by O'Brien et al. (2014); O'Brien et al. (

Figure 2 :

 2 Figure 2: Estimated marginal densities using two types of scaling functions

  and extended by O'Brien et al. (2014); O'Brien et al. (2016).

e

  Φ := F[f ], usually called the characteristic function, and Φ := F[ f ]. By noticing that f can be seen as the convolution between the kernel and a sum of Dirac functions centered on iz k •s ∈ C.

  leading to the optimal transformed kernelκ opt (s) := N N -1 + |Φ(s)| -2 , ∀s ∈ R d . (124)Hence the optimal density estimator relative to the MISE is given by fopt(z) = F -1 [ Φopt ](z),(125)C.0.3. Practical Considerations Heuristics for choosing A N were proposed in Bernacchia and Pigolotti (2011); O'Brien et al. (2014) for the univariate case and in O'Brien et al. (2016) in a multivariate setting.

  with K G (z) := exp -z 2 σ 2 , ∀z ∈ R d , and σ ∈ R * + . The FFT can than be used to approximate Φ(s) := F[ f ](s), and by dividing it by the transformed Gaussian kernel κ G (s) := F[K G ](s), we obtain the ECF evaluation:∆(s) = Φ(s) • [κ G (s)] -1 . (137)As computing { f (z j )} Ñ j=1 still takes O( Ñ • N ) operations, Dutt and Rokhlin (1993) suggested to use only N c < N surrounding points from {z k } N k=1 to evaluate each new grid point zj . We obtain an overall complexity of O(N c • Ñ + M log M ) which, in the case where N c < M ≤ N , is better than the original DFT formulation O(N • M

Table 1 :

 1 Average and standard deviation of the Mean Marginal Likelihood scores using confidence level

	and normalized density for 50 real flights (Real), 50 optimized flights with operational constraints (Opt1 )
	and 50 optimized flights without constraints (Opt2 ).
			(a) 1-dimensional case
	Var.		Confidence level		Normalized Density
		Real	Opt1	Opt2	Real	Opt1	Opt2
	V	0.52 ± 0.16 0.38 ± 0.14 0.15 ± 0.09 0.63 ± 0.16 0.45 ± 0.16 0.17 ± 0.10
	γ	0.54 ± 0.09 0.24 ± 0.12 0.22 ± 0.09 0.67 ± 0.09 0.33 ± 0.17 0.29 ± 0.11
	α	0.53 ± 0.06 0.08 ± 0.05 0.02 ± 0.01 0.65 ± 0.07 0.10 ± 0.06 0.02 ± 0.01
	N1	0.47 ± 0.24 0.71 ± 0.00 0.03 ± 0.01 0.57 ± 0.27 0.83 ± 0.01 0.04 ± 0.02
	Mean 0.52 ± 0.07 0.35 ± 0.06 0.10 ± 0.02 0.63 ± 0.07 0.43 ± 0.08 0.13 ± 0.02
			(b) 2-dimensional case
	Var.		Confidence level			Normalized Density
		Real	Opt1	Opt2	Real	Opt1	Opt2
	(V, γ)

Table 2 :

 2 Leave-one-out cross-validated MML scores of the training trajectories.

	Var.	Confidence level Normalized Density
	V	0.50 ± 0.19	0.60 ± 0.22
	γ	0.51 ± 0.13	0.63 ± 0.15
	α	0.51 ± 0.15	0.62 ± 0.17
	N1	0.51 ± 0.22	0.61 ± 0.24
	Mean	0.51 ± 0.21	0.62 ± 0.22

Table 3 :

 3 Average and standard deviation of the normalized density scores using Functional PCA and Least-Squares Conditional Density Estimation of 50 real flights (Real), 50 optimized flights with operational constraints (Opt1 ) and 50 optimized flights without constraints (Opt2 ).

				(a) FPCA
	Var.		Real	Opt1	Opt2
	V	0.15 ± 0.22	4.9e-04 ± 9.0e-04 2.1e-05 ± 8.1e-05
	γ	0.20 ± 0.22	9.3e-03 ± 1.5e-02 1.4e-02 ± 2.2e-02
	α	0.28 ± 0.28	1.2e-05 ± 1.8e-05 7.0e-08 ± 1.7e-07
	N1	7.6e-03 ± 6.1e-03 1.6e-02 ± 2.3e-04 1.3e-06 ± 6.7e-07
	Mean	0.16 ± 0.12	6.4e-03 ± 3.8e-03 3.6e-03 ± 5.4e-03
				(b) LS-CDE
		Var.	Real	Opt1	Opt2
		V	0.81 ± 0.13 0.63 ± 0.11 0.40 ± 0.23
		γ	0.65 ± 0.05 0.55 ± 0.10 0.53 ± 0.08
		α	0.91 ± 0.02 0.74 ± 0.03 0.68 ± 0.01
		N1	0.72 ± 0.10 0.79 ± 0.01 0.35 ± 0.05
		Mean 0.77 ± 0.05 0.68 ± 0.04 0.49 ± 0.06

where Φopt (s) = κ opt (s)∆(s) = N ∆(s) N -1 + |Φ(s)| -2 .

(126)

C.0.2. Self-Consistent Estimator

The practical problem with estimator ( 126) is that it depends on the true characteristic function Φ, which is unknown. Hence, the solutions to the fixed-point equation ( 127) was suggested by [START_REF] Bernacchia | Self-consistent method for density estimation[END_REF] to approximate Φopt :

This is justified by the fact that the optimal estimator Φopt should be very close to the true characteristic function Φ, as illustrated by the MISE criterion (108). Equation ( 127) can be transformed into a second order equation in

which admits a solution in R + provided that

When inequality (129) holds, the two possible solutions of ( 128) are

After some analysis, we can show that Φ+ is a stable fixed-point, while Φis unstable. [START_REF] Bernacchia | Self-consistent method for density estimation[END_REF] suggest to keep only the stable one, which brings us to the self-consistent estimator of the characteristic function: