Benjamin Piwowarski

Pierre Gançarski

Oihana, Olivier, Tianyi, Daria Oublier Wafa

Omar Je Veux Aussi Remercier

Lionel Rigal

Gaelle Benjamin

Florence

Bernard, Guy Vincent David

Lotka-Volterra Lv

Lorenz Lat

Attractor

MTL Multi-Task Learning SCAD Smoothly Clipped Absolute Deviation GIST Generative Iterative Sequential Thresholding DOC Damped Oscillator with Cubic dynamic

Keywords:

Tout d'abord, je tiens à remercier mes directeurs de thèse, Olivier Teste et Josiane Mothe, de m'avoir donné l'opportunité de réaliser cette thèse. Je les remercie pour leur accompagnement le long de mes trois premières années de recherche, notamment Josiane pour ses précieux conseils, de présentation de mes travaux, et d'écriture lorsqu'il s'agissait de rédiger le premier article.

Je remercie chaleuresement Olivier qui, de façon générale, a très souvent su trouver une solution en ma faveur et s'est montré encourageant du début à la fin.

Je tiens également à remercier mon encadrant insutriel, Adil Soubki, sans qui cette thèse de doctorat n'aurait pas pu exister. Il m'a laissé une grande autonomie, a su me recadrer sur ma problématique générale quand il le fallait et a su me faire prendre du recul sur mes les idées que j'ai proposées. Je

variables of a time series, the outlyingness can result in its curve shape. To highlight this shape outlyingness, we proposed to aggregate the variables of a multivariate time series in several geometric manners. Then, we used the output representation, provided by these aggregation operators, of the initial time series, as input of existing outlier detection algorithms. We empirically showed that our approach outperforms state-of-the-art methods.

Our second contribution tackles the data-driven discovery of a deterministic model underlying the dynamics between variables of a multivariate time series. We focused on the case where this (unknown) model is a system of ordinary differential equations whose solution is the function representing the time series itself. We proposed a multi-task learning algorithm to discover such a model. Each task aims at learning a single differential equation which can be coupled to the other ones. We addressed the multi-task aspect with a regularizer that enforces both sparsity within tasks and similarity between tasks. In addition, contrary to state-of-the-art multi-task regularizers, ours is nonconvex and thus provides accurate estimate of the model parameters.

We empirically showed on a benchmark of systems of ordinary differential equations that learning in a multi-task way with nonconvex sparsity outperforms state-of-the-art approaches in terms of bias and reconstruction error from the model learnt. • a p = ( k |a k | p ) 1/p : the p norm of a ∈ R d with p ≥ 1

Mathematical notations

• A p,q = k a k• q p
1/q : the p,q norm of A with p, q ≥ 1 • T : a closed interval of R

• T = t 1 , . . . , t m : an arbitrary discretization of T with m values

I.1.2 Why analyzing flight test data ?

Flight test data analysis aims both at fixing design problems and validating the aircraft performances. Once these goals have been achieved, the manufacturer has to write documentation for government certification and customer acceptance. Such a documentation entails, for instance, outlines, tables and flow-charts of aircraft performances under nominal operating conditions and emergency procedures as well. As operating conditions can vary with respect to external conditions (e.g. weather), every flight is different. Therefore the documentation cannot cover every nominal operating condition but a reference domain of known ones from which the pilot can easily infer unknown ones. "Interpolating" means adapting and inferring the aircraft behavior to a situation that does not exactly correspond to one described in the documentation but is a priori close to it. For that reason, it is necessary to understand, under various operating conditions, the behavior of the aircraft through flight test data. Furthermore, such an understanding helps flight test engineers to classify the nominal operating conditions and thereby to detect abnormal/unexpected aircraft behaviours.

I.1.3 Multiple-flight data analysis

When the behavior of the flight parameters is known by flight test engineers, multiple-flight data analysis is in principle fast. However, on the first hand, to understand both the similarities and specificities of multiple flights, the data analyses must be done in a short time (e.g. propulsive system engineers can be asked to analyze 200 flight test data in four days). On the other hand, when there is a lack of knowledge about the behavior of some flight parameters involved in a specific system, their analysis result harder and longer. As an example, propulsive systems are designed by engine manufacturers who cannot share information on the engine operating logic, making the analysis of engine parameters harder. Taking benefit from the data variability of multiple flights can help to understand the operating modes of such a system. Hence, there is a need of interpretable data-driven based methods to perform multiple flight analyses. Such analyses can reveal, for instance, an abnormal behavior of a given system which thus requires some correction on ground. Since the flight test data take the form of time series (see example in ). Grey curves correspond to normal behaviors whereas red ones correspond to abnormal behaviors. For data privacy, the vertical axes have been scaled.

I.2 Scientific motivations I.2.1 General definition of time series

In many fields (e.g. engineering, chemometrics, finance), measurement systems generate data under the form of time series [Dos Santos et al. 2016[START_REF] Boulfani | Anomaly detection for aircraft electrical generator using machine learning in a functional data framework[END_REF], Ramsay & Silverman 2006]. A time series is a sequence of random objects x indexed by some time steps t ∈ T = t 1 , . . . , t m , either naturally discrete (e.g. month, year), either resulting from a discretization of a time-continuous variable, i.e. a sampling process. When the indexed objects are univariate, and lie in R, the dataset is a collection of real-valued univariate time series (see example in Figure I.3). By extension, when the indexed objects are multivariate, and lie in R p (i.e. x is vector-valued of dimension p), the dataset is a collection of multivariate time series (see examples in Figure I.4 with p = 2), [Lafabregue et al. 2019]. When substituting R with N, one can have time series of categorical data. Hence, the notion of time series can be easily extended to non-numerical types of data like text [Baril et al. 2020].

I.2.2 The variability of univariate time series

In this thesis, we deal with datasets of time series whose indexes are the discretization of a continuous-variable. More specifically, we consider the time series as random realizations of functions depending on a continuous variable of T ⊂ R. Representing time series as functions is the building block of Functional Data Analysis (FDA). Such a representation aims to extract features containing information on the functional variability of the time series and can serve for dimensionality reduction, smoothing, feature extraction, prediction as well as visualization [Ramsay & Silverman 2006, Ferraty & Vieu 2006]. In time steps in T = t 1 = 0, . . . , t 1000 = 10 , generated by x(t) = t sin(t) + t 2 cos(t)+ for t ∈ T = [0, 10], where is a random variable with distribution N (0, σ = 0.25).

I.2.3 The second variability of multivariate time series

Although, the only difference between univariate and multivariate time series is their dimension at each time step, in multivariate time series the p variables can be correlated between them along time. This correlation can result from a mechanism that is deterministic (e.g. a system of differential equations) or stochastic (e.g. a vector auto-regressive model, stochastic differential equations). Therefore, in addition to the temporal variability of each variable individually, multivariate time series have a variability across its p variables.

Example I.1. In Figure I.4, we illustrate the correlation between the variables with a synthetic dataset of bivariate time series generated as follows:

• For each sample i ∈ {1, . . . , n = 10}, the time steps Ti = {t 1i , . . . , t m i } are randomly chosen from U(0, 10) with m i = 250.

• The first variable for sample i is generated as x 1i (t) = t sin(t) + t 2 cos(t) + where is a random under the distribution N (0, σ = 0.25) and for t ∈ Ti .

• The second variable is generated according to x 2i = A i sin(2πf 0 t) + where is independently generated as in x 2i , and A i is drawn from U(2, 3).

Simple algebraic manipulations show that x 1 and x 2 are related by x 1 (t) = tx 2 (t) + t 2 x 2 (t + π 2 ). Thus there is a nonlinear correlation between x 1 and x 2 . Furthermore, one can show, by computing analytically ẋ1 and ẋ2 , that the non-noisy version of the functions x 1 and x 2 (without ) is the solution of the following system of ordinary differential equations:

     ẋ1 (t) = A(1 -t 2 )x 1 (t) + 3Atx 1 (t + π 2 ) ẋ2 (t) = Ax 1 (t + π 2 ) (I.1)
with A = 0 and the initial condition x 1 (0) = x 2 (0) = 0.

This example serves as an illustration of a (functional data based) generating process of multivariate time series. Of course, in practice the generating process is unknown but it shows how the correlation between the p variables along time can be underlied. This example also motivates to retrieve a functional mechanism that underlies a time series dataset, this is the first motivation of FDA. Example I.1 illustrates how a deterministic model can characterize the correlation between the variables of a multivariate time series. Accessing to the system of differential equations underlying a real dataset can be of huge interest for engineers in real life applications. For instance it can serve as a simulation model or as a representation of the phenomenon dynamics. For real datasets, such kind of model is often unknown by the practitioner. Recently the machine learning community has proposed data driven methods to retrieve the model underlying the dynamics of a noisy dataset of time series [START_REF] Brunton | [END_REF][START_REF] Schaeffer | [END_REF], Schaeffer 2017, Mangan et al. 2017].

In Chapter III, we interest in the data-driven discovery, in closed form, of a system of differential equations that models the dynamics of a multivariate time series. Such an approach is built upon the recent framework of sparse identification on nonlinear dynamics [START_REF] Brunton | [END_REF].

I.2.4 General definition of outliers

The concept of outlier has been defined in [START_REF] Hawkins | [END_REF] as "an observation which deviates so much from other observations as to arouse suspicions that it was generated by a different mechanism". This definition is not specific to time series and is now well accepted. According to the definition in [START_REF] Aggarwal | [END_REF], an outlier is defined as a data point that is highly different from the others, based on some measure. Therefore, outliers can be seen as observations of a dataset that do not follow an expected

CHAPTER I: GENERAL INTRODUCTION

behavior. Depending on the context, outliers are also referred as anomalies in [START_REF] Erfani | [END_REF], Chandola et al. 2009, Rousseeuw & Hubert 2018, Liu et al. 2008], novelty in [Markou & Singh 2003, Pimentel et al. 2014] or more rarely singularities. In certain cases, outliers are of analysis interest since they can reveal a design fault in the system. For instance, detecting outliers can help engineers working on a system requiring a high safety level (e.g. aircraft) to detect unexpected system operating modes during the system test phase.

I.2.5 Outliers in time series

Outlier in time series has different meanings in the literature. In many papers, it is assumed that the dataset is made up of a single multivariate or univariate time series. Hence, detecting outliers in such a dataset amounts to find points of the time series that deviates from a regular pattern in a given sense. Oppositely, a dataset can entail multiple time series. Detecting outliers in such dataset amounts to find the time series that deviate from the other ones. In this thesis, we deal with datasets containing multiple multivariate time series.

I.2.5.1 Outliers in a univariate time series dataset

In a dataset of univariate time series, a sample can result as an outlier in several ways. Assuming that a function u of a temporal continuous variable, t ∈ T (see example Figure I.3) underlies the dataset, according to the taxonomy of [Hubert et al. 2015], the deviation of an outlier can be:

• horizontal-shift: the time series was generated by the same process as inliers up to a time translation, i.e. u horizontal = u(t -τ ) where τ is random with nonzero mean.

• in the magnitude: the time series was generated by the same process as inliers up to a magnitude shift, i.e. u magnitude = u(t) + a where a is random with nonzero mean.

• in shape: the time series was generated by a different process as inliers and does not stand out at any time point.

In Figure I.5, we give an illustration for each of these types of outlier. with τ drawn from N (6, 0.1). The magnitude outlier, in red, was generated by u magnitude with a drawn from N (-30, 0.1). The shape outlier, in green, was generated according to u(t) + sin(2πf 0 t) with f 0 = 7.

I.2.5.2 Outliers in a multivariate time series dataset

In order to reduce the analysis of a dataset of p-dimensional multivariate time series, one can split it into p datasets of univariate time series. This boils down to omit the correlation between the p variables of each time series.

However, such a correlation can be of great matter to detect outliers. Indeed, when considering each x 1 , . . . , x p independently, the deviation of an outlier in, say x k , can be irrelevant compared to its deviation in the joint representation of, say x k and x k . This is mentioned in [Hubert et al. 2015, section 2.2] as "None of the outliers [. . . ] will necessarily be outlying in one of the marginals, hence it is not sufficient to analyze the components separately" and illustrated in Then, comes the question of how to consider the p variables jointly across the time variable. To answer such a question, in Chapter II, we propose to represent multivariate time series as multidimensional path from which we extract geometric features that highlight their outlyingness.

I.3 Contributions

To tackle the detection of outliers in multivariate time series, we made two contributions. The common idea of our two contributions is to represent time series as functions over a time-continuous variable. In the first one, a multivariate time series is represented as a multivariate path (or a trajectory) from which shape-based features are extracted, under the form of a univariate time series, and then input into an outlier detection algorithm. In this manner, multiple kinds of outliers can be detected. In our second contribution, a multivariate time series is seen as the solution of an unknown system of differential equations that is recovered by a multi-task learning algorithm.

Our algorithm returns the analytic form of the differential equation which is useful to model and understand the underlying mechanism of a multivariate time series.

I.3.1 Overview

I.3.1.1 Outlier detection with shape-based features

Our first contribution is to propose a shape-based feature extraction from multivariate time series and is based on the multivariate functional data framework [Lejeune et al. 2020a, Lejeune et al. 2020b]. Multivariate functional data refer to a population of multivariate functions generated by a system involving dynamic variables depending on a continuous variable (e.g. time). Outlier detection in such a context is a challenging problem because, as we mentioned along Chapter I, both the individual behavior of the variables and the dynamic correlation between them, are important. To identify the outliers, recent work has focused on multivariate functional depth [START_REF] Claeskens | [END_REF], Dai & Genton 2019, Kuhnt & Rehage 2016] which assigns an outlyingness score to each variable independently and then sums over the scores. However, such an approach fails when the outlyingness manifests in the curve shape rather than in the curve magnitude. To remedy it, the curve geometry has to be considered across its variables, jointly, rather in each Data-driven approaches to automatically discover the underlying SODE have been made possible thanks to sensor technology which can collect large amount of data [START_REF] Brunton | [END_REF], Long et al. 2018, Schaeffer et al. 2018, Bhat & Rawat 2019]. Nevertheless, state of the art approaches are based on single-task learning that means each component of the SODE is learned independently [START_REF] Argyriou | [END_REF] with sparse linear regression, where the sparsity is reflected in the convex 1,1 penalty. This leads to SODEs that weakly identify the underlying phenomenons since in reality, equations are related to each others. Moreover, it turns out that the convexity of the penalty involved in the learning criterion results in a SODE which is biased with respect to the true one [Fan & Li 2001, Zhang 2010] -Gil & Romo 2014, Hubert et al. 2015]), in chemometrics (to determine the nature of an active substance produced by a chemical process based on near-infrared spectra data [Hubert et al. 2015]), or in air pollution studies (to detect highly contaminated locations in urban areas [Torres et al. 2011]). In these fields, the data are typically functional and exhibit outlying behavior. Moreover, several parameters should be simultaneously recorded to accurately understand the studied process.

Hence, outlier-detection methods should be specifically designed for multivariate functional data. When the variables are cautiously selected by a domain expert, the outlying behavior can be detected through the potential correlation between them.

The correlation between the p variables is important in multivariate functional data because it can reveal an outlying behavior of the underlying process, as discussed in [Hubert et al. 2015]. In analyzing each variable implies that the potential correlation between the variables is not considered, According to the definition by [START_REF] Aggarwal | [END_REF], an outlier is defined as a data point that is highly different from the others, based on some measure.

Such a point often contains useful information regarding the abnormal behavior of the system described by the data. Moreover, if the data dimension is high, the data are likely to be more scattered in the space (i.e. curse of dimensionality), and therefore, the probability that the outliers are scattered is higher. Hence, outlier-detection is inclined to the curse of dimensionality as other classification tasks that assume well-balanced classes. However, regarding some typical algorithms for classification (e.g. logistic regression)

and clustering (e.g. K-means and mean-shift), the rarity and scattering of

CHAPTER II: SHAPE-BASED OUTLIER DETECTION (a) -2 -1 0 1 2 0.0 0.2 0.4 0.6 0.8 1.0 -2 -1 0 1 2 x 1 x 2 t -2 -1 0 1 2 -2 -1 0 1 2 (b)
x the curve variables by using differential analysis [START_REF] Claeskens | [END_REF]]. This refers to adding derivatives or integrals (computed with respect to t) for each initial variable. Hence, curve shape provides information regarding "hidden outlying features" of the curve variables and the outlying relationship between them. However, as mentioned previously, the joint analysis of the p variables becomes complex as p increases. In the present study, we address CHAPTER II: SHAPE-BASED OUTLIER DETECTION this problem by using differential geometry. Specifically, we use aggregation functions (here termed as mapping functions) of the variables. Thereby, we implicitly consider the correlation of the variables through geometrical characterizations of curve shape. In contrast with current functional-outlier detection methods, which consider curve shape differently and only base the final detection on the resulting depth values (Section II.2), we use both functional curve-shape features and state-of-the-art outlier-detection algorithms e.g. Isolation Forest [START_REF] Liu | [END_REF]]. Thus, the originality of the proposed approach lies in the shape characterization of the initial curves through the proposed mapping functions whose output representation is input to stateof-the-art outlier-detection algorithms.

The mapping functions that we use refer to aggregation functions that enable capturing curve-shape features, such as curvature, length (i.e. perimeter of a shape), or tangential velocity, and consider all the variables, as a curve is viewed as a path. More precisely, a mapping function aggregates the variables through different interpretable combinations of the derivatives of the variables. Mapping functions have been used in shape analysis [Srivastava & Klassen 2016], that is, for curves lying in a two-or threedimensional space (e.g. extracted from images), but not in the detection of multivariate functional outliers.

To capture the potential outlyingness of the curves through their shape, we propose mapping functions among those used in differential geometricmethod in shape analysis [Srivastava & Klassen 2016]. These functions map multivariate to univariate curves; however, for accurate computation, they require the curves to be smooth. Although this is the case for multivariate functional data, raw data are often noisy when sampled, and we use the functional-data representation to recover the smooth version of the curves.

Then, using the proposed mapping functions, we map the functional representation (in the form of a curve) so that some of its shape features capture curve outlyingness. Finally, based on this new representation, we use outlierdetection algorithms to assess the outlyingness of each sample and determine a threshold for flagging outliers.

Our contributions are summarized as follows:

(i) We propose an end-to-end method for detecting outliers through their curve shape, which is characterized by geometrical transformations.

The method is based on the functional representation of the data.

( (ii) It decreases monotonically and vanishes as the sample moves away from the center (up to infinity).

(iii) It does not depend on the dataset scale.

Therefore, given an outlyingness threshold, samples with a depth value close to 0 can be flagged as outliers. This type of measure has been extended

to functional data and used for classification [START_REF] Cuevas | [END_REF], ranking [Fraiman & Muniz 2001, Cuevas et al. 2006], as well as outlier detection [Febrero et al. 2008 with depth value lower than the 5%-quantile of this distribution is likely an outlier). [Febrero et al. 2008]

II.2.2 Depth-based multivariate functional outlier detection

Depth-based outlier detection methods for multivariate functional data are more recent. [START_REF] Claeskens | [END_REF] generalized any given univariate functional depth to the case of multivariate functional data. This corresponds to a weighted sum of a given univariate functional depth applied to each variable (x 1 (t), ...x k (t), ...x p (t)) pointwise and then integrated over T . The authors also discussed the selection of the weight function. As a special case, in [Ieva & Paganoni 2013] proposed the multivariate band depth by using the modified band depth as the given univariate functional depth; the weights associated to the variables are constant with respect to t. [Hubert et al. 2015] noted that the generalization proposed in [START_REF] Claeskens | [END_REF] does not allow to detect all shape outliers since low-depth samples stand near the boundary of the dataset but may not be outliers. Conversely, high-depth samples may present outlyingness in their curve shape because, pointwise, the curve does not exhibit any significant deviance in each variable, as this generalization is the sum of the individual univariate functional depths. To address this, the entire shape of the curve should be considered.

A few studies incorporate curve shape into a multivariate functional depth measure. Recently, [Kuhnt & Rehage 2016] proposed the functional tangential-angle (F U N T A) pseudo-depth, which considers curve shape based on the intersection angles of the centered variables (i.e. the variables are scaled so that their integral over T values is 0). More precisely, for each variable, F U N T A computes the intersection angles of a given sample x ik with all the other samples x jk ∀j = i, and then averages these angles over the number of intersection angles of x ik and over the variables k = 1...p. Thus, F U N T A separately considers the shape for each variable with respect to t, but not the shape between the p variables.

More recently, [Dai & Genton 2019] proposed the directional outlyingness measure (Dir.out), which considers curve shape through the weighted pointwise direction in R p of the vector X(t) toward the median of the distribution of X(t). The purpose of the weights is the up-weighting of the directions in which the outlyingness of X(t) is likely to appear. In contrast with the aforementioned multivariate functional depths, which provide a score in [0, 1], the Dir.out depth returns a vector in R p ×R + corresponding to the concatenation of the mean directional outlyingness (in R p ) and the total variance of the directional outlyingness (in R + ). A final outlyingness score is computed as the robust Mahalanobis distance between this vector and a mean vector of the same type computed on a subset of independent samples. Then, the upper tail of this distance distribution is approximated by an F -distribution, and the outlyingness threshold is defined as a high-probability quantile of this F -distribution. Hence, unlike in other multivariate functional depths, the outlyingness threshold provided by the Dir.out approach is not data-driven,

CHAPTER II: SHAPE-BASED OUTLIER DETECTION

as it is based on the (approximately) true distribution of the outlyingness scores. However, in this approach, the parameters should be tuned by simulation and are difficult to interpret beyond the statistical framework.

Multivariate functional depths are related to curve shape through the individual behavior of the curve variables. In our contribution, we adopt a different approach, as we view a curve as a path in R p and process it as a geometrical shape with mapping functions. Indeed, the curve shape is not only reflected by the curve variables, individually along t, but also by the relationship between them. Our mapping functions consider the curve shape along t thus highlighting outlying correlation features between the variables.

As all the aforementioned multivariate functional depths yield an outlyingness score with unknown distribution (except for Dir.out), an outlyingness threshold can be computed from the resulting empirical distribution of the depth values through a bootstrap procedure as in the univariate case [Febrero et al. 2008]. It can also be computed from a training dataset based on the receiver operating characteristic (ROC) curve.

In the experimental study (Section II.5), we use the F U N T A and Dir.out functional depths as baselines because they have been demonstrated to be promising for outlier detection in multivariate functional data by regarding outlyingness as a curve-shape feature.

II.2.3 Geometry-based functional outlier detection

Representing functional data in a geometric framework is a recent idea, and few studies have considered such representations for outlier detection. Recently, [START_REF] Xie | [END_REF] they computed on the original dataset (although the computational methods are quite different), the authors proposed a method for constructing a functional boxplot for each of the three components so that outlying features may be identified, and outliers may therefore be detected. [Xie et al. 2019] extended this method to multivariate functional data and added other components such as shape orientation (reflecting rotational variability). They additionally provided useful visualization techniques for identifying outlying features, in fact, they only focused on the bivariate, p = 2, and trivariate, p = 3, cases, which are shape data extracted from images. However, when the size of the dataset and the number of variables p increase, this method is computationally costly, as the shape-based component-extraction procedures include several continuous optimization problems. Moreover, in these studies, the outlier-detection methods are based purely on the empirical distribution (through the functional boxplot) of the proposed geometrical features, whereas in our work we geometrically aggregate the dimensions of the curves, resulting in univariate functional data and subsequently detect outliers from this new data representation. The latter can be seen as implicit non-parametric learning of the inlier distribution based on the functional data mapped to a geometric curve feature. Hence, we take advantage of both the geometrical mapping and the outlier-detection algorithm.

CHAPTER II: SHAPE-BASED OUTLIER DETECTION

II.3 Background in functional data

This section is concerned with the handling of high-dimensional vectors of discrete noisy measurements that can be represented as smooth continuous functions; moreover, we discuss how such representations can be achieved.

We use the notation introduced in Mathematical notations. The functional data representation is twofold:

(i) As the xi s are smooth functions, the reconstructed data are noiseless.

(ii) The reconstructed data are "aligned" in the sense that two reconstructed sample values x1 (t j ) and x2 (t j ) t j at t j are comparable, as they refer to the same evaluation point t j ∈ T . This is not the case in raw data because one can have t m 1 = t m 2 (the curves can be sampled on different grids).

II.3.1 Functional-data representation

The first step in FDA is to approximate an unknown smooth function x i : t → R p , which underlies the sample i, by another smooth approximation function xi (t), ∀t ∈ T , through m i discrete noisy measurements x i (t 1 ), ..., x i (t m i ); this is referred to as the functional approximation step. Its purpose is to remove the noise, thus allowing accurate evaluations of some derived functions, such as combinations of derivatives and integral functions. This is necessary in our case, as the proposed mapping functions correspond to combinations of derivatives and integrals.

We should first select a functional representation as an approximation function. As a function is intrinsically infinite-dimensional, in FDA, it is commonly assumed that the underlying function can be approximated by a finite linear combination of non-linear basis functions. Such an approximation is called a basis expansion function [Ramsay & Silverman 2006]. We assume that x ik , the k-th variable (hence a univariate function) of x i , is to be approximated. The intuition behind the basis expansion is to combine a small number of "specific functions" (a set of given functions), each of which can capture some local features of the underlying function x ik , so that x ik could be recovered with a small approximation error. This approximation function can be formulated as

∀t ∈ T , xik (t) = L ik l=1 α ikl φ l (t) = α ik φ(t) (II.1)
where φ(t) = {φ l (t)} 1≤l≤L ik is a vector of orthonormal basis functions at t for some L ik ∈ N * (referred to as the basis size) with fewer basis functions than sampled observation points (L ik m i ), and α ik = {α ikl } 1≤l≤L ik is the coefficient vector, the element α ikl of which is the importance of the l-th basis function.

Another choice of functional representation in FDA is to use non-parametric smoothing [Ferraty & Vieu 2006], which achieves a similar approximation, but its form is less tractable than that of the basis expansion function, for instance, to compute derivatives.

According to Equation (II.1), one should select (i) the basis {φ l } 1≤l≤L ik and (ii) the basis size L ik .

The choice of the basis is data-dependent. As suggested by Ramsay and Silverman [Ramsay & Silverman 2006], when the data are smooth and periodic, the Fourier basis should be selected; when the data are smooth, a B-spline basis is suitable. A B-spline is a piecewise-polynomial function of order at least three [De Boor 1978]. It is not exactly an orthonormal basis but since we only exploit the differentiability that B-splines induce, here we omit the non-orthonormality. If one requires orthonormality and smoothness, the functional principal components basis is a good solution [Ramsay & Silverman 2006].If the data have irregularities, a wavelet basis CHAPTER II: SHAPE-BASED OUTLIER DETECTION should be preferred [Nason 2008]. See [START_REF] Ramsay | [END_REF] for other examples and details on the choice of the basis according to the data. In turn, the choice of the basis size parameter L ik depends on the selected basis. An inappropriate choice of the basis results in requiring a large L ik because each basis function will focus on an irrelevant part of the data variability (low bias and high variance or, high bias and low variance); the worst case is to capture the noise, leading to over-fitting [Ramsay & Silverman 2006]. By contrast, an appropriate choice of the basis functions results in a small L ik , that is, the basis is sufficiently rich to approximate an unknown function using few functions. Subsequently, once a suitable basis is selected, the biasvariance trade-off should be considered. This refers to the balance between the approximation error and a reasonable L ik [Ramsay & Silverman 2006]. Such a balance is generally achieved by a grid search by cross-validation for each sample i and variable k. When φ(t) and L ik are specified, a computing method is required to estimate the coefficient vector α ik , which is introduced in the next paragraph.

II.3.2 Functional-data fitting

The linearity of the basis expansion function with respect to the coefficient vector α ik enables its efficient estimation (assuming the data were sampled with a Gaussian noise ij , that is,

x ik (t ij ) = xik (t ij ) + ij , where ij is inde- pendent of xik (t ij )
) by minimizing the least-squares criteria:

J (α ik ) = m i j=1 (x ik (t ij ) -xik (t ij )) 2 (II.2)
or equivalently, with vector notation,

J (α ik ) = x ik (t i• ) -Φ ik α ik 2 2 (II.3)
where t i• is, by abuse of notation, the vector containing all the samples

{t ij } 1≤j≤m i and Φ ik = (φ l (t ij )) 1≤j≤m i ,1≤l≤L ik is the m i × L ik matrix contain-
ing all the L ik basis functions evaluated at the observation points. Thus, Φ ik is a discretization over t i• of the vector of orthonormal basis functions φ(t) in Equation (II.1). As L ik m i and Φ ik has all its columns linearly independent, by the orthonormality of the basis functions (and thus orthonormality of Φ ik ), Φ ik Φ ik is invertible. Hence, equating the gradient of J to 0 with respect to α ik leads to the following minimizer:

α * ik = arg min α ik J (α ik ) = (Φ ik Φ ik ) -1 Φ ik x ik (t i• ) (II.4)
which is known as the classical least-squares solution [START_REF] Hastie | [END_REF]].

However, as the data are fitted according to the basis functions, the smoothness of xik depends greatly on the noise influence on the basis functions. Consequently, xik may lack smoothness and overfit the data. To analyze such a noise influence, one can compute the derivative of xik , which is "excessively" variable if a large amount of noise remains in the approximation function. To ensure smoothness, the least-squares criteria should be minimized by penalizing the derivative(s) of xik with an amount λ > 0 as follows:

J λ (α ik ) = m i j=1 (x ik (t ij ) -xik (t ij )) 2 + λ T (D q xik (t)) 2 dt (II.5)
where D q = d q (•) dt q is the q-th derivative of xik (t). More generally, D q can be any linear combination of derivatives of x ik , that is, a linear differential operator [Ramsay & Silverman 2006]. A penalization term including derivatives is also known as a roughness penalty. The parameter λ is arbitrary and can be computed by cross-validation. This is detailed in Section II.5.2. Equation (II.5) can be written using vector notation as follows:

CHAPTER II: SHAPE-BASED OUTLIER DETECTION J λ (α ik ) = x ik (t i• ) -Φ ik α ik 2 + λα ik R ik α ik (II.6)
where

R ik = ( T D q φ j (t)D q φ m (t)dt) 1≤j≤L ik ,1≤m≤L ik is a L ik × L ik positive
semi-definite matrix. The matrix R ik contains the inner products of the q-th derivative of the L ik basis functions. This matrix can be computed provided that the q-th derivative of the basis functions exists. In practice, it is common to choose q = 1 or q = 2 (i.e. to penalize the velocity or acceleration of xik , or a combination of both).

As J λ remains quadratic with respect to α ik , approximating xik with a roughness penalty is equivalent to ridge regression [Hoerl & Kennard 1970, Hastie et al. 2009]. Thus, the penalty term allows xik to (i) be smooth, as defined by the operator D q and, (ii) avoid over-fitting by pushing the coefficient vector near 0.

Equating the gradient of J λ to 0 with respect to α ik leads to the following minimizer [START_REF] Hastie | [END_REF], Ramsay & Silverman 2006]:

α * ik,λ = arg min α ik,λ J λ (α ik,λ ) = (Φ ik Φ ik + λR ik ) -1 Φ ik x ik (t i• ) (II.7)

II.3.3 Approximation functions as building blocks

Once the coefficient vectors have been estimated for the p variables of the n samples (with or without penalization), we can consider the approximations xik to be smooth multivariate functions that well recover the underlying functions. Although these functions can be theoretically evaluated at an infinite number of points in T , in practice, there are two methods to handle the approximations computationally (e.g. to compute derived functions such as derivatives and integrals):

(i) The first method is to compute the derived functions based on the basis functions. As the basis functions are known analytically, their derived functions can also be obtained analytically. Thus, by the linearity of the basis expansion, one can easily obtain the derived functions of the approximation functions (the integral and derivative are linear operators). We illustrate this using the k-th derivative of the approximation function. We assume that an unknown function x is approximated by

x through a basis expansion with a basis size L (in Equation (II.1)), provided that the k-th derivative {D k φ l (t)} 1≤l≤L of the basis functions exists, and the coefficient vector {α l } 1≤l≤L is available (or has been estimated as in Equation II.4). The k-th derivative of x with respect to

t is D k x, where ∀t ∈ T , D k x(t) = D k L l=1 α l φ l (t) = L l=1 α l D k φ l (t) (II.8) (ii)
The second method is to estimate the underlying functions by evaluating all the approximation functions on the same grid T . Thus, from these estimates, one can compute derived functions, such as integral or derivatives, using numerical methods, such as quadrature or finite difference schemes, respectively [Stoer & Bulirsch 2013]. These methods are easy to implement, but they do not consider the basis functions and require that the arbitrary grid be sufficiently fine (so that the approximation functions are evaluated at a large number of observation points).

Thus, if the derivatives of the basis functions are known analytically (as is the case for B-splines, Fourier basis functions, etc.), the derivatives of x are also known and do not need to be estimated from the raw data or from the smooth reconstructions by a data-driven method such as finite differences.

Equation II.8 demonstrates the flexibility of the linear basis expansion for computing derived functions in FDA. Then, a derived function, for instance D 1 x, can be evaluated on an arbitrary grid. Such an approach is different from estimating the derivatives from an evaluation of x on the grid by using finite differences.

The first method is safer than the second because the analytic form of the basis functions is fully considered, and therefore the corresponding derived functions can be obtained accordingly. For instance, if the basis functions φ l are B-splines (which are piecewise polynomial), we know the analytic form of D 1 x, as D 1 φ l results in a piecewise polynomial as well. Thus, the evaluation of D 1 x by the first method provides more accurate estimates of D 1 x (which is unknown) than numerical methods applied to x evaluated on a fine grid of T .

In the following section, we suggest some mapping functions for capturing functional outlyingness in the detection process. These mapping functions may have a complex analytical form because they involve several derivative (first and second order derivatives, as well as integral functions). Therefore, it is mandatory to have accurate evaluations of derivative functions, and we follow the first method in the computational experiments since we use B-splines and Fourier basis functions, whose derivatives are known.

II.4 Shape-based representation of curves

We regard a multivariate curve as a path lying in a p-dimensional space, R p , and derive mapping functions (aggregation functions over the variables), established in differential geometry, to capture shape features of the curves (e.g. length, velocity, or curvature) so that outlying features may be detected. These mapping functions have been used in shape analysis to extract features based on the edge (a bivariate curve) of an object in an image [Srivastava & Klassen 2016].

In this section, we investigate several mapping functions that enable the detection of multivariate functional outliers from the shape they exhibit in R p . Such mappings jointly consider the p variables, as they aggregate, in several ways, some derivatives (with respect to t) of the curve variables.

Hence, the individual and collective variations of the variables are considered.

These mapping functions take each data sample, represented by its smooth approximation function xi , as input and return a univariate curve (i.e. the resulting aggregation) reflecting certain shape features. Hence, they provide a means to "summarize" the shape of a multivariate curve, in the sense given by the mapping function, and reduce the number of functional variables to one. The univariate function returned by a mapping function is then fed into an outlier-detection algorithm; this is detailed in Section II.5. In the sequel, we simplify the notations by referring to a functional-data sample as an arbitrary curve x = x 1 , . . . , x p instead of xi = xi1 , . . . , xip .

II.4.1 Arc-length mapping

The arc-length mapping function enables analyzing the length of a curve between two points in The length of the curve between two observation points t 0 (darkgrey dot) and t (white dot) is defined as the sum of infinitesimal length elements ds(t 0 )...ds(t) along the curve (red diagonal arrows) for all t. The crossed-circle dots represent such points between t 0 and t.

Hence, the arc-length maps an original functional-data sample to univariate functional data that represent the increases in the cumulative length of the underlying curve from the starting-point x(t 0 ) = ((x 1 (t 0 ), x 2 (t 0 )) to an arbitrary point x(t) = ((x 1 (t), x 2 (t)) for t > t 0 . Figure II.2 shows that the length of a bivariate curve between x(t 0 ) and x(t) is the infinite sum from t 0 to t of infinitesimal length elements ds(•) (aka integral), corresponding to an infinitesimal length element in each direction (x 1 and x 2 ) in R 2 . This mapping always returns a positive increasing function, as it computes the cumulative length of the initial curve. Moreover, the arc-length mapping function is not influenced by a warping (i.e. a horizontal deformation) of the curve 2 which is usual in FDA to compare curves sampled on different grids [Srivastava & Klassen 2016]. This mapping function can discern functional samples with a shape of different size, which is a global shape feature.

2 Let α(•) be a differentiable warping function i.e. a monotone nondecreasing function defined in T → T .

The arc-length mapping function on a warped functional datum x is equal to the arc-length mapping function on the initial unwrapped functional datum: Thereby, the detection of functional outliers can be improved when their underlying curve is longer or shorter than those of the other samples. For instance, an isolated outlier, which exhibits a peak for a small part of T , induces a sharp increase in its curve length, whereas the length of other curves increases more slowly.

s(α(t)) = t t0 D 1 x(α(u)) du = t t0 D 1 x(α(u)), D 1 x(α(u)) 1/2 du = t t0 D 1 α(u) D 1 x(α), D 1 x(α) 1/2 du, and as D 1 α(u) = dα du , we have t t0 D 1 α(u) D 1 x(α), D 1 x(α) 1/2 du = t t0 D 1 x(α), D 1 x(α) 1/2 dα = t t0 D 1 x(α) 2 dα, which implies that s(α(t)) = s(t).

II.4.2 Velocity mapping

The velocity mapping function enables analyzing the instantaneous variations of the curve with respect to t. It has a simple interpretation when t corresponds to a time instant. In this case, velocity measures how fast a point moves on the curve. More generally, it can be interpreted as the norm of the projection of the curve onto D 1 Y (t), the tangent vector to the curve at t. In Figure II.3, the velocity mapping at t of a bivariate curve is shown as the 2 -norm of the tangent vector D 1 x(t) (vector of the first-order derivatives of the curve variables x 1 and x 2 ). It is defined as Curvature is defined to be the inverse of the radius of the osculating circle. In this example, in a neighborhood of the curve at t 1 (dark-grey dot), the tangent vector D 1 x(t 1 ) has almost the same direction; hence, the osculating circle has a large radius (r(t 1 ) = 1 κ(t 1 ) ), resulting in a small curvature. In a neighborhood of the curve at t (white dot), the tangent vector D 1 x(t) quickly changes direction; hence, the osculating circle has a lower radius, that is, a higher curvature than at t 1 . and is related to the arc-length mapping by ψ(t) = ds dt , or conversely, by s(t) = t t 0 ψ(t)dt; however, these mappings capture different features. Indeed, the arc-length mapping outputs an increasing function and thus "memorizes" the local variations of the curve as t increases, whereas the velocity mapping characterizes the local variations (i.e. pointwise) with respect to t. The function returned by the velocity mapping may be regarded as a measure of the variation of the arc-length mapping. Thus, the velocity mapping can be used to identify the local outlyingness of a sample that is to say an isolated outlier.

ψ(t) = D 1 x(t

II.4.3 Curvature mapping

Curvature is a notion that relates to how "bended" a curve is, or geometrically, the degree to which a curve deviates from the tangent line at a given point. An alternative interpretation concerns the radius of the osculating circles. At a given point t, a smaller radius of the osculating circle implies larger curvature. In fact, the radius of the osculating circle is equal to the inverse of the curvature at this point. The bivariate curve in Figure II.4 shows that at a neighborhood of t 1 , where the tangent vector D 1 x(t 1 ) has almost constant direction, the osculating circle has a larger radius r(t 1 ) than the radius of the osculating circle at a neighborhood of t where the direction of the tangent vector D 1 x(t) changes quickly. Thus, the curvature mapping function allows analyzing the change of direction of the curve with respect to t. Indeed, if the curve is a line, curvature is constant, and the curve directions remain constant as well. Curvature is defined as [Srivastava & Klassen 2016]

κ(t) = D 1 ( D 1 x(t) D 1 x(t) 2 ) 2 D 1 x(t) 2 (II.11)
or equivalently,

κ(t) = D 1 x(t) 2 2 D 2 x(t) 2 2 -D 1 x(t), D 2 x(t) 2 D 1 x(t) 3 2 (II.12)
We now provide insight into the definition of κ in Equation (II.11). D 1 x(t) D 1 x(t) 2 is the direction vector (i.e. the normalized tangent vector); therefore,

D 1 D 1 x(t)
D 1 x(t) 2 is the rate of change of the direction vector, and the normalization D 1 x(t) 2 relates to the rate of change of the direction with respect to the tangent vector. Consequently, the curvature mapping can detect functional outliers with a curve that exhibits a differently bended shape than those of the other samples.

II.5 Experimental study

We conducted an experimental study on real and synthetic datasets to demonstrate the effectiveness of the proposed mapping functions in improving outlier detection in multivariate functional data. The detection performance was evaluated in terms of the true detection rate (i.e. the proportion of outliers correctly detected), false detection rate (i.e. the proportion of outliers falsely detected), and area under the ROC curve (AU C).

II.5.1 Real data

II.5.1.1 ECG data

We tested the proposed approach on the real dataset used in [Dai & Genton 2019]. The dataset consists of electrocardiogram (ECG) time series of the electrical activity (voltage) of heart changes [Goldberger et al. 2000]. Such data can reveal abnormalities in heart activity. The time series are univariate and were labeled by cardiologists as abnormal or normal. This dataset has been used for time-series classification [START_REF] Wei | [END_REF].

The ECG dataset consists in n = 810 time series including 208 abnormal and 602 normal cases. All the time series have an equal size of m i = 86. In contrast with [Dai & Genton 2019], who only considered the time series between the time stamps t = 6 and t = 80 to avoid boundary effects, we considered the entire time series to demonstrate the robustness and applicability of the proposed approach. Dai and Genton also augmented the univariate time series to multivariate by adding the first and the second derivatives. We did not follow this, as in our proposal, these aspects are considered (e.g. velocity mapping in Equation (II.10) or curvature mapping in Equation (II.11)); rather, we added the squared time series. Indeed, power is proportional to the square of voltage. Thus, in terms of interpretability, this data augmentation appears to be more relevant than the second derivative of voltage. We applied the same multivariate functional data augmentation to all ECG data II.5. EXPERIMENTAL STUDY experiments and for all methods. We did not apply the derivative augmentation, as this would bias the interpretation of the results, that is, it would not be possible to discern whether the results were due to the specific augmentation or to the method. This would be of interest if the focus was specifically on the ECG data, but here, we use it as a real dataset example.

As in [Dai & Genton 2019], to obtain a rare class of samples representing outliers, we randomly created a partition of 400 samples (i.e. the training set) out of the 810 samples by parameterizing the contamination level (i.e. the rate of abnormal samples) in this partition to 5%, 10%, 15%, 20%, and 25%. Then, for each contamination level, we evaluated the proposed method on the 410 remaining samples (i.e. the test set).

II.5.1.2 Pen-digits data

We tested the proposed method on the real dataset consisting of n = 10, 992 bivariate time series representing pen digits (PenDig) [START_REF] Dua | [END_REF]].

The digits are labeled according to their class, from '0' to '9'. Each digit has m i = 8 observation points regularly sampled on both the horizontal and vertical coordinates. As this initial dataset cannot be considered highdimensional, we upsampled it by linear interpolation to m = 200 on the two coordinates before fitting the approximation functions.

To simulate the outlier classes, we considered a single digit to be the outlier class, and the nine other classes to be the inlier class, as in [Ruff et al. 2018].

The training set was generated using 75% of the entire dataset with a contamination level equal to c = 5% (i.e. 5% of the training set are outliers).

Each digit was separately considered the outlier class, and thus the experiment was conducted in 10 independent ways. Then, for each case of outlier class, we assessed the proposed method on the test set.

II.5.1.3 Airbus flight test data

We also tested our method on a real dataset provided by an Airbus flight test expert, of n = 7, 118 four dimensional multivariate time series wherein The ECG and PenDig datasets were not used to assess the same properties of the proposed method. The ECG data were used to demonstrate the robustness of the proposed method with respect to different contamination levels for some given outliers, whereas the PenDig data were used to assess the detection performance for different outliers and a given contamination level.

Thus, we only compared these two in terms of performance, in the comparison of the various methods in Section II.5.8. The AFT dataset served as an application example of outlier detection in flight test data.

II.5.1.4 Synthetic data

We simulated five multivariate functional datasets according to different models proposed in [Dai & Genton 2019]. To the best of our knowledge, this is the most recent study concerned with outlier detection in multivariate functional data providing detection rates. For each synthetic dataset, n = 150 curves are generated by a common simulation model with a continuous variable discretized on a regular grid of size m = 200 in the interval [0, 1].

Among the n curves, c = 10%, referred to as the contamination level, were outliers generated by specific contamination models x c1 (t) . . . x c5 (t). Each contamination model generates one type of outlier. Testing the proposed approach and the baselines on with different contamination models enables assessing the efficiency of each mapping function for a given type of outlier.

The uncontaminated model is a bivariate Gaussian process GP(µ(t), Σ(s, t)) [START_REF] Rasmussen | [END_REF]], with a constant mean function µ(t) = 0, and a crosscovariance function C kr between the two variables indexed by k and r, as follows:

C kr (s, t) = ρ kr σ k σ r M(|s -t|; ν kr , β kr ) k, r = 1, 2 and s, t ∈ [0, 1] (II.13)
where ρ 12 is the correlation between the variables x 1 and x 2 , ρ 11 = ρ 22 is the variance of each variable, σ 1 and σ 2 are the marginal variances,

M(h; ν kr , β kr ) = 2 1-ν Γ(ν) -1 (β|h|) ν K ν (β|h|) is the Matérn class function [Matérn 2013] (K ν is a modified Bessel function [Bowman 2012]), ν kr > 0
is a smoothness parameter, and β kr > 0 is a range parameter. For this simulation, we used the same parameter setting as in [Dai & Genton 2019]:

ρ 12 = 0.6, ρ 11 = ρ 22 = 1, σ 1 = σ 2 = 1, ν 11 = 1.2, ν 22 = 0, 6, ν 12 = ν 21 = 1,
β 11 = 0.02, β 22 = 0.01, and β 12 = β 21 = 0.016. This covariance function is implemented in the R package [Schlather et al. 2015]. We summarize the uncontaminated model u(t) = (u 1 (t), u 2 (t)) as follows:

u(t) ∼ GP =   µ(t) = 0; Σ(s, t) =   C 11 C 12 C 21 C 22     (II.14)
The five contamination models are (we annotate the variables with an index c referring to "contamination"):

1. Model 1 persistent magnitude outlier: x c1 (t) = 4u(t).

2. Model 2 isolated outlier:

x c2 (t) = u(t)(1 + 11I z<t<z+0.1
), where I is the indicator returning 1 if the indexed condition is true, and 0 otherwise, and z is a uniform random variable in [0, 0.9].

3. Model 3 persistent magnitude outlier, the contamination model is different for the two variables:

x c3 (t) = (x 1,c3 (t), x 2,c3 (t)) , with x 1,c3 (t) = 1.7u 1 (t)

II.5.2 Experimental protocol

II.5.2.1 Functional-data fitting

Without loss of generality, we selected T = [0, 1] as the domain (closed interval) of t for all the data sets. We recall that we represent all the curves in the common interval T because we assume that the functional samples were generated by a random function depending on t relating to the same event in R p . For instance, when the samples are measurements of a given process depending on t, which represents time, T can be viewed as the relative temporal range of the process (i.e. from the beginning at t = 0 to the end at t = 1), and t ∈ T = [0, 1] can be interpreted as the progress rate of the process.

Choice of the basis of functions. For the ECG, PenDig and the Airbus flight test datasets, we approximated each variable of the bivariate time series by a basis consisting of B-splines of order eight (B-splines are piecewisepolynomial functions of order at least three, and are located at a given observation point t ∈ T ). Indeed, we noticed that in this dataset, the curves exhibit a smooth pattern without periodicity; hence, the B-spline basis is a suitable choice, as recommended in [Ramsay & Silverman 2006].

For the synthetic dataset, we approximated each variable of the bivariate time series by a Fourier (sine and cosine functions) basis with a fundamental period of T = 1 F = 1 (i.e. the length of T ). The Fourier basis was suitable because we noticed low-frequency periodicity (induced by the covariance function

C kr (s, t)) over T .
Application of the functional-data fitting procedure. We now provide the computational details of the functional-data fitting.

Following the recommendations in [Ramsay et al. 2009, Febrero-bande & Oviedo de la Fuente 2012], for all datasets, we selected both the penalization λ k and the basis size L ik for the variable k of sample i through a leave-one-out cross-validation procedure over a given grid search for λ k and L ik . We penalized both the first-and second-order derivatives of xik to gain smoothness in the mapping-function output. We note that for all the samples of a given variable k, we equally penalized the approximations xik by the same λ k to compute the coefficient vector α * ik . Then, by computing α * ik,λ k according to Equation (II.7), we selected the value of λ k and L ik < m i that minimize the leave-one-out cross-validation

score CV λ k (L ik ), CV λ k (L ik ) = m i j=1 x ik (t j ) -x-j ik (t j ) 2 (II.15)
where x-j ik corresponds to the approximation of x ik by L ik basis functions by omitting the pair (t j , x ik (t j )) in the functional-fitting step, as in Equation (II.5), where the penalization is λ k .

For the ECG, PenDig and AFT datasets, the grid search of λ 1 and λ 2 was fixed on logarithmic scale in [-9, -1], with a thickness of 0.1. The grid search of L ik was fixed as {35, . . . , 60}, that is, for a given L ik , the L ik Bspline functions are regularly located in T .

For the synthetic datasets, the grid search of λ 1 and λ 2 was fixed on logarithmic scale in [-9, -4], with a thickness of 0.1. The grid search of L ik was fixed in {20, . . . , 25}, that is, for a given L ik , the synthetic data were approximated by the first L ik frequencies 2πF × . . . 2πF L ik . Then, for each variable, we retained the coefficient vector associated with both the optimal regularization and basis-size parameters to recover the smooth approximation function

CHAPTER II: SHAPE-BASED OUTLIER DETECTION xi (t) = (x i1 , xi2 ).
Finally, we used the coefficient vector associated with both the optimal regularization and basis-size parameters to recover the smooth approximation functions xi (t) on a given grid and applied a mapping function to them.

II.5.2.2 Applying the mapping functions

We now explain the computational application of the mapping functions and then how their output was fed to an outlier-detection algorithm.

After computing the approximation functions xi (t), we centered and scaled each variable x ik with the empirical mean and standard deviation functions computed from the training set (see [Ramsay & Silverman 2006] for details on the computation of mean standard deviation functions). This scaling prevents the mapping functions from overweighting some variables with a wider range than others. Indeed, (i) The variables require to be scaled since the unit of the output value of the arc-length mapping function (Len map in Equation (II.9)) is intrinsically a length. Then, we applied the three mapping functions introduced in Section II.4. As the arc-length mapping is the integral function of the velocity mapping, the arc-length mapping in Eq (II.9) was computed from the minimum of T (i.e. t = 0) and was then integrated up to t for all t ∈ T . In these experiments, the integral was efficiently estimated by a Riemann sum, as in this study, all the observation points were regularly sampled in T , and therefore the sum converges to the integral. We note that if the observation points had been irregularly sampled, the integral could have not been approximated by a Riemannian sum, and numerical techniques, such as Simpson's or the trapezoidal rule, should have been used instead [Ramsay & Silverman 2006].

(ii) Regarding the velocity mapping V map in Equation (II.10), the firstorder derivative of each variable of xi (t) was computed according to Equation (II.8).

(iii) The curvature mapping (Curv map ) requires the computation of both first-and second-order derivatives. Thus, we computed them as in Equation (II.8) and combined them as in Equation (II.11).

The approximation functions recover the functional data on the entire domain T . Thus, the approximation functions can be computed on an irregular grid, and therefore the computation of the mapping functions should be carefully performed (e.g. (i) in the computation of an integral function).

For both V map and Curv map , which are based on derivative functions only, simple and efficient derivative estimation methods can be used, as mentioned in Section II.3.3.

Each mapping function returns a univariate function. Thus, applying a mapping function to all n approximation functions xi (t) results in n univariate functional-data samples. We used the resulting univariate functional data in several outlier-detection algorithms. In practice, the functions returned by a mapping function should be evaluated over a grid of observation points in T to obtain the output samples in vector form. As we selected T = [0, 1] for all datasets and the observation points are regular, the grid is a regular discretization {t 1 ...t j ...t J } of T with a thickness of 1 J (t 1 = 0 and for j > 1, t j = j J ). Hence, for the outlier-detection algorithms, the data correspond to J-dimensional numerical vectors that, in turn, correspond to univariate functional data output by a mapping function. We selected the thickness of the grid as the original size of the time series for both the synthetic and ECG datasets (ECG dataset: J = m i = m = 86, PenDig dataset: J = m = 200, AFT dataset: J = 1000, synthetic data sets: J = m = 200). An irregular grid can also be used to evaluate the approximated functions, but the computation of the mappings should be performed cautiously, as mentioned in (i) for Len map .

II.5.3 Outlier detection from the functional data output by a mapping function

We detect outliers in the functional data returned by a mapping function using a state-of-the-art outlier-detection algorithm. To this end, we selected isolation forest (iFor) [START_REF] Liu | [END_REF]] and a one-class support vector machine (OCSVM) [Schölkopf et al. 2001]. iFor is a bagging model that generates a large number of decision trees grown on random subspaces. A subspace corresponds to a subsample of features randomly selected from the full feature space (here, {1 . . . j . . . J}). Each tree isolates the data samples based on a random split value of a randomly selected feature from the subspace until all the data samples have been isolated, or all the features of the subspace have been selected. The sample outlyingness score returned by a tree is based on the path length between the root node and the terminal node of a tree. Outliers are samples that are easy to isolate and thus have short path length in the trees. The path length is normalized in [0, 1] so that if We then computed the outlyingness score of the test samples and achieved detection using the previously computed outlyingness threshold. Regarding OCSVM, we finetuned the bandwidth hyperparameter of the radial basis kernel on the training set through a 20-fold cross-validation procedure on the grid {2 -25 ...2 -5 } for the real datasets as well as the synthetic ones.

II.5.4 Result assessment

We assessed the results with the correct detection rate ρ c (i.e. number of correctly detected outlier divided by the total number of outlying curves) and the false detection rate ρ f (i.e. number of falsely detected outliers divided by the total number of inliers). In addition, as a measure of discrimination between outliers and inliers by the proposed approach, we also computed Area Under the ROC (AUC) from the labels of the test set. It is a standard performance measure in outlier detection [START_REF] Erfani | [END_REF], Liu et al. 2008] and demonstrates that the proposed method can outperform the baselines regardless of the computed outlyingness threshold.

The threshold-selection step is simple and is not part of iFor [START_REF] Liu | [END_REF] or OCSVM [Schölkopf et al. 2001], which are both unsupervised. We assume that the training data is labeled even if there are few outlier samples. In realworld applications, the user has some knowledge about the training data and can thus label inliers and some outliers.

If the training set surely has no outlier, the proposed method only requires the modification of the threshold selection rule. This modification is easy because both iFor and OCSVM are unsupervised methods and output a normalized score. Using the threshold that achieve the highest AUC, we computed the correct detection rate (i.e. number of correctly detected outlier divided by CHAPTER II: SHAPE-BASED OUTLIER DETECTION the total number of outlying curves) ρ c and ρ f (i.e. number of falsely detected outliers divided by the total number of inliers) to demonstrate the complete application of the proposed method and compare it with the baselines. There are other methods for learning an outlyingness threshold, such as using a specific decision rule involving, e.g. an empirical quantile associated with a reference distribution of the outlyingness scores [Dai & Genton 2019], or threshold selection from the mass-volume curve [Clémençon & Thomas 2017] when no outlier label is available, but this is beyond the scope of the present study, as we assume that the training set has low non-zero contamination level.

II.5.5 Baseline comparisons

We compared the proposed method with two recent outlier-detection ones based on multivariate functional depth.

The first baseline method is FUNTA, proposed by [Kuhnt & Rehage 2016] (see Section II.2). It only requires centering each variable x ik of each sample to a zero mean. As FUNTA has been demonstrated to be robust to noise and can handle curves of different size, we used it on the raw data without any functional data approximation. For the computation of the outlyingness threshold, we applied the same procedure as in the proposed method, that is, we selected the best outlyingness threshold for the training set using ROC and applied it to the test set. We used the R implementation proposed in [Rehage 2016].

The second baseline method is Dir.out proposed in [Dai & Genton 2019] (see Section II.2). We used the same parameter setting as in [Dai & Genton 2019] and did not perform any functional-data approximation. In this method, the outlyingness score is based on the robust Mahalanobis distance of the directional outlyingness vector computed on a subset of the data; in the present case, we computed it using the training data to obtain comparable results and to assess the performance measures on the test set. The tail of the distribution of the distances is approximated by an F -distribution with degrees of freedom (p+1, m-p), where p is the number of curve variables, and m is calculated through a simulation procedure (see [Dai & Genton 2019], p. 7 for details). Consequently, the outlyingness threshold is not data-driven and is computed as a quantile of probability 99, 3% of an F -distribution.

Then, we used the outlyingness threshold on the test set to asses performance. We used the R implementation provided by the authors.

II.5.6 Experimental protocol application

The performance of the proposed approach was evaluated by simulation for both the real and the synthetic data. The simulation settings for the ECG and synthetic data were as in [Dai & Genton 2019]. We proceeded as follows:

(i) We randomly generated a train/test split.

(ii) We then applied the proposed and the baseline methods. Except for Dir.out (baseline), which does not require outlyingness-threshold learning because the outlyingness score follows a known distribution (see Section II.5.5), the outlyingness threshold was learnt on the training set based on the ROC curve.

(iii) We evaluated the performance in terms of the true detection rate (ρ c ), false detection rate (ρ f ), and AU C on the test set.

For the ECG dataset (resp., PenDig datasets), steps (i) to (iii) were repeated 50 times for each case of the five contamination levels (resp., for the 10 outlier classes) (see end of Section II.5.1), and 500 times for the synthetic data for each of the five models (Section II.5.1.4). For the AFT dataset, steps (i) to (iii) were repeated 50 times.

II.5.7 Results and discussion

We report the results for the ECG dataset in Table II.1, where for each contamination level c (columns) and for each method (rows), we provide ρ c , ρ f , and AU C (sub-columns). The results for the PenDig dataset are shown 54 CHAPTER II: SHAPE-BASED OUTLIER DETECTION in Table II.2. The results for the synthetic data are reported in Table II.4. In these tables, the value in a cell is the average of a performance measure over the number of simulations. We conducted a hypothesis testing procedure in Section II.5.8 to validate the relevance of these averages. for each method (proposed in blue, and baselines in red). We notice that when the contamination level c increases, the proposed method (except for iFor(Len map ) and OCSVM(Len map )) outperforms the baselines in terms of ρ c , ρ f and AU C.

II

The results for the ECG data set (Table II.1) show that the proposed method outperforms the baselines with V map and Curv map with iFor and OCSVM. We notice that when the contamination level c increases, the proposed method (except for iFor(Len map ) and OCSVM(Len map )) outperforms the baselines in terms of ρ c , ρ f and AU C. Moreover, performance does not degrade as the contamination level increases, in contrast with that of the baselines. In terms of ρ c , F U N T A performs as well as V map and Curv map when used with both iFor and OCSVM but significantly degrades in terms of ρ f (i.e. it falsely detects outliers) for low contamination levels. Dir.out performs as well as the proposed method in terms of ρ f but degrades in terms of ρ c for high contamination levels. Hence, F U N T A performs well when the contamination level II.5. EXPERIMENTAL STUDY

55

Table II.1: Results on the ECG dataset. For all the contamination levels, the proposed method achieves the best results with Vmap and the Curvmap. For our proposed method and a given mapping function and outlier-detection algorithm, performance does not degrade when c varies, whereas for F U N T A and Dir.out it does.

c = 5% c = 10% c = 15% c = 20% c = 25% Methods ρc ρ f AU C ρc ρ f AU C ρc ρ f AU C ρc ρ f AU C ρc ρ f AU C F U N T A (
Our method outperforms state-of-the-art methods when there are few outliers.

is high, and Dir.out performs well when the contamination level is low. This shows that the outlying features captured by these mapping functions are more robust to the contamination level than those captured by the baselines.

II.5.7.2 PenDig data

From the results on the PenDig dataset in Table II.2, it can be seen that the proposed method always outperforms the baselines in terms of AU C. This implies that the baselines are not as effective in capturing shape outlying features. When the outliers are '0' digits, the results by the baselines are consistent with the results on the synthetic data when some shape outliers are simulated (Model 5 in Table II.4). This is not surprising, as Model 5

generates bivariate functional outliers with an elliptic shape in R 2 ; hence, a zero-like shape ('0'). As an AU C value close to 0.50 implies that the detector performs as efficiently as a random method, we note that the '0' outlier case is the only in which the baselines are effective. The baseline methods cannot distinguish different shape outliers with abrupt shape irregularities such as (smooth) right angles, for example, when the outlier is the '1', '4', or '5' digit.

In such cases, we obtain the best results in terms of AU C with V map . For more regular shapes, such as '3', '6', '8', and '9', the best results are achieved by Curv map . We conclude similarly for the ρ c results. However, baselines are better in terms of ρ f which shows that our method can confuse inliers with outliers and lack robustness of this dataset. Outlier detection results for the PenDig dataset when each of the 10 classes ('0'...'9') is considered an outlier (columns), and the nine other classes inliers. For each case of outlier class and each performance measure (sub-columns), we marked the best results in bold. We note that our method achieve the best results in terms of AU C, showing that it better discriminates outliers. Outlier detection results for the AFT dataset with five contamination levels c. We report the results on the AFT dataset in Table II.3 where one can see that our proposal outperforms the baselines with V map and Len map in terms of ρ c , ρ f and AU C. These results thus suggests that the outlying multivariate time series of the AFT dataset are underlied by a process whose duration and velocity are both abnormal. Also, since the Curv map does not enhance the discrimination between outliers and inliers, we can also say that they are underlied by a path whose change in direction (in R p=4 ) are similar.

II.5.7.3 AFT data

II.5.7.4 Synthetic data

For all types of outliers (generated by a distinct contamination model), iFor(Len map ) achieves similar results to the baselines. For shape outliers (Model 5),the baselines are outperformed.

For isolated outliers (Table II.4, Model 2 and Model 4 columns), the results on the synthetic datasets show that the baseline methods perform as well as Len map , and V map with iFor. Moreover, since in this case the outlyingness is a short peak, the velocity quickly changes when the peak occurs; thus, the V map function is an appropriate candidate for detecting isolated outliers. Curv map Table II.4: Results on the five synthetic datasets

CHAPTER II: SHAPE-BASED OUTLIER DETECTION

shows poor performance for the two models. This implies that it is ineffective in detecting isolated outliers. Indeed, the contamination models (Model 2 and Model 4, Section II.5.1.4) generate stationary functional data (constant mean and only lag-dependent covariance) except in the part of T where the outlyingness (low short peak) occurs. Thus, considering the second-order variations (second-order derivatives in Equation (II.11)) is irrelevant and leads to high ρ f values (ρ f columns and Curv map rows).

For persistent magnitude outliers (Table II.4, Model 1 and Model 3 columns),

Dir.out and F U N T A yield the best results in terms of both ρ c and ρ f . We obtain similar results for Model 1 with V map , and Len map with iFor.

Nevertheless, V map is not as efficient for Model 3 as for Model 1. Indeed, Model 1 has high contamination (high short peak), resulting in high velocity mapping values, and we recall that velocity and curvature relate to local variations of the curves. Consequently, as here the magnitude outlyingness is a global feature, Len map is better than V map and Curv map . This indicates that for detecting persistent magnitude outliers, the proposed approach is more reliable with Len map than Curv map and V map .

For persistent shape outliers (Table II.4, Model 5 column), the proposed method outperforms the baselines with iFor on Len map . Furthermore, V map yields results similar to those of Dir.out in terms of ρ c and AU C. Table II.4

shows that the state-of-the-art F U N T A totally fails to capture shape outlyingness because it is based on the intersection angles between the samples and is computed for each variable separately. Thus, it fails to consider the correlation between them (as explained in Section II.1).

For Models 1-5, among the proposed mapping functions, Len map achieves the best results and appears to be suitable for outlier detection if the variables are weakly correlated, whereas V map and Curv map are preferable if the correlation between the variables is strong.

As V map and Len map achieve satisfactory results, the geometric characterization (velocity and length) of the samples provides a different shape features in the outlier detection. We note that functional-data approxima-tion affects the geometric characterization. Indeed, functional approximation enables smoothing out a curve and properly extracting derivative-based features because the induced smoothing renders the samples differentiable (see Section II.3.3); this is not a required property for the baselines Dir.out and F U N T A. Here, we carefully monitor the functional-approximation step using leave-one-out cross-validation (Equation 14). Thus, in contrast with the approximation step, the outlier-detection step depends greatly on the mapping-function computation.

Overall, considering the results on synthetic datasets, we recommend using

Len map in the case of (potential) persistent magnitude or shape outliers. In practice, Len map does not directly indicate whether a sample is a shape or magnitude outlier. However, as shape and magnitude are quite distinctive outlyingness classes, the class of such an outlier can be known a posteriori by visual inspection or by setting a magnitude threshold with respect to the magnitude of the outliers detected. If the outliers are suspected to be isolated, we recommend using V map and Curv map , as both mapping functions extract local curve features in R p . In the case of a low contamination level, both OCSVM and iFor are suitable (even though on the ECG data, OCSVM is better for small c), whereas for high contamination levels, iFor is better.

We demonstrated that each mapping function can detect multiple classes of outliers. The results are consistent on both synthetic and real data. However, knowing the class of a outlier detected by a given mapping function is not easy, and this issue will be addressed in future work.

II.5.8 Statistical assessment of the results

We followed the hypothesis-testing procedure recommended in [START_REF] Demsar | Statistical Comparisons of Classifiers over Multiple Data Sets[END_REF], García et al. 2010] (also used in [START_REF] Erfani | [END_REF]) to assess the statistical significance of the results introduced above. Demsar provided an evaluation protocol for a more general assessment of the difference between several classifiers used on multiple benchmark datasets. The protocol consists of two steps: First, a global significance test is conducted to determine whether

CHAPTER II: SHAPE-BASED OUTLIER DETECTION

there is a difference among the evaluated methods. Second, if this is the case, the methods are pairwise compared to evaluate the gain of one over another.

We applied Demsar's protocol because the present detection task reduces to a two-class classification in the evaluation step (outliers/inliers). We applied the protocol for the three performance measures ρ c , ρ f , and AU C separately.

As described in [START_REF] Demsar | Statistical Comparisons of Classifiers over Multiple Data Sets[END_REF], García et al. 2010], there are several ways of conducting the tests and we primarily applied it as in [START_REF] Erfani | [END_REF].

We applied the protocol as follows:

(i) First, the Friedman test [Sheskin 2003] was applied to detect the global statistical significance for each of the three performance measures among all the methods on all the datasets. The Friedman test can be viewed as the non-parametric version of ANOVA (where, here, a group refers to a method, and the samples in the group refer to the performance of the method on the datasets), as it is based on the ranks and thus does not make the Gaussian assumption for the performance measures for each method [START_REF] Demsar | Statistical Comparisons of Classifiers over Multiple Data Sets[END_REF]]. We conducted the Friedman test with the Iman-Davenport correction [Sheskin 2003], as recommended in [START_REF] Demsar | Statistical Comparisons of Classifiers over Multiple Data Sets[END_REF]], to handle the well-known family-wise error rate, which can bias the p-value in a multiple-hypothesis test.

We recall that in the present context, the family-wise error rate refers to the probability of erroneously asserting that one method is more reliable for detecting outliers than some of the others.

(ii) Second, if statistical significance was detected by the Friedman test, we performed a post-hoc test to determine which methods are different.

More precisely, the post-hoc test is based on the p-values returned by a pairwise-comparison test applied to all pairwise comparisons of the methods. A nonparametric test can be selected for the pairwise comparisons (owing to the absence of the Gaussian assumption), such as the post-hoc Friedman's aligned ranked test [García et al. 2010]. As the all-pairwise-comparisons test is a special case of a multiple-hypothesis test, it also suffers from the family error rate and requires a correction procedure. Thus, we used the Finner correction as recommended in [García et al. 2010].

We separately applied this evaluation protocol to the three performance measures for the five contamination levels of the ECG data, the 10 outlier classes of the PenDig dataset, and the synthetic data to compare the methods on two distinct types of data and to demonstrate the benefit of the proposed approach on real data. Moreover, this enables assessing the difference of the methods in a given context (i.e. when the outlier class is known). For all datasets, we used a significance level of 0.1, as in [START_REF] Erfani | [END_REF].

We report the average ranking (vertical axis) of all methods (horizontal axis) applied to the ECG and PenDig datasets (resp., synthetic data) for each performance measure (colors) in a vertical-bar plot in ECG data. The Friedman test rejects the null hypothesis of equivalence of the methods for the three performance measures at a significance level of 0.1.

The p-values are 3.0 × 10 -10 for the correction detection rate ρ c , 3.0 × 10 -10 is consistency with respect to the ECG data except for the false detection rate ρ f . Indeed, both V map and Curv map outperform the baselines in terms of ρ c and AU C (Tables II.5 andII.7). Moreover, among the three mapping functions, Len map yields the worst results and is not different from Dir.out.

However, there is an inconsistency ranking regarding ρ f in the PenDig data with respect to the ECG data (Figure II.6 and Table II.6). Indeed, as the proposed method is not ranked first in terms of the false detection rate, it may be claimed that it recognizes the outliers but tends to be excessively severe.

We note that this conclusion regarding the correct and false detection rates is drawn according to the adopted outlyingness thresholding rule, which can be modified, as discussed at the end of Section II.5.3.

From the global ranking (Figure II.6) and the pairwise comparison tests, it may be concluded that the proposed method outperforms the baselines on both the ECG and PenDig datasets at the significance level 0.1.

Synthetic data.

Regarding the synthetic data, the Friedman test rejects the null hypothesis of equivalence of the methods for the three performances measures at a significance level of 0.1. The p-value is 2.4 × 10 -10 for the correct detection rate, 2.4 × 10 -10 for the false detection rate, and 1.0 × 10 -6 for AU C. As these p-values are below 0.1, we conducted a post-hoc test to compare the methods pairwise and assess the gain of one over another. The p-values of each pairwise comparison test is given in Tables II.8, II.9, and II.10 for ρ c , ρ f , and AU C, respectively. We notice that Dir.out is significantly equivalent to iFor(Len map ), OCSVM(Len map ), F U N T A, and iFor(V map), and these methods are ranked first, second, and third on average, respectively (Figure II.7). Thus, on the synthetic dataset, the baseline methods are slightly better than the proposed method; however, based on the pairwise comparison tests, the best methods (iF or(Len map ) and OCSV M (Len map )) II.5. EXPERIMENTAL STUDY 69 are statistically equivalent. As discussed in the two previous paragraphs, the proposed method is superior on real datasets. Moreover, in the iFor rows and OCSVM columns, it can be seen that there is a pairwise equivalence between iFor and OCSVM for (Len map ) and (V map ), that is, these two outlier-detection algorithms are empirically consistent for a given mapping function. Therefore, we have equivalent methods to achieve state-of-the-art results (which cannot be improved, except for MODEL 5) for the synthetic data.

Overall assessment. Tables II.5, II.6, and II.7 (in the iFor rows and OCSVM columns) show the pairwise consistency between the iFor and OCSVM algorithms for each mapping function. The same holds for the synthetic data. Thus, for a given dataset and mapping function, iFor and OCSVM achieve statistically the same performance results. This implies that the detection performance relies more on the outlying features provided by the mapping function than on the capacity of the outlier-detection algorithm to discover outlying features itself.

Table II.8: Statistical significance of the pairwise comparisons for the correct de-tection rate ρ c on the synthetic datasets. For ρ c and AU C, the ranking is given in decreasing order (i.e. for high ρ c and AU C values, the rank tends to 1); for ρ f , the ranking is given in increasing order (i.e. for low ρ f values, the rank tends to 1). The y-axis represents the average ranking over the five models, and the integers on the top of the bars represent the final ranking. If there are ties, we take the average ranking. (1 is the best, 8 is the worst) for ρ c , ρ f , and AU C based on the post-hoc Friedman's aligned rank test. For ρ c and AU C, the ranking is given in decreasing order (i.e. for high ρ c and AU C values, the rank is close to 1), and for ρ f , the ranking is given in increasing order (i.e. for low ρ f values, the rank is close to 1). The y-axis represents the average ranking over the five models, and the integers on the top of the bars represent the final ranking. If there are ties, we take the average ranking. The baseline methods are slightly better than the proposed method, but the best results by the proposed method (iF or(Len map ) and OCSV M (Len map )) are statistically equivalent to those by the baseline, as showed by the pairwise comparison tests in Tables 8,9, and 10.

II.6 Conclusion

In this chapter, we proposed a method to improve the detection of different types of outliers in multivariate functional data, based on curve shape.

We assumed that the original discrete curves can be well approximated by finite functional basis expansions, where the basis is specified. Based on the smooth reconstruction provided by the fitted basis expansion, we used the arc-length, velocity, and curvature mapping functions to capture latent shape features. Then, we detected the outliers from the mapped curves using outlier-detection algorithms.

Through an experimental study on real and synthetic datasets, we showed that the proposed approach outperforms multivariate functional depth baselines on real data and can perform similarly on synthetic data (except for persistent shape outliers, where the proposed method performs better). We demonstrated that, compared with the baselines, the proposed approach is robust to the variation of the contamination level.

We did not assume any weighting of the curve variables in the mapping functions; this is left as future work. This weighting could be user-driven, as proposed for functional depth in [START_REF] Claeskens | [END_REF], or data-driven. It is conceivable that this can enhance outlier detection in the presence of nonoutlying curve variables (when p increases). Another possible improvement would be to combine mapping functions in the same detector so that multiple outlier classes may be detected in the same dataset.

II.7 Publications

We shall mention that the work of Chapter II has been published in two venues. We first published the short paper [Lejeune et al. 2020b] in the international Conference Extended Data Base Technology. In this paper, we used the curvature mapping function on the ECG data set and showed empirically that our method is more robust than some baselines in terms of contamina-tion rate. From repeated trials, we also showed that our method provides predictions with lower uncertainties than the baselines. We secondly published the journal paper [Lejeune et al. 2020a] in the Knowledge-based Systems journal. In this paper, we extended our work in [Lejeune et al. 2020b] by proposing the arc-length and the velocity mapping functions. We validated our method on synthetic datasets and on the ECG and PenDig real datasets.

Chapter III

Data-driven discovery of systems of ordinary differential equations with nonconvex multi-task learning

In this chapter, we tackle the problem of analytical discovery of systems of ordinary differential equations (SODE) from a multivariate time series.

Once known, the analytical form of such model provides a quantitative representation of the dynamics underlying an observed phenomenon. The problem is encompassed by the recent framework known as the sparse identification of nonlinear dynamics [START_REF] Brunton | [END_REF]. In state of the art, the problem is formalized as a linear regression whose weights are learned with the sparse convex LASSO penalty. However, this penalty does not take into account the coupling between the equations of the SODE and, due to its convexity, gives a biased estimate of the weights. To address these two limitations, we re-cast the problem as a multi-task learning one involving a nonconvex penalty. Numerical experiments, on known SODEs,

show that both the multi-task and the nonconvexity features of our method outperform state of the art. We also apply our algorithm on an Airbus flight test dataset.

Abstract.

III.1 Introduction

Governing equations are mathematical models, such as partial differential equations (e.g. Navier-Stokes equation for fluid dynamic modeling) or systems of ordinary differential equations (SODE), widely used in sciences and engineering to model dynamic systems [Butcher 2016]. A governing equation models the dependency relationship between several parameters, named state-variables, (e.g. velocity, chemical concentration, temperature, etc.), in a dynamic system.

A SODE is a system of equations where each one models the dependency relationship between several scalar-dependent state-variables and one state- Governing equations are traditionally derived from principled rules themselves formalized from general empirical observations under certain hypotheses. For instance, the SODE of the damped harmonic oscillator is derived from the Newton's second law under the constant-mass hypothesis [Greiner 2006].

However, there remains some complex systems (e.g. in mechanical engineering, fluid mechanics) whose dynamic behavior is poorly understood and quite hard to be modeled within a governing equation derived from existing principled rules. Hence, the SODE underlying such systems is unknown. Accessing to the model, like a SODE, that governs an unknown dynamic is a challenging task of scientific interest to improve the understanding of a physical system as well as of practical interest to get a simulation model e.g. for prototype design [START_REF] Brunton | [END_REF][START_REF] Schaeffer | [END_REF]. In this chapter, we address such a challenge and propose a data-driven method to solve it based on multi-task learning (MTL). In Thanks to the development of sensor technology, many data can be sampled from a dynamic system. Such an amount of data gives new opportunities to extract important knowledge on the underlying physical behavior of a dynamic system. Consequently, there has been a growing interest over the last years in the development of data-driven methods for the discovery of governing equations [START_REF] Brunton | [END_REF], Long et al. 2018, Bhat & Rawat 2019, Schaeffer 2017, Zhang & Schaeffer 2019].

State-of-the-art methods for the discovery of SODE are based on the matrix-extended LASSO (the • 1,1 norm) sparsity-promoting learning methods [START_REF] Brunton | [END_REF][START_REF] Schaeffer | [END_REF], Rudy et al. 2019, Tishbirani 1996]. In Figure III.1, we illustrate and generalize the core idea of these methods: such kind of methods takes as input samples of state-variables time derivative and a dictionary of resultants of (arbitrary) nonlinear transformations of state-variables samples. Then the algorithm outputs the relevant elements of the dictionary that best model the relationship between the (samples of) state-variables and their associated first order derivatives.

Actually, learning in this way, i.e. with the • 1,1 norm penalty, reduces to single-task learning, where a task refers to the discovery of a single equation.

However, (i) single-task learning does not accommodate the relatedness between the equations of the SODE [START_REF] Obozinski | [END_REF]. We give an example of such relatedness in Figure III.1 bottom-left where the occurrence of x 3 1 and x 3 2 , within both equations, makes the SODE coupled. Consequently, issue (i) results in an over-complete (i.e. not sparse enough) or under-complete (too sparse) set of selected dictionary elements to identify the SODE.

Moreover, since the LASSO penalty is convex, (ii) the learned coefficients associated to each element of the dictionary result biased [Zhang 2010, Fan & Li 2001].

To remedy both (i) and (ii), we propose to cast the problem of discovering a SODE as a MTL problem [START_REF] Argyriou | [END_REF], Caruana 1997, Obozinski et al. 2010] involving nonconvex regularizer that promotes sparsity, relatedness between the equations of the SODE and encourages unbiasedness.

Our contribution is to learn a SODE from data samples of a dynamic system with sparsity-promoting algorithm in the same spirit as [START_REF] Brunton | [END_REF][START_REF] Schaeffer | [END_REF], Rudy et al. 2019], while preventing both inconsistent sparsity and bias issues. Actually, learning a SODE can be seen as a problem made up of p regression tasks. The k-th learning task consists to learn the k-th component of the SODE from a large dictionary of linear and nonlinear functions built from the raw data i.e. the samples of the state-variables of the dynamic system under study. Since, by definition of a SODE, there are as many tasks as state-variables, the regression problem results in a MTL problem. Contrary to single-task learning, MTL takes benefit from the correlation between the tasks to learn a consistent set of predictors across tasks in order to improve predictive performances of each of of them [START_REF] Argyriou | [END_REF], Obozinski et al. 2010]. In other words, when applied to the discovery of a SODE, MTL amounts to learn a sparse set of elements from a dictionary of candidate functions for each equation by considering the coupling feature of the SODE.

We sum up our two-fold contribution:

• We cast the discovery of SODE as a nonconvex MTL problem.

We formalize the learning as an optimization problem involving a matrix-structured, sparse and nonconvex regularizer to account both for task relatedness and unbiasedness in the learned SODE. To perform the learning, we instantiate an efficient generic algorithm [Gong et al. 2013].

• Through experiments on benchmark of reference SODEs, we show that learning with a nonconvex multi-task penalty enables a better recovery of the underlying equations than learning with a convex single-task penalty.

combinatorial approach based on genetic-programming to select the parsimonious model, among a large set of candidate models, that best recovers the data. As mentioned in [START_REF] Brunton | [END_REF], genetic-programming meth-

III.2.1 Building block for sparse learning of a SODE

State-of-the-art methods for the discovery of SODE, [START_REF] Brunton | [END_REF], Rudy et al. 2019], assume that each one of ẋ1 , . . . , ẋp in Equation (III.1) are independent targets which can be predicted by a sparse combination of elements comprised in a dictionary of candidate functions. See Figure III.1 for an example of a two-dimensional SODE and where the dictionary is denoted Θ Xn . The dictionary is first specified by the user by building linear as well as nonlinear candidate functions (e.g. polynomials) from the noisy samples X n of x. This dictionary reflects the prior knowledge on the observed phenomenon and is possibly over-complete. In [START_REF] Brunton | [END_REF][START_REF] Schaeffer | [END_REF], Rudy et al. 2019], f is assumed to be linear with respect to (w.r.t) the elements of the dictionary (and not w.r.t t). The linear assumption on f w.r.t x 1 , . . . , x p makes it easy to learn and to interpret. Then to learn the SODE, a sparsity-promoting algorithm, e.g. LASSO [Tishbirani 1996], elastic-net [Zou & Hastie 2005], learns f 1 , • • • , f p separately by both selecting the relevant elements (i.e. a small set of candidate functions) of the dictionary and estimating their associated coefficient in the linear model.

III.2.2 Discovery of a SODE by sparse linear regression

The first step of the approach introduced above can be described as follows [START_REF] Brunton | [END_REF]. Starting from n noisy samples of a p-dimensional state-variable comprised in X n and the associated time-derivatives samples Ẋn = ẋ•1 , ẋ•2 , . . . , ẋ•p (which can be computed numerically if they were not sampled), one first builds an arbitrary dictionary of m candidate func-

tions Θ Xn = x •1 , x •2 , . . . , x 2 •1 , x 2 •2 , . . . , cos x •1 , . . . ∈ R n×m .
Then from the linear assumption on f , i.e. Ẋn = f X n = Θ Xn β where β ∈ R m×p is a matrix wherein the q-th column refers to the coefficient-vector associated to the candidate functions of the q-th SODE component (see Equation (III.1)), one can find a sparse β * by minimizing a loss (data fidelity term) plus a sparsity-promoting term:

β * = arg min β∈R m×p ( Ẋn , Θ Xn β) + λR(β) (III.2)
λ > 0 is the sparsity amount hyperparameter. We illustrate our approach in Figure III.1 with a two-dimensional SODE.

To learn the SODE, [START_REF] Brunton | [END_REF], Rudy et al. 2019, Schaeffer 2017] instantiate Problem (III.2) by choosing R as the 1,1 norm i.e. β * is the solution of the following problem:

β * = arg min β∈R m×p 1 2 Ẋn -Θ Xn β 2 F + λ β 1,1 (III.3) which is a special case of Problem (III.2) where • 2 F = • 2 2,2 = p j n i (•) 2 is the squared Frobenius norm that serves in the loss (•, •) and • 1,1 is the sparsity-promoting regularizer R(•). Note that, since β = β •1 , • • • , β •p and
• 1,1 acts independently on each entry of β, solving Problem (III.3) reduces to p independent LASSO [Tishbirani 1996] sub-problems where each one is a learning task that consists to estimate, for k

= 1 . . . p, β •k from ẋ•k , x •k
with the 1 norm weighted by λ. In this way, learning the SODE is performed in a single-task learning manner. [Schaeffer 2017] formalized the discovery of nonlinear dynamics in the case of partial differential equations and thus the dictionary used in his algorithm is different to our Θ Xn as it involves partial derivatives. It turns out that the learning criterion that the author used was formalized similarly to Problem (III.3) for SODE. Since Problem (III.3) is convex in β, Schaeffer solved it with the Douglas-Rachford algorithm [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF]] (Algorithm III.1) which is a proximal algorithm [Parikh & Boyd 2013]. One can see that the main steps of Algorithm III.1 are in lines 4 and 5. These steps consist to iterate between the shrinkage of W (sparsity promotion) followed by the regularized inversion of Z (loss minimization). Step 6 maintains sparsity in W . The algorithm can be applied in low as well as in high-dimension settings, m n, i.e. when there are more candidate functions than samples in Θ Xn . The notion of proximal operator is at the core of this learning algorithm as well as in our proposal for MTL, so we recall its formal definition.

III.3 Learning algorithm

Algorithm III.1: Douglas-Rachford algorithm for Problem (III.3) Input: data samples Ẋn , X n , sparsity amount λ > 0, initial guess β 0 , 0 < µ < 2 1: build Θ Xn arbitrarily from X n 2: W ← β 0 , β ← W 3: while β has not converged do 4: thresholding operator [Tishbirani 1996].

Z ← 2prox λ • 1,1 (W ) -W 5: W ← µ(I + λΘ Xn Θ Xn ) -1 (Z -λΘ Xn Ẋn ) -Z +(1 -µ 2 )W 6: β ← prox λ • 1,1 (
Roughly speaking, in the learning algorithm of a linear model, the proximal operator of a convex non-differentiable function (e.g. • 1,1 ) serves as a shrinkage operator which assigns zero to coefficients in β that do not decrease the learning criterion enough. Actually, the exact shrinking to zero, of the coefficients associated to irrelevant candidate functions for each task, is due to the non-differentiability of the regularizer. Note that • 1,1 is used as a sparsity-promoting regularizer in Problem (III.3), is separable and thus does not consider any matrix structure. Thus, for MTL, the separability of regularizer is not desirable. Indeed, one can permute any element within the coefficient matrix β to learn, the resulting 1,1 norm remains unchanged.

Consequently, the regularizer acts as if β were a vector in R mp .

III.3.1 Shortcomings

Considering a simple vector structure on β through a fully separable norm regularizer, rather than a matrix-structured one, amounts to omit task relatedness i.e. correlations between columns β

•1 , • • • , β •p . Furthermore,
the convexity of a norm in R m×p induces a bias in the learned coefficients [Fan & Li 2001[START_REF] Boyd | [END_REF]

III.3.2 Building-block of MTL of linear regressions

Since we address the discovery of a SODE with multiple linear regressions, herbelow we describe the building block of MTL in this framework. The core idea of MTL consists to regularize all the tasks jointly during the learning.

Here a task refers to the discovery of a single equation of the SODE. We discuss stat of the art sparsity-promoting regularizers that can deal with sparsity across tasks.

MTL consists in learning p functions f 1 , • • • , f p jointly by assuming that they are close to each other in some similarity metric [START_REF] Argyriou | [END_REF] and share a common set of features. For the k-th task, one is given a data set {y ik ; z ik } i≤n k of n k samples with m features. Therefore, the p regression coefficient vectors can be represented in a matrix

β = β 1 , • • • , β p ∈ R m×p .
The task similarity is reflected within the learning criterion by a regularizer applied on this matrix. For p linear regressions, MTL can be formulated as computing β * as:

β * = arg min β=[β 1 ,••• ,β p ] p k=1 n k i=1 1 2n k (z ik ; y ik β k ) 2 + λR(β) (III.4)
where R is the regularizer that may take into account for task relatedness and λ > 0 is the regularization amount. Note that when n k = n and R is • 1,1 , Problem (III.4) reduces to Problem (III.3). Thus by choosing a regularizer, more appropriate than the 1,1 , to considers relatedness between tasks, the discovery of a SODE can be reformulated as a MTL problem. Hence solving for Problem (III.3) corresponds to learn f 1 , • • • , f p independently.

III.3.2.1 Considering task relatedness

To account for task relatedness in the learning and not to perform p independent single-task learning, the regularizer has to consider a matrix structure on the coefficients. For instance, solving Problem (III.4) with R chosen as R 2,1 , i.e. β 2,1 = m i β i• 2 matrix norm (the so-called group-lasso [Yuan & Lin 2006]) makes β row-sparse i.e. some of its rows are nonzero and all the others exactly equal zero. Such a regularizer enforces all the components of the SODE to have the same nonzero coefficients and thus to share the same candidate functions specified in Θ Xn [START_REF] Obozinski | [END_REF].

III.3.2.2 Considering task specific elements

Rather than only considering specific candidate functions with the entrywise • 1,1 norm, or oppositely, only common candidate functions across the components of the SODE with the • 2,1 norm, it is more realistic to consider both. Such a compromise can be achieved by taking R as a convex combination of these two norms

i.e. R 2,1 + 1,1 = α • 2,1 + (1 -α) • 1,1 with α ∈ [0, 1].
For single-task linear regression, it is formulated as an extension of the group-lasso when groups of features are known and some of the features do not act at the group level [START_REF] Simon | [END_REF]. When α > 0.5, the learning criterion attaches more importance on common candidate functions across tasks than task specific ones, and conversely when α < 0.5. In practice, α is set around 0.9 to allow few task specific components but it can be fine-tuned by cross-validation.

III.3.2.3 Shortcomings

Despite being able to select relevant candidate functions, the norms R 2,1 and R 2,1 + 1,1 induce a bias in the learned coefficients within β. Similarly to • 1,1 , this bias is induced by the convexity of such regularizers through their associated proximal operator which are available analytically (no optimization step is needed to evaluate them) and given in [START_REF] Simon | [END_REF].

We consider these regularizers as baselines in our numerical experiments in Section III.5 and show that the nonconvexity is important to get an accurate estimate of the SODE.

In Section III.2.1 and Section III.3.2, task relatedness and nonconvexity are showed to be two important weaknesses not addressed by the 1,1 regularizer used in state-of-the-art algorithms. Our contribution leverages both the nonconvexity and task relatedness within a single regularizer to improve the discovery of a SODE. One can note that such an extension results, by construction, in a regularizer having the separability property. Thus, the proximal operator of R SCAD λ,θ is obtained by evaluating prox r SCAD λ,θ entrywisely. However, since a separable regularizer cannot take into account relatedness between tasks, like the 1,1 norm, the benefit gained by learning with R SCAD λ,θ over the latter is only unbiasedness. Indeed, the first summand in R SCAD λ,θ acts both at the row and column levels of β. In other words, it acts for the coefficient β ij associated to i-th candidate function (i-th column of Θ Xn ) of the j-th task, and therefore independently across the p equations of the SODE. We propose to "un-separate" R SCAD λ,θ to both leverage from non-separability as well as nonconvexity of the SCAD penalty.

III.4.2 Nonconvex non-separable regularizer

We propose to learn

β = β 1• , • • • , β m• in Ẋn = Θ Xn β,
by solving Problem III.4 with our proposed regularizer R SCAD-1 . To make the proposed regularizer non-separable, such that it accounts for task relatedness as well as unbiasedness, the key point is to replace the second summation of R SCAD λ,θ by the evaluation of r SCAD λ,θ of the 1 norm of β i• for every i:

R SCAD-1 λ,θ (β) = m i=1 r SCAD λ,θ β i• 1 (III.7)
In this manner, since the regularizer acts for the i-th candidate function onto the coefficient vector β i• ∈ R p through the SCAD penalty of the 1 norm, γ ← γ min 12: end while Output: β loop (lines 4-12) consists in performing a loss gradient descent step (line 9) followed by a shrinkage operator step (which here learns the coefficients of the underlying SODE). In practice, the outer loop (line 4) is stopped if the 2,2 norm of the relative change of β between two consecutive iterations is less than a low value like 10 -5 or the total number of iterations is greater than 5000. The inner loop (lines 5-8) consists to perform a line search to compute a gradient step size γ so that it ensures a sufficient decrease of the loss along the negative gradient direction. Here we used the backtracking line search criterion with parameter value 0.8 (line 7), which is commonly used, and serves as the 'slow-rate' of the line search (must lie in ]0, 1[) [START_REF] Boyd | [END_REF]. Note that learning with R 1,1 in Algorithm III.2 approximately reduces to learning with Algorithm III.1. In the next section we instantiate Algorithm III.2 with different regularizers and compare their ability to recover a SODE from noisy data. minimizing the Bayesian information criterion (no need for train/test splits with cross-validation) 1 in the logarithmic grid {λ 0 = 10 3 , . . . , λ N -1 = 10 -2 } containing N = 10 3 values. As this requires one learning for each λ, we performed warm-start to estimate the models along λ [Friedman et al. 2010].

Warm-start consists to first estimate βλ 0 with the largest sparsity amount λ 0 , such that βλ 0 = 0, and then to (sequentially) estimate βλ i by initialising GISTA with β 0 = βλ i-1 .

III.5.2 Comparison with baseline regularizers

We compared our approach by learning with: R 1,1 , R 2,1 , R 2,1 + 1,1 (α = 0.9), for which the proximal operators are available analytically [Chierchia et al. 16 ]. For R 1,1 and R 2,1 , we used Algorithm III.1(µ = 1.5) and the multi-task-LASSO (MTLa) from Scikit-Learn, respectively. We used

III.5.3 Results

III.5.3.1 Synthetic SODEs

For each SODE, we repeated the experiment ten times. We report the average and standard deviation of β , T and M IS in Table III.1. We show the simulated SODEs, with their closed form, from their best estimate over the ten trials in Fig III .6. The results show that on the DOC, LV and LAT datasets, learning with R SCAD-1 outperforms the convex baselines (Table III.1 R 1,1 , R 2,1 , R 2,1 + 1,1 rows) both in terms of SODE identification (smallest M IS ) as 1 However, without training set, by abuse of language we keep on using "learn" that refers to compute an estimate of β R SCAD-1 (our) 10 -4 ± 10 -4 3.5 ± 3.5 8.10 -2 ± 10 -3 2.10 -3 ± 1.10 -3 7.8 ± 2.1 1.10 -4 ± 10 -3 4.10 -5 ± 3.10 -5 0.0 ± 0.0 13 ± 10 -3 β measures unbiasedness. M IS is the misidentification error of the SODE. T is the error between the solution of a learnt SODE and the ground truth. For each SODE, the two smallest errors are in bold.

well as unbiasedness (smallest β ). Fig III.6 (LAT simulations) shows that even if the SODE is well identified (small M IS ) with a convex regularizer, the bias of the coefficients leads to a degradation of the forecasting performances of the SODE. Moreover, despite the a priori inability of R SCAD to consider the coupling within an SODE, the results nevertheless (Table III.1 penultimate row) show that it is good for variable selection (except for LV whose dictionary Θ Xn involves more correlations than for DOC and LAT).

Such a result can be attributed to the (semi) concavity of R SCAD that leads to a penalty with stronger sparsity (i.e. nearer to the 0 oracle) than convex regularizers. As a global result, our experiments show that leveraging both from MTL and nonconvexity improves the discovery of the SODE in closed form. We report the error T , for each SODE discovered by the five regularizers, in Table III.2. The results show that learning with R SCAD-1 outperforms the baselines. We note that the second best result is achieved by learning with R 2,1 , meaning that considering the coupling within the SODE, regardless of the bias induced by convexity, is a first benefit of our MTL-based approach.

III.5.3.2 Airbus flight test data

Thus, the results are consistent with the ones of the synthetic SODEs, which

III.6 Conclusion

We re-casted the learning of an unknown SODE from noisy data as a MTL problem. We proposed a nonconvex regularizer that better (i) accommodates the coupling of the equations of the underlying SODE than convex regularizers and (ii) results in unbiased coefficients. To learn the SODE with our nonconvex regularizer, we instantiate a generic algorithm from the literature.

Numerical experiments on real and synthetic datasets confirm that both the MTL feature and the nonconvexity of our proposal outperform learning with state of the arts regularizers.

General conclusion and perspectives

In this thesis, we tackled two problems in multivariate time series analysis: the detection of outliers and the discovery of SODE. To address both problems, we represented a time series as function over time.

In Chapter II, we used the functional data representation as a pre-processing step. Then, we proposed to extract shape features with the arc-length, velocity and curvature mapping functions which are well established in differential geometry. We used the obtained shape-based data representation as input of an outlier detection algorithm. Through numerical experiments on synthetic and real datasets, we showed the improvement of our approach, in terms of detection and robustness to the contamination level, w.r.t to state of the art.

As future work, a possible improvement is to combine mapping functions so that multiple classes of outliers in the same dataset might be detected. For In Chapter III, we represented the time-dependent function underlying a noisy time series as the solution of an unknown SODE. Accessing to the analytical form of the SODE is of a great interest since it can give an explicit, and functional, understanding on a particular dynamic system. This inverse problem is formulated as a sparse regression problem. However, state of the art algorithms are single-task based and thus omit the coupling feature within the SODE. Furthermore, they involve the convex 1 sparse penalty that results in biased estimate of the regression weights. To remedy these two limitations, we re-casted the problem with MTL involving a nonconvex penalty. Our numerical experiments showed the improvement of our algorithm w.r.t to state of the art.

A limitation of our approach, also not addressed by state of the art, is its restriction to discover the SODE underlying a single multivariate time series.

Therefore, for a dataset comprising multiple multivariate time series, this is required to learn each SODE independently. In this way, as the size of the dataset increases, the learning can be costly and does not consider the potential similarities between the time series. A future work, would be to jointly discover multiple SODEs with a constraint of similarity between the regression weights across all the SODEs.

•

  R d : the set of real valued d-dimensional column vectors • R n×d : the set of real valued n × d matrices • a: an element of R • a = a 1 , . . . , a d : an element of R d • A = a 1 , . . . , a n : an element of R n×d • a k,• : the k-th row of A • a •,l : the l-th column of A • a, b : the standard (Euclidean) inner product between a and b

•

  x(t): the value of the function x at t. If x(t) ∈ R d=1 , then x(t) = x(t) • x: an approximation of the function x • N (µ, Σ): the Gaussian distribution with mean µ and (positive-LIST OF TABLES definite) covariance Σ • U(a, b): the uniform distribution in the a × b rectangle with a i > b i Chapter IGeneral introductionIn this chapter, we first present the industrial context that has motivated our research exposed in this thesis. For that, we define what a flight test is, and introduce the issues involved by the analysis of time series generated by flight tests, i.e. the detection and understanding of abnormal dynamic behaviors. Secondly, we both define the notions of time series and outliers, which we encompass in the functional data framework. Thirdly, we sum up our two contributions: the first one addresses the detection of outliers in functional data, and the second one tackles the discovery of systems of differential equations from a time series.In civil aviation, to be allowed to operate, an aircraft must be certified by a government certifying agency (e.g. European Aviation Safety Agency) independent from the manufacturer (e.g. Airbus). To achieve the certification level, the manufacturer needs to perform multiple flight tests to ensure that all the safety and performance requirements are met: this is referred as the flight test campaign as illustrated in Figure I.1. The certification can vary from a whole aircraft to a single system. During a flight test, a large number of embedded functions in the aircraft is used and most of them require the monitoring and/or the analysis by flight test experts in-flight (online) as well as on ground (offline). These functions involve physical flight parameters that are sampled by specific sensors along time: such kind of measurements gives rise to time series data. These time series depict events like engine start, compression pumping engine or an automatic approach, to name a few. For a specific system and a set of flight parameters given by flight test experts, the time series relative to the flight parameters are often correlated with each other. Indeed, their behavior are induced by a control logic and constrained by a set of physical laws (e.g. from flight mechanics, combustion laws) which, depending on the system of interest, are fully or partially understood.

Figure I. 1 :

 1 Figure I.1: High-level flow-chart description of a flight test campaign (from [Bunge et al. 2015])

Figure I. 2 )

 2 Figure I.2), in this thesis we propose several time series based methods and algorithms which aim to extract knowledge e.g. to assist flight test analysts.

Figure I. 2 :

 2 Figure I.2: Example of an Airbus dataset of multivariate time series relating the behavior, along time (vertical axis), of a system involving four parameters (P 1 , . . . , P 4 ). Grey curves correspond to normal behaviors whereas red ones correspond to abnormal behaviors. For data privacy, the vertical axes have been scaled.

Figure

  Figure I.3: Synthetic dataset of n = 10 univariate functional time series, with time steps in T = t 1 = 0, . . . , t 1000 = 10 , generated by x(t) = t sin(t) + t 2 cos(t)+ for t ∈ T = [0, 10], where is a random variable with distribution N (0, σ = 0.25).

Figure I. 4 :

 4 Figure I.4: Synthetic dataset of n = 10 bivariate time series dataset, generated according to Example I.1, among which three were randomly chosen and colored not to overload the bottom plot. Top-left: the samples of the variable x 1i . Top-right: the samples of the variable x 2i . Bottom: the three (i = 1, 2, 3) colored samples plotted in the space of variables to visualize the correlation between x 1i and x 2i along t.

Figure I. 5 :

 5 Figure I.5: Same synthetic dataset than in Figure I.3 (grey) with outliers (not grey) of three types. The shift outlier, in yellow, was generated by u hoz-shif twith τ drawn from N (6, 0.1). The magnitude outlier, in red, was generated by u magnitude with a drawn from N (-30, 0.1). The shape outlier, in green, was generated according to u(t) + sin(2πf 0 t) with f 0 = 7.

  Figure I.4. For instance, in Figure I.4 the samples i = 1 (blue) and i = 3 (green) are quite similar in their individual variable representation (top) but their distinction is much more obvious in their joint variable representation (bottom).

  Figure II.1, we show a multivariate functional dataset contaminated by one outlier whose variables are non-linearly correlated where: in (a) and (b), the variables x 1 and x 2 are plotted independently, whereas in (c) and (d), they are plotted one versus the other, thus highlighting their correlation along t. Thus, independently

  proposed detecting outliers in univariate functional data by decomposing each univariate functional sample into three features: translation, phase, and amplitude. The authors defined the translation of a functional sample by its mean over the observation interval T . Both the amplitude and phase components are functional data extracted from the original samples. The amplitude component reflects the vertical variability of the functional data, whereas the phase component reflects the horizontal variability. Analogously to the functional boxplot of [Sun & Genton 2011],

  FigureII.3: The norm of the tangent vector D 1 x(t) (red diagonal arrow), the components of which are infinitesimal variations (dx 1 (t), dx 2 (t)), shown by the horizontal and vertical red arrows of the variables of the curve allows the computation of the speed at which the curve "progresses."

  CHAPTER II: SHAPE-BASED OUTLIER DETECTION each variable refers to a flight parameter. This dataset contains 85 outliers annotated by the expert. For data privacy, we cannot share the Airbus flight test data (AFT) neither give technical details on its physical meaning. We show a plot of the scaled time series in Figure I.2 where red curves corresponds to outliers.

  the score is close to 1, then the sample is likely an outlier. OCSVM is a distance-based model formulated as a constrained quadratic minimization problem, the variables of which correspond to the radius and the center of the smallest hypersphere containing the data. To allow flexibility on the hypersphere boundary owing to the presence of outliers in the training data, slack variables are introduced in the objective function in addition to the two other variables. The hyperparameter ν corresponds to an upper bound on the a priori proportion of outliers in the training set. A sample is declared as an outlier if it lies outside the fitted hypersphere. We used the radial-basiskernel version of OCSVM with ν equal to the exact proportion of outliers in the training set. The bandwidth hyperparameter of the radial basis kernel was optimized by a 20-fold cross-validation procedure. For the ECG, PenDig, AFT and synthetic datasets, we set the number of trees to 1000, and the subsampling size to 32 which are recommended values in [Liu et al. 2008]. We randomly split each dataset into a training set and a test set. As in [Dai & Genton 2019], the training set represents 50% of the data for the ECG dataset. The training set for the PenDig and AFT datasets consists in 75% of the entire dataset. The training set contains 60% of the data for the synthetic data. The training set was used to both fit the model (iFor and OCSVM) and select an outlyingness threshold from the Receive Operating Curve (ROC) that discriminates inliers from outliers.

  Figure II.5: The three performance measures ρ c , ρ f , and AU C, averaged over the simulations versus the contamination level (c = 5%, 10%, 15%, 20%, 25%) for each method (proposed in blue, and baselines in red). We notice that when the contamination level c increases, the proposed method (except for iFor(Len map ) and OCSVM(Len map )) outperforms the baselines in terms of ρ c , ρ f and AU C.

  for the synthetic data generated by the five models (columns), as described in Section II.5.1.4. We compared the proposed methods, iFor(•) and OCSVM(•), with the two baselines, F U N T A and Dir.out, in terms of three performance measures (in sub-columns): correct detection rate (ρc), false detection rate (ρ f ), and AU C. For each model and each performance metric, we marked in bold the best results (i.e.highest ρc and AU C, and lowest ρ f ).

  FigureII.7). Each bar has a height equal to its average rank (1 is the best, 8 is the worst) based on the post-hoc Friedman's aligned rank test across the five contamination levels (resp., five models). For ρ c and AU C, the ranking is given in decreasing order, and for ρ f , the ranking is given in increasing order. The above number of bars refers to the global ranking (i.e. ranks from the average ranks).The Friedman test yielded to the rejection of equality of the methods for the ECG, PenDig and synthetic datasets, for each performance measure (p-values are given in the discussion). Therefore, we conducted pairwise comparison (post-hoc) tests. The significance of the pairwise comparison tests of ρ c , ρ f , and AU C for the ECG and PenDig datasets is given in Tables II.5, II.6, and II.7, and for the synthetic data, in Tables II.8, II.9, and II.10. The significance (at level 0.1) of a test is indicated by = * , and non-significance is indicated by = .

F

  corresponding methods in the row and the column of the cell are significantly different at a level of 0.1, and = indicates that they are not. The lower triangular part was replaced by dashes because it is equal to the upper part.

F

F

Figure

  Figure II.7 shows the average ranking of the methods according to the posthoc Friedman's aligned rank test.

F

  Significance of the pairwise comparisons for the false detection rate ρ f on the synthetic dataset.

FF

  U N T A Dir.out iFor(Vmap) iFor(Curvmap) iFor(Lenmap) OCSVM(Vmap) OCSVM(Curvmap) OCSVM(Lenmap) same as in Table II.5 72 CHAPTER II: SHAPE-BASED OUTLIER DETECTION Table II.10: Significance of the pairwise comparisons for AU C on the synthetic dataset.

FFFigure

  FigureII.6: (1 is the best, 8 the worst) for ρ c , ρ f , and AU C based on the post-hoc Friedman's aligned rank test, considering the five contamination levels in the ECG data (upper bar plot) and the PenDig data (lower bar plot). For ρ c and AU C, the ranking is given in decreasing order (i.e. for high ρ c and AU C values, the rank tends to 1); for ρ f , the ranking is given in increasing order (i.e. for low ρ f values, the rank tends to 1). The y-axis represents the average ranking over the five models, and the integers on the top of the bars represent the final ranking. If there are ties, we take the average ranking.

  Figure II.7: (1 is the best, 8 is the worst) for ρ c , ρ f , and AU C based on the post-hoc Friedman's aligned rank test. For ρ c and AU C, the ranking is given in decreasing order (i.e. for high ρ c and AU C values, the rank is close to 1), and for ρ f , the ranking is given in increasing order (i.e. for low ρ f values, the rank is close to 1). The y-axis represents the average ranking over the five models, and the integers on the top of the bars represent the final ranking. If there are ties, we take the average ranking. The baseline methods are slightly better than the proposed method, but the best results by the proposed method (iF or(Len map ) and OCSV M (Len map )) are statistically equivalent to those by the baseline, as showed by the pairwise comparison tests in Tables8,9, and 10.

  variable first order derivative. To illustrate a SODE, we give an example of a two dimensional SODE in Figure III.1 bottom-left wherein one can see that the two equations entail both state variables x 1 and x 2 through cubic monomials and the first equation has in addition a linear monomial. The solution of a governing equation depicts the temporal and/or spatial evolution of the state-variables in a dynamic system. Famous examples of SODEs are for example the FitzHugh-Nagumo equation to model neural excitement in biology, the damped harmonic oscillator equation in mechanics or the Lotka-Volterra equation to model population dynamics [Ramsay & Hooker 2017].

  Figure III.1, we illustrate the core of both our and state of the art methods: (top) based on data ( Ẋn and Θ Xn ) sam-pled from a dynamic system, solving for the optimization problem involving a sparse penalty R (bottom-center) leads to the identification of a SODE (bottom-left). The learned SODE strongly relies on the chosen penalty. We propose a nonconvex sparse penalty to learn the matrix coefficient β that accounts for correlations within the SODE. The nonconvexity of such a regularizer results in a penalty with stronger sparsity than convex regularizers leading to better selection of the candidate functions. Furthermore, this penalty provides in a SODE learned more accurately than with state of the art penalty.

  ods do not scale to large data sets and are inclined to over-fitting. To remedy it, Brunton et al. re-casted the learning of a SODE as sparse regression problem. Such a modern problem formulation has recently formed a new general framework referred as the sparse identification of nonlinear dynamics in the literature.

  CHAPTER III:DATA DRIVEN DISCOVERY OF SODE sion(s) i.e. for any vector or matrix W , R SCAD λ,θ (W ) = i j r SCAD λ,θ (w ij ).

  functions M IS i.e. the number of misidentified functions divided by the total number of candidates functions for every tasks p × m. β measures the unbiasedness w.r.t the true coefficient matrix β * . T measures the relative total error, along t, between the (numerical) solution of the learnt SODE w.r.t the true one. The lower β , T and M IS , the better the recovery of the SODE.

98

  CHAPTER III:DATA DRIVEN DISCOVERY OF SODETableIII.1: Results (average (%) ± standard deviation (%), over ten trials) of the learning of the DOC, LV and LAT SODEs with five regularizers.

100Figure

  Figure III.3: Left: plot in R 2 of the clean solution (i.e. ground truth) (black) of the DOC SODE and the noisy version (blue). The black circle is the initial state x DOC (0). Right: same data plotted in the temporal domain.

FigureFigureFigureFigure

  Figure III.4: Left: plot in R 2 of the clean solution (i.e. ground truth) (black) of the LV SODE and the noisy version (blue). The black circle is the initial state x LV (0). Right: same data plotted in the temporal domain.

  instance, such a combination can be made with ensemble models, wherein each model has to detect outliers from the representation returned by a single mapping function. Another line of research is the outlier class recognition (in addition to the detection), which would enhance the interpretability of our approach. Recognizing the type of outlyingness can inform the user on the abnormality of the behavior e.g. a delay if the outlier is horizontally 107 108 CHAPTER III: GENERAL CONCLUSION shifted w.r.t to the inlier class. This could also be addressed by ensemble models wherein each model is trained on a single outlier class. Then, one can decompose the global outlyingness score in multiple sub-scores, where each one refers to an outlier class. Such a decomposition would inform on the contribution of each outlier class to an outlier.

II.2 Related work II.2.1 Depth-based univariate functional outlier detec- tion

	CHAPTER II: SHAPE-BASED OUTLIER DETECTION
	ii) We propose different mapping functions to capture different types of The detection of outliers in functional data is a recent topic and has pri-outlyingness based on curve shape. marily been addressed by extending statistical depth 1 to functional depth.
	(iii) We demonstrate that the proposed method is superior to previous Statistical depth measures the centrality of a sample relative to a dataset by
	outlier-detection algorithms and, in contrast to baseline methods, per-providing an outward-center ordering of the samples through a score lying in
	forms well regardless of the proportion of outliers. [0, 1]. A value close to zero implies that the sample is more likely to be an
	outlier [Tukey 1975]. Statistical depth has several theoretical properties (see This chapter is organized as follows. In Section II.2, we review related work on outlier detection in both univariate and multivariate functional data. In [Zuo & Serfling 2000] for details):
	Section II.3, we discuss curve representations in the functional-data frame-(i) It attains its maximum value for the most centered (i.e. most repre-
	work. In Section II.4, we present the mapping functions that can capture sentative) sample.
	shape outlyingness from the obtained functional representation. The exper-
	imental results are presented and discussed in Section II.5. Finally, Sec-
	tion II.6 concludes this part.

  (t), y 2 (t)) with y 1 (t) = u 1 (t) + z 11 cos(4πt) and y 2 (t) = u 2 (t) + z 12 sin(4πt), where z 11 and z 12 are independent uniform random variables in[2, 3]. The contamination model x(t) c5 is x 1,c5 (t) = u 1 (t) + z 21 cos(4πt) and x 2,c5 (t) = u 2 (t) + z 22 sin(4πt), where z 21 , z 22 are uniform random variables on[4, 5].

	CHAPTER II: SHAPE-BASED OUTLIER DETECTION
	5. Model 5 persistent shape outlier, the new uncontaminated model is
	referred to as Y , and the contamination model as x(t) c5 : y(t) =
	(y 1
	and x 2,c3 (t) = 1.5u 2 (t).
	4. Model 4 isolated outlier: x c4 (t) = u(t)(1 + 4I z<t<z+0.1 ), with z as in
	model 2.

Table II .

 II 2: Results for the PenDig dataset

		AU C	0.53	0.58	0.77	0.63	0.56	0.74	0.66	0.69		AU C	0.53	0.61	0.64	0.73	0.63	0.70	0.72	0.67
	Outlier '4'	ρ f	0.21	0.02	0.09	0.33	0.45	0.18	0.15	0.25	Outlier '9'	ρ f	0.21	0.34	0.51	0.26	0.40	0.25	0.14	0.46
		ρc	0.23	0.00	0.74	0.81	0.78	0.78	0.61	0.83		ρc	0.23	0.43	0.92	0.79	0.74	0.70	0.62	0.88
		AU C	0.52	0.55	0.66	0.69	0.64	0.66	0.68	0.67		AU C	0.51	0.60	0.67	0.77	0.77	0.63	0.66	0.60
	Outlier '3'	ρ f	0.20	0.02	0.45	0.38	0.23	0.41	0.44	0.42	Outlier '8'	ρ f	0.01	0.17	0.30	0.20	0.08	0.45	0.44	0.24
		ρc	0.23	0.00	0.61	0.57	0.29	0.53	0.56	0.54		ρc	0.23	0.43	0.47	0.48	0.40	0.58	0.56	0.51
		AU C	0.58	0.60	0.63	0.57	0.65	0.60	0.59	0.76		AU C	0.58	0.65	0.60	0.68	0.65	0.60	0.62	0.61
	Outlier '2'	ρ f	0.19	0.42	0.15	0.47	0.12	0.35	0.23	0.18	Outlier '7'	ρ f	0.00	0.16	0.30	0.21	0.29	0.22	0.19	0.13
		ρc	0.22	0.75	0.86	0.87	0.59	0.77	0.55	0.84		ρc	0.22	0.43	0.93	0.93	0.97	0.87	0.86	0.79
		AU C	0.51	0.52	0.79	0.61	0.64	0.75	0.70	0.70		AU C	0.51	0.52	0.61	0.63	0.64	0.56	0.61	0.57
	Outlier '1'	ρ f	0.21	0.02	0.38	0.60	0.56	0.51	0.60	0.42	Outlier '6'	ρ f	0.02	0.17	0.36	0.28	0.21	0.38	0.40	0.47
		ρc	0.01	0.24	0.44	0.43	0.46	0.50	0.50	0.37		ρc	0.01	0.43	0.56	0.54	0.47	0.55	0.61	0.62
		AU C	0.60	0.82	0.87	0.92	0.59	0.85	0.91	0.75		AU C	0.60	0.59	0.69	0.61	0.61	0.73	0.64	0.62
	Outlier '0'	ρ f	0.22	0.01	0.05	0.12	0.26	0.02	0.11	0.10	Outlier '5'	ρ f	0.22	0.34	0.26	0.29	0.13	0.04	0.18	0.30
		ρc	0.49	0.72	0.78	0.82	0.63	0.82	0.80	0.81		ρc	0.49	0.43	0.69	0.62	0.42	0.59	0.58	0.67
		Methods	F U N T A (baseline)	Dir.out (baseline)	iFor(Vmap)	iFor(Curvmap)	iFor(Lenmap)	OCSVM(Vmap)	OCSVM(Curvmap)	OCSVM(Lenmap)		Methods	F U N T A (baseline)	Dir.out (baseline)	iFor(Vmap)	iFor(Curvmap)	iFor(Lenmap)	OCSVM(Vmap)	OCSVM(Curvmap)	OCSVM(Lenmap)

  Table II.3: Results on the AFT dataset.

	Methods	ρ c	ρ f	AU C
	F U N T A (baseline) 0.85 0.10 0.84
	Dir.out (baseline)	0.90 0.12 0.90
	iFor(V map )	0.98 0.01 0.97
	iFor(Curv map )	0.5	0.31 0.60
	iFor(Len map )	0.97 0.02 0.96
	OCSVM(V map )	0.95 0.03 0.96
	OCSVM(Curv map ) 0.52 0.34 0.55
	OCSVM(Len map )	0.95 0.05 0.92

  Table II.5: Statistical significance of the pairwise comparisons for the correct detection rate ρ

c on the ECG (upper table) and PenDig (lower table) datasets.

  Table II.6: Statistical significance of the pairwise comparisons for the false detection rate ρ f on the ECG (upper table) and PenDig datasets (lower table).

  Table II.5.

	II.5. EXPERIMENTAL STUDY	67
	on the
	Table II.7: Statistical significance of the pairwise comparisons for the AU C	ECG (upper table) and PenDig (lower table) datasets.

  CHAPTER III:DATA DRIVEN DISCOVERY OF SODEh : H → R, is defined for every y ∈ H, with λ > 0 as: Here, depending on the context, H reduces either to R p or R n×p . By the strong convexity of the two terms in the optimization problem involved in Definition III.1, the optimizer prox λh (u) is unique and thus is single-valued in H[Parikh & Boyd 2013].

	prox λh (y) = arg min u∈H	1 2	y -u 2 H + λh(u)
	Remark III.3.0.1.
	W )
	7: end while
	Output: β
	Definition III.1. (proximal operator [Parikh & Boyd 2013]) The proximal
	operator associated to a closed proper convex function of a Hilbert space,

Remark III.3.0.2. When h is separable, i.e. for any vector or matrix W , h(W ) = i j h ij (w ij ) with h ij : R → R, computing prox λh (W ) reduces to compute the proximal operator of h ij for every i, j and then to concatenate prox λh ij (w ij ) ij according to the dimensions of W . In other words, in the case of a separable vector or matrix function, evaluating the proximal operator on a given element boils down to evaluate the proximal operator for each of the separable parts. Hence, since the 1,1 matrix norm is separable, h ij (w ij ) = |w ij |, its proximal operator acts entrywisely and prox λ|•| (w ij ) = sign(w ij ) max (0, |w ij | -λ) which is known as the soft-

  and can degrade the identification of the SODE. Indeed, the soft-thresholding operator evaluated ona regression coefficient, i.e. prox λ|•| (β ik ) = βik = sign(β ik )(|β ik | -λ) if |β ik | -λ > 0.Hence, if the i-th candidate function is relevant for the k-th equation of the SODE and if the sparsity amount λ must be large, the learned coefficient βik is underestimated (in absolute value) with a bias amount of λ sign β ik .

R

  SCAD-1 enforces the coefficients of the i-th candidate function to be sparse, unbiased and correlated across the p tasks. Note the analytical similarity with the R 2,1 = i β i• 2 norm (Section III.3.2) which does not allow each equation of the SODE to have a specific candidate function. Contrary to R 2,1 , our proposal R SCAD-1 enforces each β i• to have a small 1 norm and thereby enables the components of the SODE to have specific candidate functions (this is due to the sparsity induced by the 1 norm). RSCAD-1 Algorithm III.2: GIST for Problem (III.4) with R possibly nonconvex Input: data samples Ẋn and X n , sparsity amount λ > 0, unbiasedness level θ > 2, initial guess β 0 , step size γ, step size minimum γ min . 1: build Θ Xn arbitrarily from X n 2: β ← β 0 , γ ← γ min 3: W ← β -γ min Θ Xn (Θ Xn β -Ẋn ) 4: while β has not converged do

	5:	while line search criterion unsatisfied do
	6:	β ← prox γλR (β)
	7:	γ ← 0.8γ
	8:	end while
	9:
	10:	β ← prox γλR (W )
	11:

W ← β -γΘ Xn (Θ Xn β -Ẋn )

  Table III.2: Results (%) of the SODE discovery experiment on an Airbus flight test dataset with five regularizers. The two smallest errors are in bold.

	SODE	Flight test dataset (p = 2)
	Reg/Error
		67.41
	R 2,1 (MTLa)	55.87
	R 2,1 + 1,1	55.29
	R SCAD	62.56
	R SCAD-1 (our)	44.96

T R 1,1 (LASSO)

statistical depth was not specifically proposed for functional but for multivariate data[Tukey 1975, Fraiman & Muniz

2001]. However, we distinguish between univariate functional depth and multivariate functional depth, which were proposed specifically for functional data.

Remerciements

CHAPTER II: SHAPE-BASED OUTLIER DETECTION

for the false detection rate ρ f , and 2.2 × 10 -16 for AU C. Thus, we conducted a post-hoc test. Figure II.6 shows the average ranking of the methods based on the Friedman's aligned rank test (from the best 1 to the worst 8). The pvalue of each pairwise comparison in the post-hoc test is given in Tables II.5, II.6,and II.7 for the correction detection rate, the false detection rate, and AU C, respectively, where a cell indicates whether the resulting p-value of the pairwise comparison test of the methods in the corresponding row and column is below of above the significance level. The symbol = indicates a p-value greater than the significance level of 0.1, allowing the acceptance of the null hypothesis of equivalence of the two methods; rejection is indicated by = * .

Based on the results in Figure II.6 and Tables II.5,II.6,and II.7, it is concluded that both V map and Curv map outperform the baselines in terms of the three performance measures. We notice that Dir.out is not significantly better than the methods with the worst performance (i.e. iFor(Len map ), F U N T A, and OCSVM(Len map )). F U N T A is not significantly different from iFor(Curv map ) and OCSVM(Curv map ) (Tables II for the false detection rate, and 1.1 × 10 -4 for AU C. We note that there

III.2 Related work

We focus on the case where the state-variables are sampled from a dynamic system along a scalar variable t ∈ R that, without loss of generality, refers to time. Thus each one of the sampled state-variables form a time series.

To discover the dynamic relationship between the state-variables, one can assume that the data in underlied by the solution of an unknown SODE as proposed by [START_REF] Brunton | [END_REF]. Discovering the SODE from the data enables the practitioner to better understand the underlying dynamic of the phenomenon. We emphasize that, contrary to recent work on deep learning algorithm driven by SODE solver [Chen et al. 2018], we aim to discover a SODE from data samples and not to solve a given SODE since the later is unknown by definition of the problem.

For clarity, let's recall how a SODE is formalized. Let f = f 1 , . . . , f p : R p → R p be a Lipschitz continuous map defining the evolution of a statevariable

, where ẋ refers to the first-order time derivative of x (we use the 'dot' notation instead of D for a better readiness), or equivalently:

III.4 Nonconvex matrix-structured regularizer

In this section, we expose our contribution and propose a nonconvex regularizer which considers both relatedness between the tasks, through sparsity, and unbiasedness. To learn a SODE with such a regularizer, we instantiate a generative iterative thresholding algorithm [Gong et al. 2013] which can be used with convex as well with nonconvex regularizers which can be expressed as the difference of two convex functions.

III.4.1 Nonconvex separable regularizer

We first introduce the single task learning version of the nonconvex regularizer that we extend to MTL by applying it entrywisely on a matrix, i.e. as the 1 norm to 1,1 norm. This nonconvex regularizer, the Smoothly Clipped Absolute Deviation (SCAD) [Fan & Li 2001], does not account for task relatedness but serves as a building block for the regularizer we propose and incorporates unbiasedness.

Definition III.2. Let λ > 0 and θ > 2 be two hyperparameters that serve as the sparsity amount and unbiasedness level respectively. The SCAD penalty is defined for any w ∈ R as:

is not differentiable for x = 0 and nonconvex (actually concave in R + ). Like any convex regularizer used in a learning algorithm, the non-differentiability at 0 involves sparsity which is a desirable property. Also, when θ increases, the SCAD regularizer approximates the 1 norm [Fan & Li 2001]. given in [Fan & Li 2001] by: can be easily built by applying it entrywisely and summing the resulting outputs along the input dimen-was reformulated as a sparse group penalty to learn single-task linear models when groups of correlated predictive variables are known but some of the groups may have been misspecified [Jiang & Huang 2015]. To the best of our knowledge, this regularizer has never been extended to MTL.

Proximal operators

Like R SCAD λ,θ

, the regularizer R SCAD-1 λ,θ is nonconvex but its proximal operator can be computed according the 1 norm of the rows of

as, for i ∈ {1, . . . , m}:

where sign, max and | • | act entrywisely. We will see in the next section that knowing, analytically, this proximal operator is essential to perform the learning.

III.4.3 MTL with a nonconvex regularizer

In this section, we instantiate the Generative Iterative Sequential Thresholding (GIST) algorithm of [Gong et al. 2013] to learn β ∈ R m×p with nonconvex regularizers that can be expressed as a difference of two convex functions, see [Gasso et al. 2009] for examples of such regularizers. We present the GIST algorithm in Algorithm III.2 with a regularizer R that can be either convex, either nonconvex and can be written a difference of two convex functions, see [Gong et al. 2013] for convergence guarantee of GIST. This To show the efficiency of our approach, we evaluated it on three reference two-dimensional SODEs: the Damped Oscillator with Cubic dynamic (DOC) (used to model damped behaviors with non-linearity), the Lotka-Volterra (LV) system (used for predator-prey interactions modeling) and the Lorenz Attractor (LAT) system (used for excitation systems modeling like neurons). Each one of these systems has common functions as well as specific functions across their two equations. We used the same settings as [START_REF] Schaeffer | [END_REF]. We generated the time series by solving numerically the true SODEs with the explicit Runge-Kutta-45 method. We corrupted the time series and their time derivative with a Gaussian noise (σ = 0.1). Based on these noisy samples, we built the dictionary with monomial candidate functions (up to degree five, with first and second order interactions e.g. x 1 x 2 , x 2 1 x 2 ).

For clarity, we give the analytic form of the DOC, LV and LAT SODE:

LV

LAT

III.5.1.2 Airbus flight test data

We experimented the baselines and our method on a multivariate time series provided by an Airbus flight test expert. For data confidentiality, we cannot share details on the data meaning. The dataset includes the time series of the sampled state-variables (Figure III.7, blue curves), as well as their sampled time derivatives (Figure III.8 blue curves). We built the dictionary with monomial candidate functions (up to degree five, with first, second and third order interactions terms e.g. x 1 x 2 , x 2 1 x 2 , x 3 1 x 2 ). Based on the knowledge of a flight test expert, we also included cos(2πx 1 ), cos(2π0.5x 1 ), sin(2πx 1 ) and sin(2π0.5x 1 ). Also, to account for temporal variability of the SODE, regardless of the state-variables, we included exponential terms e -0.5t , e -t , e -1.5t , e -2t in the dictionary. Since it is a real dataset, the true SODE is unknown thus we can only report the T error between the sampled state-variables and the numerical solution of the learnt SODE. For data confidentiality, we cannot report the analytical SODE found but a plot of the numerical solution. To solve the SODE, we used the explicit Runge-Kutta-45 method implemented in the SciPy library [Oliphant et al. 2001].

III.5.1.3 Implementation

We learnt the SODEs with GISTA with R SCAD-1 and R SCAD . We set θ = 2.01, near its lower bound. The best value of λ was computed as the one III.5. NUMERICAL EXPERIMENTS