Piet Mondrian
email: @

Multi-Block representation

Keywords: . Nullam nec mi et neque pharetra sollicitudin. Praesent imperdiet mi nec ante. Donec ullamcorper, felis non sodales commodo, lectus Bayesian Nonparametric Multi-Coclustering Bayesian Nonparametric Coclustering Parametric Classification Non Supervisée, Classification par blocs, Classification Multiple, Séries temporelles, Modèles de mélange, Modèles Bayésiens Non Paramétriques, Validation de Systèmes d Clustering, Coclustering, Multi-Clustering, Time Series, Mixture Models, Bayesian Non-Parametric Models, Advanced Driver-Assistance Systems validation

Advanced driver-assistance systems (ADAS) development remains one of the biggest challenges car manufacturers must tackle to provide safe driverless cars. The reliable validation of these systems requires assessing their reaction's quality and consistency to a broad spectrum of driving scenarios. In this context, Groupe Renault uses large-scale simulation systems, which accurately reproduce the physical driving conditions and produce large quantities of high-dimensional time series data. The role of the ADAS developer is to explore these datasets to determine precisely the capabilities of the driving assistance system studied and, if necessary, to refine its design.

A simulated data set can contain up to several hundred thousand simulations for a given use case, described by several hundred variables. The expert's work requires a great deal of business knowledge and a meticulous time-consuming investigation with datasets of such dimensions to discriminate. This thesis aims to produce algorithms and tools to help explore, structure and analyze these massive simulation data sets. We propose four probabilistic approaches for time series clustering, each based on a specific datasets structure hypothesis.

The first contribution of this thesis is built on a scenario-based construction hypothesis and performs a dictionary-based classification of univariate time series. This method uses an existing piecewise polynomial segmentation model to construct a driving pattern dictionary and recode the time series into chains of driving phases. The subsequent contributions focus on the analysis of multivariate time series datasets with multi-blocks methods that regroups variables based on their distribution and partitions.

Applications from Renault industrial use cases illustrate the usefulness of these contributions : validating an emergency braking system, a lane-keeping system, and an emergency avoidance system.

Résumé

Le développement de systèmes d'aide à la conduite demeure un défi technique pour les constructeurs automobiles. La validation de ces systèmes nécessite de les éprouver dans un nombre considérable de contextes de conduites. Pour ce faire, le Groupe Renault a recourt à la simulation massive, qui permet de reproduire précisément la complexité des conditions physiques de conduite et produit une grande quantité de séries temporelles multivariées. Le rôle de l'expert métier est alors d'explorer ces données pour déterminer précisément les capacités du système d'aide à la conduite étudié et, si nécessaire, de raffiner sa conception.

Un jeu de données simulées peut contenir jusqu'à plusieurs centaines de milliers de simulations, décrites par plusieurs centaines de variables. Face à de telles dimensions, le travail du développeur requiert une grande expertise, ainsi qu'un travail d'investigation minutieux et chronophage. L'objectif de cette thèse est de produire des algorithmes et outils d'aide à l'exploration, la structuration et l'analyse de ces jeux de données issues de simulation massive. Basées sur des hypothèses de construction argumentées, nous proposons trois approches probabilistes de classification non supervisée de séries temporelles, adaptées à des jeux de données temporelles univariées et multivariées.

La première contribution de cette thèse est construite sur une hypothèse de construction par scénario, et réalise une classification de séries temporelles univariées par dictionnaire. Partant d'un modèle existant de segmentation polynomiale par morceaux, cette méthode crée un dictionnaire de régimes de conduite pour recoder les séries temporelles en phases similaires. Les contributions suivantes se concentrent sur la classification de jeux de données multivariés par bloc, qui regroupent les variables en se basant sur leur distribution et leur partitionnement.

La suite des contributions se concentre sur l'analyse de jeux de données multivariés, grâce à des modélisations multi-blocs qui regroupent les différentes variables en fonction de leur distribution.

La pertinence de ces contributions est illustrée par des applications issues des cas d'usage industriels Renault : la validation d'un système de freinage d'urgence, d'un système d'aide au maintien dans la voie, et d'un système d'évitement d'urgence.

viii

Table des matières

Introduction

Introduction

" De plus en plus, les technologies connectées et la conduite autonome permettront aux usagers de faire un choix entre conduire et être conduits. Ce constat ouvre la voie à de nouveaux scénarios de mobilité. Voici notre prochaine révolution.

-Guillaume Eurin Directeur du Développement des ADAS et du Véhicule Autonome . Groupe Renault Rapport Annuel (2019-2020) Cette thèse a été organisée en convention entre le constructeur automobile Groupe Renault 1 et le laboratoire d'informatique LIPN 2 (Laboratoire dInformatique de Paris Nord, UMR CNRS 7030) de l'Université Sorbonne Paris Nord 3 . Ce chapitre présente le contexte industriel qui a motivé cette convention de recherche : la validation de systèmes d'aide à la conduite (ADAS) et ses contraintes opérationnelles. Cette description permet de détailler le rôle des ADAS dans le fonctionnement d'un véhicule, l'interêt de leur développement pour le constructeur. Elle présente également la nécessité de l'utilisation de la simulation massive pour la validation et la nature du travail d'analyse de données associée aux jeux de données simulées. Enfin, cette partie présente également les contributions, en réponses aux contraintes opérationnelles et les publications associées, avant de conclure par le plan global du manuscript.

Systèmes d'aide à la conduite (ADAS)

Les systèmes d'aide à la conduite (ADAS) sont des systèmes électroniques embarqués qui fournissent une assistance pour certaines tâches associées à la conduite. Depuis les premiers systèmes d'Anti-Blocage des Roues (ABS) et de Vitesse de Croisière développés dans les années 1950, jusqu'à la très récente fonction de Pilote Autoroutier Automatique, le développement de ces systèmes a accompagné le marché automobile durant toute son histoire. Toute voiture récente en est aujourd'hui équipée.

Tandis que certains systèmes ne sont que produire une information ou un avertissement (par exemple les systèmes de navigation ou le rappel de non-bouclage des ceintures de sécurité, . . . ), certains fournissent une aide active à la conduite quotidienne (par exemple, la vitesse de croisière adaptative, l'allumage des feux automatique, . . . ). D'autres, enfin, fournissent une aide vitale lors de situations d'urgence (par exemple le freinage d'urgence, le maintien dans la voie, . . . ). Les ADAS peuvent être classés en fonction de leur niveau d'automatisation de conduite, sur une échelle de 0 à 5 développée par la Society of Automotive Engineers (SAE). Cette échelle, representée sur la Fig. .1, peut également être vue comme une feuille de route menant au développement de la voiture autonome, et dont les ADAS constitueraient les fondations. La liste suivante présente trois exemples d'ADAS existants ou en développement :

-Le Freinage d'Urgence Automatique (AEB) -Niveau 0. Durant la conduite, ce freinage se déclenche lorsque la voiture équippée (appelée Ego) détecte une cible (piéton, voiture, camion, vélo, . . . ) se rapprochant à grande vitesse.

-L'Aide au Maintien dans la Voie (LKA) -Niveau 1. Durant la conduite, ce système détecte une dérive non intentionnelle d'Ego en direction d'une autre voie ou du bord de la route et corrige sa trajectoire.

-Valet Parking -Niveau 4 (En développement). Dans une aire de stationnement, ce système dirige le véhicle sans intervention extérieure, cherche et trouve un espace de stationnement disponible et réalise la manoeuvre de stationnement.

Pour le conducteur, ces ADAS sont des gages de sécurité, de confort et d'efficacité. D'un point de vue plus global, les ADAS sont pensés comme une forme d'aide à la mobilité, d'assistance aux personnes agées, et une solution pour réduire le traffic automobile [START_REF] Miranda | Autonomous driving-political, legal, social, and sustainability dimensions[END_REF]. D'un point de vue commercial, ils constituent des arguments de vente décisifs qui permettent aux constructeurs de prouver leur savoir-faire et de se démarquer sur un marché de l'automobile extrêmement compétitif. Dans un récent rapport de recherche [START_REF]Advanced Driver Assistance System (ADAS) Market Size, Share, Growth, Trends, Regional Outlook and Forecasts[END_REF], le cabinet d'étude marketing Precedence Research estimait la valeur du marché mondial de l'ADAS à 32.

Table des matières

d'un ADAS. Pour s'en assurer, il est donc nécessaire de mettre en place un système de validation rigoureux qui vérifiera les performances de l'ADAS dans tous les contextes de conduite.

Validation

L'objectif de l'étape de validation est de confirmer que les réactions d'un ADAS sont conformes au cahier des charges fonctionnel, rédigé par le Groupe Renault, en accord avec les recommandations des agences gouvernementales de sécurité (par exemple Euro NCAP en Europe)

Dans ce contexte, pour la validation d'un ADAS hautement autonome (Niveau 4), le Groupe Renault requiert un niveau minimum de fiabilité de 10 -8 décès par heure, ce qui correspond à 10 fois moins que la moyenne du nombre d'accidents sur autoroute en Europe, et 100 fois moins que la moyenne d'accident sur tous types de route. Vérifier un tel niveau de fiabilité demanderait de réaliser au moins 2 × 10 10 kilomètres de conduite [START_REF] Didier | Autonomous Vehicles Digital Simulation Challenges[END_REF], soit plus de 100 années sur la route.

Ce niveau de fiabilité doit être atteint pour chaque véhicule de la gamme, pour chaque ADAS, et dans tous les cadres législatifs mondiaux existants. De plus, les développeurs doivent également s'assurer que les ADAS sont inter-opérables, et que leur combinaison ne provoque pas d'effets de bord (par exemple : comment combiner le freinage d'urgence et l'évitement d'urgence). L'accumulation de tous ces facteurs de complexité rend impossible la validation basée uniquement sur des tests physiques sur piste ou sur route.

Pour résoudre ce problème et parvenir à tenir ces engagements de fiabilité, Le Groupe Renault a fait le choix d'investir dans la technologie de simulation massive de conduite. Ces systèmes de simulation reproduisent les conditions de conduite en se basant sur un scénario, la physique du véhicule, le comportement du conducteur, et les interactions avec un environnement lui-même configurable. Le procédé de simulation consiste en un processus distribué qui génére des scènes physiques en parallèle. La qualité de la simulation, c'est-à-dire sa capacité à reproduire fidèlement les conditions de conduite réelles, est méthodiquement monitorée afin de s'assurer de la pertinence du procédé de simulation. Le système de simulation est representé sur la Fig. La simulation massive produit de grandes quantités d'information provenant de l'observation de variables simulés, et aggrégée en un unique jeu de données de séries temporelles multivariées. Le rôle de l'expert métier est alors d'analyser ce jeu de données pour discriminer des comportements de conduite, détecter des situations aberrantes, et de comparer les performances de l'ADAS au cahier des charges.

Les dimensions et la complexité des données sont considérables ; pour un cas d'usage donné, le nombre de simulations peut atteindre O (10 6 ), être décrit par O(10 3 ) variables, chacune en O (10 4 ) points temporels. Au total, plus de O (10 13 ) observations peuvent être produits pour tester le paramétrage d'un ADAS (ce qui représente jusqu'à 1 téraoctet). Les simulations ont des durées différentes liées aux conditions de scénario de conduite testée, par conséquent les tailles des séries temporelles simulées sont également différentes. Le jeu de données final peut être représenté sous la forme d'une matrice, contenant une série temporelle dans chaque cellule, de taille spécifique à chaque simulation. Un exemple fictif de matrice simulée est représentée sur la Fig. De plus, si certains résultats de simulation peuvent être prédits (en particulier, le comportement attendu), la plupart sont découverts après la simulation. Par conséquent, il n'est pas possible d'établir à l'avance la liste exhaustive de labels qui permettrait de réaliser une analyse supervisée. Dans ce contexte en grandes dimensions, multivarié et non-supervisé, extraire de l'information pertinente serait à la fois chronophage et laborieux pour l'expert, sans l'aide d'outils spécifiques d'analyse de données. Les méthodes qui font une classification automatique de séries temporelles dans des groupes homogènes, sans supervision, sont appelées des méthodes de classification nonsupervisée de séries temporelles, ou Clustering de séries temporelles. La Fig.

. 4 représente l'exemple d'une telle méthode appliquée à jeu de données de séries temporelles multivariées. Sur cet exemple, le graphique de droite montre les séparations entre différents groups de simulations, après une projection dans un espace 2D.

Contributions

La classification non supervisée de séries temporelles est une tâche complexe, qui dépend fortement des hypothèses formulées par l'analyste sur le jeu de données étudié. En l'absence de labels et de supervision, il n'est jamais possible de savoir si le résultat du clustering est correct, puisque, par construction, il n'y a aucune façon de savoir ce que « correct » signifie. Le problème est à considérer dans l'autre sens : lorsque l'analyste applique une méthode de clustering sur un jeu de données, c'est la compréhension de la méthode utilisée et des hypothèses sous-jacentes qui permet d'interpréter les résultats et d'en tirer des conclusions.

La première contribution de cette thèse se base sur une hypothèse de construction de séries par scénario. Observant que les simulations sont générées sur la base de scénario de conduite (c.f., Fig. .2), nous avons développé une méthode basée sur la détection de différentes phases de conduite. Durant l'analyse, les séries temporelles sont recodées sous formes de chaînes de régimes de conduite, qui sont des schémas de conduite issues de l'observation du jeu de données globales et aggrégées en dictionnaire. Cette méthode de clustering probabiliste permet également la détection d'anomalies basée sur des intervalles de confiance probabilistes, et l'inférence de valeurs manquantes.

Cependant, par construction, cette méthode ne s'applique qu'à des séries temporelles univariées, ce qui correspond, dans notre cas, à l'analyse d'un seul variable. Cependant, dans la plupart des cas, le comportement d'un ADAS influe sur plusieurs dimensions, observées par plusieurs variables. Il est donc raisonnable de vouloir analyser plusieurs variables simultanément. Les contributions suivantes répondent à cette problématique et forment un groupe cohérent de méthodes dédiées au clustering de séries temporelles multivariées. Elles partagent la même vision multi-bloc et la même racine commune : le Coclustering. Le Coclustering infère simultanément une partition de simulations (appelées partition ligne) et une partition de variables (appelées partition colonne), comme illustré sur la Fig.

. 5 -Une méthode de Coclustering Bayésien Non-Paramétrique, qui intègre nativement une étape de sélection de modèle.

-Une méthode de Multi-Clustering paramétrique qui permet d'inférer une partition colonne et une partition ligne pour chaque cluster colonne. (c.f., Fig. .5-b).

-Une méthode de Multi-Coclustering non-paramétrique. Cette méthode combine le Multi-Clustering en y intégrant des structures de Coclustering (c.f., Fig. .5-c).

Cette approche crée plusieurs vues d'un même jeu de données tout en réduisant le nombre de paramètres, le tout dans un cadre bayésien non-paramétrique qui permet de résoudre la problématique de sélection de modèle.

Le code source associé à ces distributions, majoritairement écrit en Scala, avec des visualisations en R et Python, est également mis à disposition, ainsi que les jeux de données publics de test. Le lien vers les dépôts Github sont rassemblés dans le tableau .1. Par souci de confidentialité, les données des cas d'usage Groupe Renault ne sont pas mises à disposition.

Plan de la thèse

La suite du manuscript se compose des articles suivants :

Chapter 1 : Related work

Le premier chapitre présente les travaux existants sur le sujet de la classification de séries temporelles, univariées ou multvariées. Ce chapitre introduit en particulier les notions de clustering et coclustering probabilistes, ainsi que la modélisation Bayésienne Non Paramétrique, notions qui seront ré-utilisées régulièrement dans les chapitres suivants.

Chapter 2 : Time Series Clustering with Model-Based Dictionary Representation.

Ce chapitre détaille la première contribution de cette thèse, une méthode de clustering de séries temporelles univariées construite en trois étapes : a) segmentation individuelle

Table des matières

par régression polynomiale par morceaux ; b) construction d'un dictionnaire de phases de conduites par modèle de mélange ; c) clustering de chaînes de caractères.

Chapter 3 : Multivariate Time Series And Functional Conditional Latent Block Model

Ce chapitre présente la première méthode de clustering de séries temporelles multivariées, avec une méthode de Multi-Clustering par blocs. Cette approche paramétrique probabiliste de Multi-Clustering étend le Coclustering de séries temporelles en permettant la création de plusieurs vues d'un même jeu de données multivariées.

Chapter 4 : Bayesian Nonparametric Coclustering

Cette partie présente une seconde extension du modèle de coclustering de séries temporelles multivariées. Cette extension ajoute un prior Bayésien Non Paramétrique aux composants du Coclustering, sous la forme d'un modèle de mélange de Dirichlet symétrique. Cette extension permet d'intégrer nativement une sélection de modèle durant l'inférence.

Chapter 5 : Multivariate Time Series Multi-Coclustering

Basée sur les méthodes détailles dans les chapitres 3 et 4, le Multi-Coclustering combine Multi-Clustering et Coclustering. Cette approche consiste en trois modèles de mélange de Dirichlet hiérarchiquement imbriqués et permet l'inférence d'une structure de Multi-Clustering global, puis d'un Coclustering à l'intérieur de chaque vue.

Chapter 6 : Conclusion And Perspectives

Ce dernier chapitre présente de possibles extensions des méthodes proposées. Ces propositions concernent l'extension multivariée de la méthode de clustering par dictionnaire, et la recherche de méthode de représentation alternative pour la méthode de clustering multivarié par bloc.

Publications

Ces contributions ont été l'objet de plusieurs publications au cours de la thèse, représentées sur la Fig.

.6 et listées ci-dessous :

-Goffinet, É., Lebbah, M., Azzag, H., Giraldi, L., & Renault, S. A. S. (2020). Clustering de séries temporelles par construction de dictionnaire. In EGC (pp. 181-192). -Goffinet, E., Lebbah, M., Azzag, H., & Giraldi, L. (2020). Autonomous Driving Validation with Model-Based Dictionary Clustering. In ECML/PKDD.

- -Goffinet, É., Lebbah, M., Azzag, H., Giraldi, L., & Coutant. A (2020). Multiple Coclustering de séries temporelles. Application à la validation de systèmes d'aide à la conduite. In 52emes Journées de Statistiques. Société Française de Statistique.

-Goffinet, É., Lebbah, M., Azzag, H., Giraldi, L., & Coutant. A (2021). In International Workshop on Advanced Analysis and Learning on Temporal Data @ ECML-PKDD .

-Goffinet, E., Lebbah

Introduction

" Increasingly, connected technologies and autonomous driving will allow users to choose between driving and being driven. This fact opens the way for new scenarios for mobility. Thats our next revolution.

-Guillaume Eurin Director of Autonomous Vehicle and ADAS Development. Groupe Renault Annual Report (2019)(2020) This thesis is a bi-party research project in collaboration between the car manufacturer Groupe Renault 4 and the computer science lab LIPN 5 (Laboratoire dInformatique de Paris Nord, UMR CNRS 7030) at Université Sorbonne Paris Nord 6 . This chapter describes the industrial context that motivated this research agreement : the advanced driver assistance systems (ADAS) validation and its operational constraints. This description details the ADAS role in a vehicle and why their development is essential for the car manufacturer. It also explains the necessity of using massive simulation for validation and the nature of the data analysis work associated with the simulated data sets. The last part of this introduction presents the thesis contributions that were designed to address the datasets operational constraints, and the associated publications, before concluding with the manuscript plan.

Advanced Driver-Assistance Systems

Advanced Driver assistance systems are on-board electronic systems that assist with specific driving tasks. From the first Anti-lock Braking System (ABS) and Cruise Control 4. https ://www.renaultgroup.com/ 5. https ://lipn.univ-paris13.fr/ 6. https ://www.univ-paris13.fr/ systems developed in the 1950s, to the very recent Automatic Highway Pilot functionality, the ADAS development has been associated to the automotive market growth throughout its history. Every modern contains some of these systems. Some systems are simple information or warning providers (e.g., seat belt reminder, navigation, . . . ), others actively assist the user in everyday driving tasks (e.g., adaptive cruise control, automated lighting,. . . ). Some others bring vital help in emergencies situations (e.g., emergency braking or lane-keeping, . . . ).

The ADAS can be classified in driving automation levels on a scale designed by the Society of Automotive Engineers (SAE), from level 0 to 5. This scale, displayed in Fig. .7, also represents the roadmap toward full Autonomous Driving, with ADAS serving as keystones. The following list describes three existing or in-development ADAS applications :

-Autonomous Emergency Breaking (AEB) -Level 0. While driving, an emergency brake is triggered when a target (pedestrian, car, truck, bike, . . . ) is detected moving at high speed towards the equipped car (called Ego) in the opposite direction.

-Lane Keeping Assist (LKA) -Level 1. While driving, the system detects an unintentional slow drift of Ego towards another lane or the roadside and corrects the trajectory.

-Valet Parking -Level 4 (Under development). Inside a parking lot, this automated parking system drives the vehicle, finds a parking spot, and performs the parking maneuver.

From the driver perspective, the ADAS enhance vehicle safety, comfort and efficiency. On a broader level, they are also meant to extend mobility, provide assistance for the elderly and reduce traffic congestion [START_REF] Miranda | Autonomous driving-political, legal, social, and sustainability dimensions[END_REF]. From a marketing perspective, they are selling points that allow car manufacturers to demonstrate their know-how and stand out in a highly competitive automotive market. In a recent research report [START_REF]Advanced Driver Assistance System (ADAS) Market Size, Share, Growth, Trends, Regional Outlook and Forecasts[END_REF], the market research and consulting organization Precedence Research estimated the global ADAS market size at USD 32.3 billion in 2019 and forecasted that it would reach USD 142 billion by 2027.

In the perspective of equipping millions of cars worldwide, driver and pedestrian safety is the first concern of the ADAS developers. Therefore, implementing rigorous validation processes is necessary to verify the ADAS performances in all driving situation.

SAE J3016 TM LEVELS OF DRIVING AUTOMATION TM

DRAFT-Stand alone

Validation

The validation process objective is to assess the ADAS conformity to technical specifications, established by Groupe Renault in accordance with car safety evaluation agencies recommendations (e.g., Euro NCAP in Europe).

In this context, Groupe Renault has specified a minimum reliability level of 10 -8 casualties per hour to validate a LEVEL 4 ADAS system, which is 10 times less than the average rate of fatalities per hour on European highways, and 100 times less than the average accident rate on all types of European roads. It would require at least 2 × 10 10 kilometers of driving [START_REF] Didier | Autonomous Vehicles Digital Simulation Challenges[END_REF], equivalent to 100 years of road testing. This reliability level must be reached for every vehicles, every ADAS and every traffic regulation context. In addition, the ADAS developers must also ensure that the ADAS are inter-operable and that their combination does not provoke side effects. Given this complexity, it is impossible to test every scenario with only physical driving test on-track or open roads.

Groupe Renault has made the strategical choice to invest in massive driving simulation technology to enhance the ADAS development and validation process. The simulation systems mimic car driving conditions based on a driving scenario, vehicle physics, driver behavior, and interactions with a configurable environment. The simulation quality, i.e., its capacity to reproduce real-life driving conditions faithfully, is methodically monitored to ensure the validation process pertinence. The simulation system is illustrated in Fig. .8.

The simulation outputs large amounts of information coming from the observation of the simulated variables, and aggregated as one large dataset of multivariate time series. Based on this multivariate dataset, the role of the field expert's is to discriminate the driving behaviors, detect anomalies, and compare the ADAS performances to the specifications.

Data dimensions and complexity are considerable : for a given use case, the number of simulations can be as large as O(10 6 ), described O( 103 ) of simulated sensors, each recording O(10 4 ) time steps. Overall, more than O (10 13 ) data points are produced to test one ADAS setting, which represents several gigabytes per simulation plan (up to 1 terabyte). Because the simulation duration depends on the driving scenario, the simulated time series lengths are unequal. However, because the variables describe the simulation at the same sampling rate, the time series lengths are equal for a given simulation. Even though some simulated driving behaviors are predictable (e.g., the predicted valid behavior), most are discovered after simulation, which prevents from using labels to supervise the analysis. In this high-dimensional, multivariate and unsupervised context, extracting meaningful information is a tedious and time-consuming task for the expert, and requires specific exploration tools. The action of automatically regrouping time series into homogeneous clusters, without supervision, is called time series or clustering. An instance of multivariate time series clustering method applied to a driving simulation dataset is displayed in Fig.

. 10, where the scatter plot on the right shows the multivariate observations regrouped into clusters and represented in a simpler 2D space.

Contributions

The clustering of time series is a complex task that depends heavily on hypothesis that the user decides to form on the considered dataset. Without label, it is never possible to know if the clustering result is correct, because there is no ground truth to compare the results to. On the contrary, when an user applies a clustering method on a dataset, it is the understanding of the underlying assumptions that allows to reason and draw conclusions on the results.

The first contribution of this thesis is based on a scenario-construction hypothesis. After observing that simulations are generated on the basis of driving scenarios (c.f., Fig. .8), we developed a method based on the segmentation of time series in chains of driving phases. With this method, the time series are recoded as chains of 'regimes', where regimes are driving patterns found in the entire datasets and aggregated in a dictionary. As a probabilistic method, this approach can be used for outlier detection based on probabilistic confidence interval and missing data inference. We extend this approach with three new models :

-A Bayesian non-parametric time series coclustering model that natively integrates a model selection step.

-A parametric Time Series Multi-Clustering model that infers one column-partition and one specific row-partition per column cluster (c.f., Fig. .11-b).

-A non-parametric multi-coclustering. This approach combines the multi-clustering with coclustering structures, to construct several views of the same multivariate time series dataset while reducing the number of parameters and enhancing interpretation (c.f., Fig. .11-c).

The source code associated with these distributions, mainly Scala code, with visualizations in R and Python, is also made available, with the public datasets used for test.

The github repositories links are summarized in Tab. .2. The industrial use cases cannot, however, be made available.

Summary of remaining chapters

The rest of the manuscript is organized as follows :

Chapter 1 : Related work

The first chapter focuses on the existing works related to time series unsupervised classification, univariate or multivariate. This chapter also introduces the model-based clustering, coclustering, multi-clustering, and the Bayesian non-parametric modeling, that will intensively be used in the following chapters.

Chapter 5 : Multivariate Time Series Multi-Coclustering

Built on the methods detailed in Chapter 4 and Chapter 3, the multi-coclustering takes the best of coclustering and multi-clustering. This method consists in hierarchically nested Dirichlet Process Mixtures that infer a multi-clustering view of a dataset and then a coclustering structure in each view.

Chapter 6 : Conclusion And Perspectives

This chapter summarizes possible extensions for the proposed methods. These extensions focus on multivariate extensions for the dictionary-based univariate clustering method, and different representations in the block-clustering approach.

Publications

These contributions were the subject of several publications, chronologically represented in Fig. .12, and listed below :

-Goffinet, É., Lebbah, M., Azzag, H., Giraldi, L., & Renault, S. A. S. (2020). Clustering de séries temporelles par construction de dictionnaire. In EGC (pp. 181-192).

(z t ) T Segment membership ϕ = (ϕ q ) Q Basis of polynomial func- tions β = (β q ) Q Polynomial regression co- efficients R = {r u } U Dictionary of regimes (size U ) (w k,s ) K Logistic process parame- ters in component k σ 2
covariance parameter

(ρ h ) H Multiple coclustering column proportions v = (v j ) p Multi-Clustering variable cluster mem- bership η = (η h ) H Multi-Clustering com- ponent proportions θ h k,l
Multiple coclustering block distribution parameters

µ h k,l , κ h k,l , Ψ h k,l , ν h k,l Block cluster (k, l) pos- terior distribution pa- rameters γ, (α h ) H , (β h ) H Multiple coclustering concentration parame- ters a γ , b γ a β , b β a α , b α

Beta hyper-prior parameters

Related Work

1

" Knowledge is a process of piling up facts ; wisdom lies in their simplification.

-Martin H. Fischer

In the previous chapter we described the industrial context that motivated this thesis contributions : during the ADAS validation step, the field expert deals with large quantities of driving simulations, with unequal lengths and without labels. As a consequence, he must rely on unsupervised tools to explore the datasets contents.

Time Series Clustering

This chapter describes the existing works on the topic of time series clustering, starting from the Clustering basic definition in Sect. 1.1.2. The subsequent Sect. 1.1.3 and 1.1.4 respectively introduces a review of time series dissimilarities and transformations. Section 1.1.5 describes existing clustering algorithms, with an emphasize on the methods that are used and extended in the present thesis.

Time Series Analysis

The last decades has seen an important increase of the Information Technology field, due to decrease of processors, sensors and storage costs correlated with an increase of available data and use case development.

Time series is a specific kind of data generated from successive observations of a system state. This data type can be observed in every time-depending field. In Finance, they describe the sales, stock prices and exchange rates evolution. They are also used extensively in Health use cases, for the description of a patient state through medical sensors.

In Industry applications, they can monitor machine activity, production rates, and optimize predictive maintenance. Logistic, Demography, Retail, Climate, Astronomy,. . . The application possibilities are numerous.

In correlation with the use case number, processing power and storage capacities, the time series analysis topic has seen a sharp interest increase. Three kind of usage dominates the time series data applications : missing data inference (which includes forecasting when these data are posterior to the observed time series), anomaly detection, and classification.

In our applications, we are particularly interested in two of these usages : the detection of anomalous driving behaviors, and unsupervised classification of driving patterns.

Unsupervised Classification

An informal definition of unsupervised classification, or clustering, is : "(...) to discover groups of similar examples within the data (...)" [START_REF] Christopher | Pattern recognition[END_REF]. This partitioning is performed without the help of labels to supervise the estimation. The clustering transforms a potentially complex dataset (in our case, high-dimensional multivariate time series) into a simplified information, expressed as a categorical vector of memberships.

As opposed to supervised classification (called simply classification), that seeks to correctly infer an observation class based on a ground-truth, clustering is an exploration tool that helps the user understand the content of a dataset. It can be compared to a cartography tool, whose goal is not to recognize cities and countries (as classification does), but to detect new geographical patterns in an unexplored area.

Every clustering method is associated to a set of parameters. Some of these parameters are estimated during the clustering process, but others are specified by the user before the run. These are called hyper-parameters. For instances, some methods requires to know the cluster number a priori, to specify a time series dissimilarity measure, or an initial inference state.

These hyper-parameters may have an important impact on the clustering result. For instance, the user can choose to discriminate time series based on time series trends, shape, or simply length or sampling frequency. These features have nothing in common (two time series can have similar trends but completely different length). As a consequence, depending on user hypothesis, the results can change completely. It is up to the user to choose among all the possible representations, features, dissimilarities and every possible hyper-parameter settings, according to its need and prior knowledge. It is therefore mandatory for the user to fully understand how the clustering process works, and the role and impact of each hyper-parameter.

Many time series clustering methods have been developed in the last decades, based on various construction hypothesis and application domains [Lia05 ; Keo06 ; Din+08 ; Wan+13 ; ASW15 ; For21]. In the following we detail some of the existing approaches and strategies, with an emphasizing on the methods used in the present thesis.

Time series similarity measure

The broad definition of clustering ("grouping similar elements together") outlines the importance of the similarity definition in a clustering process. One way of clustering time series is to use time series specific distance. Clustering methods based on time series in the original space. This approach is also called raw-data-based, or shape-based approach (this denomination is, for instance, given in [START_REF] Jesin | Clustering time series using unsupervised-shapelets[END_REF] or in [START_REF] Saeed | Time-series clustering-a decade review[END_REF]).

The definition of a time series dissimilarity measure is not straightforward, because there are many different ways to discriminate time series. As a consequence, many distances or similarity measures have been proposed in the literature [START_REF] Félix | Analysis of similarity measures in times series clustering for the discovery of building energy patterns[END_REF].

In the univariate time series clustering context, let denote s = (s i ) n a set of time series observations, such that each element s i = (s i (t)) T i is a time series indexed by T i .

In the following we call time series dissimilarity measure d a function of two time series with the following properties, for two time series (s 1 , s 2 ) :

d(s 1 , s 2 ) ≥ 0, d(s 1 , s 2 ) = d(s 2 , s 1 ), d(s 1 , s 1 ) = 0.
The definition of distance adds the triangle inequality property, i.e., for all (s 1 , s 2 , s 3 ),

d(s 1 , s 2 ) ≤ d(s 1 , s 3 ) + d(s 3 , s 2 ).
The most popular time series distance is theℓ 2 -norm (also called Euclidean distance) that compares time series point-to-point. This distance is defined by : However, this distance is strongly impacted by time axis distortion and requires aligned time series with equal length.

d(s 1 , s 2 ) = T (s 1 (t) -s 2 (t))
The Pearson's correlation distance is used for comparing time series based on shape similarity. However, as stated by [START_REF] Félix | Analysis of similarity measures in times series clustering for the discovery of building energy patterns[END_REF], this distance "is sensitive towards outliers and generally not suitable for other probability distributions except for the Gaussian one", which makes it impractical in many situations.

The Mahalanobis distance is based on the Euclidean distance and integrates the correlation information to form a distribution-based distance [START_REF] Zoltán | Time series classification by class-specific Mahalanobis distance measures[END_REF].

The Dynamic Time Warping (DTW) [START_REF] Hiroaki | Dynamic-programming approach to continuous speech recognition[END_REF] In addition to the high (discrete) number of existing dissimilarity measures, it is also possible to use a weighted or nested combination of them [START_REF] Krisztian | Fusion of similarity measures for time series classification[END_REF], which transforms the similarity choice from a discrete into a continuous optimization problem. However, this approach seems more adapted to a supervised context where it is easier to estimate the distances combination weights.

Because time series data are high-dimensional objects, computing distances in the original space is often expansive, as is every step of the analysis process (storage, preprocessing, analysis, visualization, . . . ). In addition, in high-dimensional space, discriminating time series is complex because of the curse of dimensionality [START_REF] Michel | The curse of dimensionality in data mining and time series prediction[END_REF], an effect that tends to amplify the observations separation when the observation space dimension increases. A preprocessing step is often applied to represent time series into smaller-dimensional spaces, where this separation is easier.

Time Series Transformation

The feature-based approach consists in designing a meaningful condensed time series representations. This process assumes that the transformation process keeps the data informative value. Two situations can be distinguished : the first when the transformation process is known, the second when it is guided by an external criterion. Either way, this method requires prior knowledge on the time series.

Explicit transformation definition

This first category of methods (also called Non Data Adaptive) is used when the user can explicitly define the transformation that he finds relevant. The mean, median, quantiles, length, sampling frequency, first or last element's values, are instances of such features. These "single-point" features are usually easy to compute, and greatly reduce the dataset dimensions [START_REF] Teemu | Feature-based clustering for electricity use time series data[END_REF], but the feature selection step may require a long prospective work.

This class of methods also includes sub-sampling, or segmenting with fixed interval windows. The Piecewise Aggregate Approximation (PAA) [START_REF] Eamonn | Dimensionality reduction for fast similarity search in large time series databases[END_REF] is a popular representation based on an uniform time step interval segmentation, then compute the mean value on each segment.

On the basis of the PAA representation, the Symbolic Aggregate approXimation (SAX) [START_REF] Jessica | Experiencing SAX : a novel symbolic representation of time series[END_REF] is one of the first developed dictionary method that regroups segment mean The Symbolic Aggregate ApproXimation (SAX) method recodes the time series with a dictionary built with the Piecewise Aggregate Approximation (PAA) representation (here, 9 segments described with a four-symbol dictionary).

values belonging to the same gaussian quantile interval. This grouping produces a dictionary that can be used to recode the original time series in chains of symbols. An illustration of the method is presented in Fig. 1.2.

In the Bag-Of-Pattern (BOP) method, the same principle is applied to subsequences of a time series to form word (i.e., chains of symbols from SAX). The final BOP representation consists in the words frequencies. In the Bag-Of-SFA Symbols (BOSS) [START_REF] Patrick | The BOSS is concerned with time series classification in the presence of noise[END_REF] approach, the authors uses a truncated Discrete Fourier Transform (DFT) to represent the segment values, instead of the mean, and the dictionary is constructed by clustering the obtained Fourier coefficients. Other variants and extensions following the same scheme have been developed (WEASEL [START_REF] Patrick | Fast and accurate time series classification with weasel[END_REF], SAXSVM [START_REF] Pavel | Sax-vsm : Interpretable time series classification using sax and vector space model[END_REF], BoTSW [START_REF] Adeline | Bag-of-temporal-sift-words for time series classification[END_REF], . . . )

The dictionary representations reduce greatly the representation space dimensions, by segmenting both the representation space (from numeric to categorical) and the time index (from T i points to several segments). These dictionary-based approaches allow, to some extent, the comparison of similar segments between time series, but require some strong hypothesis on time interval sizes, and produces a coarsed segmentation of the time series phases, which is an issue for scenario-based time series clustering with asynchronous regimes. Short-Time Fourier Transform (STFT) (used for feature extraction in [START_REF] Lin | Multi-band multi-centroid clustering based permutation alignment for frequency-domain blind speech separation[END_REF]) applies the Fourier Transform on sliding windows sub-samples, and represents the frequencies evolution over time, which allows to obtain a representation in both frequency and time.

Closely related, the Discrete Wavelet Transform (DWT) representation applies a convolution with a wavelet function, which links window duration and frequency scale.

Another popular choices of representation consists in representing time series as weighted combinations of functions from a common functional basis (e.g., polynomial, spline, Fourier, Wavelets, . . . ). Combined with a Principal Component Analysis (PCA) [START_REF] Harold | Analysis of a complex of statistical variables into principal components[END_REF] step that performs a weighted linear dimension reduction on the regression coefficients, the overall transformation is known as the functional Principal Components Analysis (fPCA) [START_REF] James | Applied functional data analysis : methods and case studies[END_REF], which is a popular representation for model-based time series clustering [START_REF] Julien | Functional data clustering : a survey[END_REF].

Closely related to the PCA, the Multi-Dimensional Scaling (MDS) uses a singular value decomposition of the distance matrix to represent the time series [START_REF] António | Analysis of temperature time-series : Embedding dynamics into the MDS method[END_REF], and finds a low-dimensional representation that preserves the similarity measure in the original space.

When the functional regression basis is different for each time series, a common representation basis can be constructed to allow the time series comparison [START_REF] Jorge | Comparison of times series with unequal length in the frequency domain[END_REF].

Latent Representation Learning

Also called Data Adaptive transformations, some transformations are guided by an internal criterion optimization such as the reconstruction error or model likelihood.

In [START_REF] Jesin | Clustering time series using unsupervised-shapelets[END_REF], the authors uses theℓ 2 norm to extract the subsequences (called shapelets that best separates the time series, and use these subsequences as discriminant features.

For instance, with specific stationarity assumptions, time series can be represented with coefficients of AutoRegressive / Integrated / Moving Average models (ARIMA), that model the time series trends [START_REF] David | Discovering characteristic actions from on-body sensor data[END_REF], and/or with the AutoRegressive Conditional Heteroskedasticity (ARCH) [START_REF] Caiado | A GARCH-based method for clustering of financial time series : International stock markets evidence[END_REF] (or one of the many extended model, c.f. [START_REF] Tim | Glossary to arch (garch)[END_REF]) that model the time series variance evolution. These models parameters are inferred with ordinary least square or likelihood maximization.

The time series segmentation is another specific transformation that seeks to segment a time series according to relevant cutpoints (e.g., based on autoregressive model [START_REF] Takeuchi | A unifying framework for detecting outliers and change points from time series[END_REF] or piecewise regression [START_REF] Samé | Model-based clustering and segmentation of time series with changes in regime[END_REF]). It can be considered a variant of the PAA representation where the time steps are not uniform but estimated from the data. This allows a more accurate representation of the time series, and also reduces the final representation space.

Several works have also addressed the reconstruction-based feature extraction based on non-linear representation. For instance, [START_REF] Minh | m-tsne : A framework for visualizing high-dimensional multivariate time series[END_REF] Combined with a segmentation objective, the same principle is used by the Timeseries Forest [START_REF] Houtao | A time series forest for classification and feature extraction[END_REF]. The Self-Organizing Map (SOM), and its NN version DeepSOM [For+19], is another instance of nonlinear representation that expresses time series on a lower-dimensional grid such that observations that are close in the projection space are also close in the original space (with respect to the Euclidean distance).

Finally, inspired by the triangulation principle, a recent contribution [START_REF] Rohit | Using dynamic time warping distances as features for improved time series classification[END_REF] has shown the interest of using the distance to other elements as features.

After this transformation step, the user obtains a representation of the time series in a simpler space. In this space, the final clustering is often obtained with the Euclidean distance [START_REF] Saeed | Time-series clustering-a decade review[END_REF] associated to a clustering algorithm.

Clustering Algorithm

Whether in the original space or in a lower-dimensional representation, the clustering objective is to form groups of similar elements and discriminate dissimilar ones. The next subsections describe the existing types of time series clustering algorithms.

Density-Based

The density-based approach consider that a cluster is equivalent to high density area in the observation space (i.e., with an important concentration of observations), separated by areas with low density. By construction, these methods creates clusters based on density contiguousness and without assumption on clusters shape.

The DBSCAN [START_REF] Martin | A densitybased algorithm for discovering clusters in large spatial databases with noise[END_REF] algorithm attributes a state to each observations, among : core member, border member, or noise of a cluster, based on its neighborhood cardinality.

The clusters are then defined as the neighborhood of the core points. This algorithm assumes that the clusters share the same density and may fail to detect clusters with varying densities [START_REF] Ahmad | Performance Evaluation of Clustering Algorithm Using different dataset[END_REF]. In the time series analysis domain, DBSCAN has been applied in the original time series space (e.g., in combination with DTW, as in [START_REF] Haiyang | Integrated approach to density-based spatial clustering of applications with noise and dynamic time warping for breakout prediction in slab continuous casting[END_REF]), or in combination with a lower-dimensional representation (e.g., after an UMAP projection as in [START_REF] Clément | Improved Time-Series Clustering with UMAP dimension reduction method[END_REF]). Several extensions of this algorithm have been developed (see [START_REF] Loh | A survey on density-based clustering algorithms[END_REF] for a complete review).

The Density Peaks [RL14 ; Du+19] algorithm works in two steps : a) first it selects the cluster centroids based on local density and distance to others density peaks ; b) it assigns the remaining observations to the closest centers with higher density. This algorithm has been recently applied to time series clustering in the original space [Put+19 ; JLR20] but has also be seen used after a feature extraction transformation [START_REF] Ruitong | Fake comment detection based on time series and density peaks clustering[END_REF]. The Density Peaks algorithm has also been extended with TADPole [START_REF] Nurjahan | Accelerating dynamic time warping clustering with a novel admissible pruning strategy[END_REF] that combines the DP with the DTW measure [START_REF] Hiroaki | Dynamic programming algorithm optimization for spoken word recognition[END_REF].

Clustering with Neural Network

The growing development of deep learning since their successful applications in image recognition [START_REF] Alex | Imagenet classification with deep convolutional neural networks[END_REF], has given birth to NN-based solutions in every field of data analysis. It has been used for time series transformation, as presented in Sect. 1.1.4, but also in time series clustering.

Neural Network clustering is based on two elements : a neural network architecture and a clustering objective (called the loss). The architecture is composed of interconnected neurons, that act individually as parameterized functions : each neuron takes a value as input and outputs a transformed value. The neurons are disposed in layers, such that each layer is a weighted combination of the precedent. The inference process consists in estimating the neurons parameters and the layers combinations weights by optimizing the loss.

In the time series clustering domain, the loss is often a weighted combination of a reconstruction error objective and a clustering loss. For instance, the Deep Temporal Clustering (DTC) [START_REF] Sai | Deep temporal clustering : Fully unsupervised learning of time-domain features[END_REF] uses a an auto-encoder reconstruction loss and a clustering assignment hardening loss [START_REF] Junyuan | Unsupervised deep embedding for clustering analysis[END_REF]. The DeepSOM method [START_REF] Florent | Deep embedded SOM : joint representation learning and self-organization[END_REF], illustrated in Fig. 1.5, uses an auto-encoder reconstruction loss, in combination with a Self Organizing Map objective One drawback of the NN-based clustering is that the NN architecture is as a black-box for the user, equivalent to a complex and hidh-dimensional hyper-parameter set. It is complex to understand the impacts of the architecture setting on the time series representations. As a consequences, the user obtains a partition, but doesn't know how the time series are transformed during the clustering process, and how to improve the NN architecture if the results is irrelevant. In addition, the NN training is often computationally expensive and requires massive quantities of data.

Hierarchical-Based

A partition P = {P k } K over s is defined as a finite collection of non-empty and nonoverlapping clusters, such that the union of every clusters contains every elements of s. Formally :

K k=1 P i = s, h, k ∈ {1, . . . , K}, P k ̸ = ∅, P k P h = ∅.
We denote z = (z i ) n the membership vector, such that z i = k means that element s i belongs to P k , (equivalently called cluster k) and n k the number of elements belonging to cluster k.

The hierarchical clustering algorithms seek to create a hierarchy of nested clusters. This tree is composed of a root containing all elements, and each branch recursively splits the previous clusters in two, until reaching the leaves, i.e., when every cluster is a singleton. As a result, the user obtains a hierarchy of partition, that in turn must be pruned to produce the final partition. The pruning threshold (or equivalently the cluster number) is specified by the user. This clustering can be performed with an agglomerative or divisive procedure.

The agglomerative (also called bottom-up) approach starts from the leaves, and recursively merges the clusters up to the one-cluster root. This deterministic approach is fully described by two elements : a dissimilarity measure (that evaluate the element-toelement proximity), and a linkage criterion (that evaluates the cluster-to-cluster proximity).

At each cluster merge step, the linkage criterion is used to select the "best" cluster couple to merge. Among the many existing criterion, the most used in the time series clustering domain are the single-linkage [START_REF] John | Fast and accurate time-series clustering[END_REF], average-linkage [HT05 ; uc16] and Ward's criterion [JLR20 ; ZA18]. This last criterion (defined in [START_REF] Jr | Hierarchical grouping to optimize an objective function[END_REF]), that we also use in our first contribution (c.f. Chapt. 2), is used to select the cluster merge that maximizes the cluster separation and minimize the dispersion inside cluster. Formally, it measures the intra-cluster inertia, expressed as a mean of pairwise dissimilarities :

K k=1 1 2n k x i ,x j ∈p k d(x j , x j ) 2 .
It can be noted that this expression differs from the standard Ward's criterion expression [START_REF] Jr | Hierarchical grouping to optimize an objective function[END_REF], as it uses pairwise distances inside clusters instead of mean distance to the cluster center. We use this pairwise-distance expression, presented by [START_REF] Marie | ClustGeo : an R package for hierarchical clustering with spatial constraints[END_REF] and advocated by [START_REF] Nathanaël | Applicability and interpretability of Wards hierarchical agglomerative clustering with or without contiguity constraints[END_REF], that extends the inertia definition to non-Euclidean dissimilarity measures, as is the case in our application (c.f. Chapt. 2). The main interest of this definition is that it does not require to estimate cluster center, which can be troublesome for some dissimilarity measures (e.g., the average DTW sequence computation is an NPhard problem).

With this hierarchical clustering approach, the user obtains a tree partition representation such that, for every value of cluster number, the obtained partition minimizes the global intra-cluster inertia. The agglomerative approach is described in Algorithm 1.

As opposed to the agglomerative approach, the divisive (also called top-down) hierarchical clustering starts from the tree's root, containing all elements, and recursively split the clusters to obtain the final singleton clusters. The complexity of this approach is large Algorithm 1 : Hierarchical Agglomerative Clustering. Input : a dissimilarity measure and a linkage criterion ℓ Set the initial partition with singleton clusters : P (0) = {{s 1 }, . . . , {s n }} for m ← 1 to n do Select the best candidate merge : (p * , q * ) = arg min p,q∈P (m) ℓ(p, q) Remove p * and q * from P (m) P (m+1) = {P (m) , p * q * } end (2 n ), due to the important number of splitting possibilities. The splitting approach is usually performed with an heuristic to reduce this complexity, for instance a partitionbased algorithm (e.g., as in [START_REF] Michael | A comparison of document clustering techniques[END_REF]). The hierarchical-based clustering is illustrated on Fig. 1.6.

Partition-based

The partition-based methods seek to estimate a membership set z, based on two inputs : a fixed number of clusters K and a similarity measure d over the observation space. This category of algorithm includes the most used and known clustering algorithm : the k-means [START_REF] James | Some methods for classification and analysis of multivariate observations[END_REF]. Starting from a random membership state, the k-means algorithm alternates two steps : a) compute {µ k } K the centers of each cluster ; b) associate each In combination with time series specific distances (e.g., with the correlation measure [START_REF] John | k-shape : Efficient and accurate clustering of time series[END_REF], Mahalanobis distance [Nel12], . . . ), the k-means has been used extensively and can be considered a "classical" approach. Its combination with DTW, where the average sequence estimation is an NP-hard problem, has been addressed by an important amount of papers [Gup+96 ; NR07 ; NR09 ; PKG11 ; AT15 ; SDG16].

Other partition-based approaches are based on the Partitioning-Around-Medoid [KR09 ; NH02] (PAM) method, which consists in using the best candidate among the observed time series set s instead of the cluster average. This approach has also been considered as a work-around for the DTW averaging problem.

The k-means algorithm can also been applied after time series transformations (c.f. Sec 1.1.4), for instance in association with shapelet representation [START_REF] Jesin | Clustering time series using unsupervised-shapelets[END_REF].

Fuzzy clustering (also called also called soft-clustering) is another approach to clustering in which observations can be assigned to several clusters. This method can also be combined with different types of distances, and has been applied to time series clustering in [START_REF] Dat | Fuzzy c-means clustering-based speaker verification[END_REF][START_REF] Yongli | Incremental fuzzy C medoids clustering of time series data using dynamic time warping distance[END_REF]. Fuzzy clustering can be considered a distance-oriented interpreta-

Model-Based Clustering

The model-based clustering approach infers a generative model over the observation space by associating a latent distribution to each cluster. In the same way partitionbased algorithms assume that a given dissimilarity measure is a suitable discriminator, model-based methods assume that a given distribution is a suitable cluster generator. An illustration is shown in Fig. 1.8.

The model-based approach can be applied in the original time series observations space, in that case the whole time series is represented. This approach depends on hypothesis on the data structure. In [XY02 ; FK08], the authors use a mixture of auto-regressive models to cluster time series generated from auto-correlated processes. Several modelbased works have also addressed the clustering based on time series shapes, for instances [START_REF] John | Probabilistic curve-aligned clustering and prediction with regression mixture models[END_REF] that uses a mixture of spline regression, or [GS99 ; Gaf04 ; BJ05 ; Cha16] that use polynomial regressions.

In some situations, it is relevant to assume that the time series are the results of chains of different states. In that case, it can be relevant to assume a piecewise-distribution modeling. In some cases, the state transition time is shared by every time series (i.e., the transitions cutpoints are synchronous), in some others it can be shared among time series in the same cluster. This last case has been addressed by the Hidden Markov Model (HMM) [START_REF] Padhraic | Clustering sequences with hidden Markov models[END_REF] and the piecewise regression mixture [Sam+11 ; CN18].

The model-based clustering can also be applied after time series transformations. The fPCA (c.f., Sect. 1.1.4) is a popular choice of representation in this context [START_REF] Charles | Model-based clustering of time series in group-specific functional subspaces[END_REF][START_REF] Julien | Funclust : A curves clustering method using functional random variables density approximation[END_REF].

The Bayesian model framework allows to add prior information on the mixture parameters under the form of prior distribution. This addition allows to automatically infer the model dimension during the inference process. This framework is used for instance in [START_REF] Marc | Functional data clustering via piecewise constant nonparametric density estimation[END_REF] to perform univariate time series clustering based on a grid-model.

As stated in this section preamble (c.f. Sect. 1.1.1), our use case specifications also include the detection of anomalous driving behaviors. It is one of the reason that motivates the choice of the model-based approach, that natively includes a probabilistic anomaly detection feature.

The model-based clustering algorithms assume the presence of a latent mixture model, where each component generates a group of elements.

Mixture Model Framework 1.2.1 Mixture Model

Definition

Let denote x = (x i ) n the set of vectors representing the time series after transformation in a d-dimensional real space (c.f., Sect. 1.1.4). We emphasize that every observation x i share the same dimension d. The variable x designates a random vector following the mixture density, and associated to a membership z.

Mixture-Modeling is a standard model-based clustering approach [START_REF] Geoffrey | Finite mixture models[END_REF]. This model assumes that the overall density on the d-dimensional space is a convex combination of densities belonging to the same distribution family F . The weights of this combination are called the mixture proportions (π k ) K , and the weighted densities are the mixture components, (F k (.)) K , with respective parameters (θ k ) K .

In the following we consider the specific case of F being the Multivariate Gaussian distribution family, as is the case in many of our applications. In this case, the components parameters are therefore

θ k = {µ k , Σ k }.
This probabilistic approach provides probabilistic noise modeling, confidence intervals and probabilistic outlier detection. With the notation defined in the introduction of this section, the Mixture Model (MM) density is given by p

(x) = K k=1 π k F (θ k ), with π k = p(z = k),

Inference

With Θ the set of all parameters, the model log-likelihood of the dataset x is given by :

ln p(x | Θ) = ln Z p(x, z | Θ), (1.1)
where Z is the set of all possible combination of the membership vector z. A direct likelihood optimization estimation would require to set the derivatives of Eq. (1.1) to zero, but the summation on Z prevents from obtaining a closed form solution [START_REF] Christopher | Pattern recognition[END_REF]. Several methods have been developed to sort out this limitation, of which the Expectation Maximization [DLR77] (EM) is the most popular.

The EM is an iterative inference algorithm, that repeats two steps until likelihood convergence.

p(z = k | x) = p(z) p(x | z) p(x) = π k p(x | z = k) K π h p(x | z = h) = π k F k (x) K π h F h (x)
.

In the following, we denote τ i,k = p(z i = k | x i ) for brevity.

M step estimating Θ by Maximizing the likelihood with fixed values of the posterior distributions.

Considering fixed values of (τ i,k ) n×K for every membership z of the dataset elements x, the parameter set Θ update is given by maximizing the expectation of the complete log-likelihood [START_REF] Christopher | Pattern recognition[END_REF] :

Q(θ old , θ) = n,K p(z i = k | x i , θ old ) ln p(x i , z i | θ),
with θ old the parameter values used for the estimation of the posterior distribution in the previous step E. The optimization of this expectation yields the following expressions [START_REF] Christopher | Pattern recognition[END_REF] in the multivariate Gaussian case :

n k = i τ i,k , π k = n k n , µ k = 1 n k i τ i,k x i , Σ k = 1 n k i τ i,k (x i -µ k )(x i -µ k ) T .
Because the posterior distribution p(z | xθ) can be computed individually, and the M step can be estimated by batch, the EM algorithm is easy to parallelize and to distribute on map-reduce frameworks [ @Spa].

With this construction, the model log-likelihood necessarily increases at each iteration [START_REF] Christopher | Pattern recognition[END_REF], or stagnates in case of convergence. However, the obtained solution is sensible to the EM initialization, and the EM can get stuck in locally optimal solutions. This sensitivity to initialization can be reduced with specific strategies [BCG03 ; BB13 ; BC15 ; SSB17], including random search, initialization with k-means, or hierarchical divisive initialization strategies. Keeping the best result among several tries is also a popular strategy.

The Stochastic EM (SEM) is a variant of the EM that adds a membership sampling step based on the posterior distributions estimated in the E step. With this version, the likelihood monotonic increase property is not conserved, but the stochastic step allows a better exploration of the solution space and reduces the risk of getting stuck in local optima [CD85 ; CD92].

The Classification EM (CEM) algorithm aims to optimize the complete likelihood p(x, z | Θ) instead of the marginal likelihood p(x | Θ). During the CEM, the τ i,k are used to select the most probable partition (as opposed to the SEM sampling). This strategy can be considered as hard clustering.

It can be noted that the CEM and the k-means are tightly connected. The k-means is in fact strictly equivalent to the CEM in the multivariate Gaussian case, when assuming identity covariances matrices (i.e., independent variables, and shared variance among cluster and variable). From this perspective, the EM can be considered a probabilistic soft-clustering generalization of the k-means. The EM, SEM and CEM are described in Algorithm 2.

Θ (0) for m ← 1 to M do for i ← 1 to n do Estimate τ i,k = p(z (m+1) i = k | x i ) given Θ (m) τ ′ i,k ← τ i,k if SEM then Draw membership z i ∼ M ult((τ i,k ) K ) Set τ ′ i,z i to 1, the others to 0 if CEM then Choose z i = arg max K τ ′ i,k Set τ ′ i,z i to 1, the others to 0 Estimate Θ (m+1) given (τ ′ i,k ) n×K Estimate ẑ.
As is the case with k-means, a model selection step is often required in association to the EM-based algorithms applications.

Model Selection

The EM requires to make several assumptions : on the number of clusters K (as the kmeans) but also on the component distribution family F . In the Multivariate Gaussian case for instance, the clusters can be assumed with equal covariance matrices (homoscedasticity), with diagonal covariance matrices (variable independency), . . . . These assumptions can be a strong limitation to real-life applications, where the number of clusters or the covariance matrices properties are rarely known.

To circumvent these limitations, model selection strategies have been developed, on the basis of selection criteria and search strategy.

Model Selection Criterion

In the parametric model-based framework, several model selection criteria have been developed. Most of them consists in a score that penalizes the model likelihood by its complexity (e.g., the Akaike Information Criterion (AIC) [START_REF] Akaike | A new look at the statistical model identification[END_REF], the Bayesian Information Criterion (BIC) [START_REF] Schwarz | Estimating the dimension of a model[END_REF], the Integrated Completed Likelihood (ICL)) in an attempt to approximate the integrated likelihood [START_REF] Christophe | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]). Under the hood, this approach aims at optimizing the bias-variance trade-off, and reducing the model over-fitting. Other non model-based approach have also been developed, as the recent Stadion [For+21] that measures the time series clustering stability against dataset transformation.

We chose the BIC in our parametric application model selection (Chapt. 2), as it has been shown in [START_REF] Christophe | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] that it well suited for density estimation tasks. Its definition is :

BIC = ln p(x | Θ) - C 2 ln(n) -2,
with C the parameter number. However, the BIC's drawback is that it tends to overestimate the number of clusters in order to obtain better density estimation, which may happen when the mixture components are not adapted to the true cluster shape. This over-estimation phenomenon has been taken into account in our application, and several measures have been implemented to mitigate the problem (e.g, Sect. 2.2.4 and Sect. 2.2.4).

Search Strategies Several search strategies exist in the parametric setting for selecting the best model structure.

The grid-search (also called exhaustive search) tries every model possibility and select the best, i.e., the model optimizing the model selection criterion value. The hierarchical search is, in fact, the equivalent of a divisive hierarchical algorithm (c.f. Sect. 1.1.5), based on a model selection criterion. The random search [START_REF] James | Random search for hyper-parameter optimization[END_REF] is another popular hyper-parameter optimization strategy, that consists in exploring the solution space randomly.

These model selection strategies make several assumptions. For instance, the grid search assumes that the true number of clusters is within the grid. The greedy search is a divisive hierarchical heuristic, and assumes that the clusters construction is based on a binary tree structure.

In some situations, the parametric model selection becomes problematic. For instance, when the number of models is extremely important.

Another solution to mixture model selection consists in estimating the model structure during the inference process. This can be done by taking the Bayesian perspective and assuming prior distributions over the model parameters, as presented in the following section.

Bayesian Non-Parametric Modeling 1.3.1 Alternative MM Representation

With the introduction of the discrete distribution G = K k=1 π k δ θ k , with δ θ the Dirac delta function, the MM admits the alternative equivalent representation :

x i | θ i ∼ F (θ i ) , θ i ∼ G, i ∈ {1, .., n}.
In this setting, each observation x i is drawn from the distribution F with specific θ i , and each θ i drawn i.i.d from distribution G. G is a linear combination of Dirac delta functions, weighted by the proportions π, and acts as a distribution over θ i , i.e., G is a distribution over distributions. An instance of G, with K = 4 is shown in Fig. 1.10, in the case where (θ k ) K are modes of bivariate normal distributions.

Because G is a finite discrete distribution and the number of cluster K is (supposedly) specified lower than n, some groups of observations share a common parameter, which

Dirichlet Process Mixture Model

The DPMM (c.f. Fig. 1.11) is a Bayesian Non-Parametric (BNP) model that assumes the presence of an infinite discrete number of latent components, of which only a fraction is observable in the dataset. In this setting, G is given an additional Dirichlet Process (DP) prior distribution that takes two hyper-parameters : a concentration α and a base

x i | θ i ∼ F (θ i ) , i ∈ {1, .., n}, θ i | G ∼ G, G ∼ DP (α, G 0 ).
Two alternative representations are often use to help interpreting the DP : the Chinese Restaurant Process (CRP), which is a distribution over partitions, and the Stick-Breaking (SB) scheme, that models the proportions construction. These two representations are strongly connected, as one can be derived from the other [START_REF] Jeffrey | An elementary derivation of the Chinese restaurant process from Sethuramans stick-breaking process[END_REF].

Chinese Restaurant Process

The CRP is a distribution over partition that depends on a parameter concentration α. The process is named after the culinar metaphor that illustrates its definition : a chinese restaurant contains several tables (each

p(z i | z -i , α) ∝      n k n -1 + α existing cluster/table k, (1.2) α n -1 + α new cluster. (1.3)
With this definition, the density of a partition P containing K clusters with cardinality n k can be obtained iteratively, with the product of conditional memberships probabilities :

p(z) = p(z 1 | z 2 , . . . , z n )p(z 2 | z 3 , . . . , z n ) . . . p(z n ) = α K i (i -1 + α) k (n k -1)!. (1.4)
The second line expression is obtained by noting that :

-Because there is K clusters in P , the "new cluster" discovery option has been chosen exactly K times, which explains the α K numerator.

-The denominator 1 n-1+α is present in both membership probabilities Eq. (1.2) and Eq. (1.3), which yields the product i (i -1 + α).

-For each cluster k, the final cardinality n k implies that the choice "entering existing cluster k" has been made exactly n k times, with n k varying from 1 (for the second customer choosing table k) to n k -1 (the last time a customer entered this table). This yields the term

1 × 2 × .. × n k -1 = (n k -1)!.
As shown by [START_REF] Charles | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF], the expression in Eq. (1.4) corresponds to the density of a partition in a DP process. Moreover, it can be noted that this density is invariant to the order of the customer arrival. This makes the With this construction in mind, the DP can be seen as a CRP prior on the memberships that also integrates the predictive distribution information.

Stick-Breaking

The stick-breaking process is another representation that focuses on the construction of the mixture proportions (π k ) ∞ . This construction comes from [START_REF] Jayaram | A constructive definition of Dirichlet priors[END_REF] which states that, if G ∼ DP (α, G 0 ), then it can be written as :

G = k π k δ θ k , 1.3 Bayesian Non-Parametric Modeling
the component distributions θ k follows the prior G 0 , and the weights (π k ) ∞ are obtained with the following stick-breaking construction : With this construction, [START_REF] Jayaram | A constructive definition of Dirichlet priors[END_REF] proved that the obtained distribution

π i (v) = v i i-1 j=1 (1 -v j ) , v = (v i ) {1,...,n} , v j i.i.d. ∼ Beta(1, α). ( 1
G = ∞ k=1 π k (v)δ θ k follows a DP(α,G 0 ) distribution.
The DPMM can be seen as a BNP extension of the MM where the number of components is infinite. Because the observation set is finite, only a finite number of clusters (noted K obs ) are observable. Moreover, the number of observed cluster grows in fact as O(log n) [START_REF] Whye | Hierarchical dirichlet processes[END_REF], which creates a clustering. This property is illustrated on Fig. 1.12.

It it possible to influence the value of K obs by tuning the concentration parameter α, which influences the probability of new cluster discovery. The role of α is also visible in the expression of K obs posterior distribution [START_REF] Charles | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF] :

p(K obs = K | α) = |S n,K |n!α K Γ(α) Γ(α + n) , ( 1.6)
where |S n,K | is the unsigned Stirling number of the first kind.

The DPMM can automatically adapt in case of dataset augmentation by adjusting the observable component number, which makes it relevant for exploration tasks where the user can use additional resources to explore new observation space areas.

Inference

Several methods have been developed to infer DPMM's parameters. For instance, based on variational inference [START_REF] David | Variational inference for Dirichlet process mixtures[END_REF] or Markov chain Monte Carlo (MCMC) [START_REF] Michael | Estimating normal means with a Dirichlet process prior[END_REF][START_REF] Radford | Markov chain sampling methods for Dirichlet process mixture models[END_REF]. For large datasets applications, variational inference methods have often been preferred over MCMC for their speed, at the cost of hypothesis on the posterior distribution structure. However, recent works [WDX13 ; Meg+19] have made MCMC processes scalable, which rehabilitates their use for industrial purposes.

In the following we present two MCMC algorithms, based on the Gibbs Sampler estimation. The Gibbs Sampler objective is to mimic the posterior distribution p(z | x, χ) by sampling values for each z i from its conditional distribution given the data and the remaining observation p(z i | x, z -i , χ).

Gibbs Sampler

The Gibbs Sampler, corresponding to algorithm 2 of [START_REF] Radford | Markov chain sampling methods for Dirichlet process mixture models[END_REF], consists in alternating the update of the clusters memberships and the clusters parameters. In a first step, the memberships z are updated one by one. Each membership z i is sampled from the conditional distribution given the observations x, the remaining memberships and the associated parameters : θ -i = (θ j ) j̸ =i , according to :

p(z i | z -i , x i , χ) ∝      n k n -1 + α p(x i | z i = k), existing cluster k, (1.7) α n -1 + α p(x i | χ), new cluster, (1.8)
where p( In Eq. (1.8), the density p(x i | χ) is the prior predictive distribution of x i , that can be obtained by marginalizing over the component parameter :

x i | z i = k) = F (x i , θ k ) is
p(x i | χ) = θ p(x i , θ | G 0 ) dθ = θ p(x i | θ)p(θ | G 0 ) dθ.
This integral has a closed-form when F is conjugate to G 0 . With F the multivariate Gaussian density and G 0 the Normal Inverse Wishart (NIW) conjugate prior with hyperparameter (µ 0 , κ 0 , Σ 0 , ν 0 ), this integral sums up to a multivariate t-distribution density (c.f. derivations in Appendix A.1) :

t ν 0 -d+1 x|µ 0 , (κ 0 + 1)Ψ 0 κ 0 (ν 0 -d + 1)
.

After this first step, i.e., after having updated every membership once, the clusters parameters are also updated by sampling according to the posterior distribution p(θ k | x k , G 0 ). This distribution is also a NIW distributions, with hyper-parameters updated with x k = {x i : z i = k} the observations belonging to cluster k. The hyper-parameters update equations are given [START_REF] Kevin | Conjugate Bayesian analysis of the Gaussian distribution[END_REF] by :

κ k = κ 0 + n k , ν k = ν 0 +n k , µ k = κ 0 µ 0 + x k x κ n , ( 1.9)
Σ k = Σ 0 + x k xx T + κ 0 µ 0 µ T 0 -κ k µ k µ T k .
(1.10)

These two steps (membership update and parameter update) are repeated until convergence or for a given iterations number.

At each iteration m, n memberships are updated, and each membership update involves the computation of K (m) existing clusters membership (Eq. (1.7)) and one new cluster memberships (Eq. (1.8)). The prior predictive distribution in Eq. (1.7) is fixed given the observations and prior G 0 , and can be computed once before the inference and cached. This first computation has a complexity O(nd 3 ), due to the multivariate t-distribution estimation cost.

The K existing membership probabilities must be computed again at each iteration, but this computation is quite cheap as it only consists in a multivariate Gaussian density estimation (O(d 2 )). The membership update step has therefore a O(n Kd 2 ) complexity, with K = max M m=1 K (m) an upper bound of the cluster number, and M the maximal iteration number.

During the parameter update step, each observation contributes to the update of the NIW posterior distribution associated with its current membership. This update step complexity is bounded by the covariance estimation term, in O(nd 2 ). The global complexity is therefore :

O(nd 3 + M n Kd 2 ).
The sequential aspect of the inference may seem a prohibitive feature at first sight, as it prevents the membership update parallelization. However, some works have addressed the problem [WDX13 ; Meg+19], and turned this algorithm into distributed and scalable versions.

One drawback of this inference process is that the clusters parameters are updated only once per iteration, i.e., after the update of every memberships. Because of this low up-date rate, the convergence can be slow. In the next algorithm, the cluster memberships state is taken into account at all time.

Collapsed Gibbs Sampler

The Collapsed Gibbs Sampler, corresponding to algorithm 3 of [START_REF] Radford | Markov chain sampling methods for Dirichlet process mixture models[END_REF], is a natively fast MCMC's method that also assumes the prior distribution G 0 conjugate to the density family F and skips the parameter sampling step. In this method, the prior conjugacy assumption enables the closed-form computation of the posterior predictive distributions, that is used to estimate posterior cluster membership probabilities. The density F in Eq. (1.7) is replaced by this posterior distribution :

p(z i = k | x i , x k , G 0 ) = θ p(x i | θ)p(θ | x k , G 0 )dθ.
(1.11)

In this expression, p(θ | x k , G 0 ) is the same distribution than the one used for parameter update in the (non-collapsed) Gibbs Sampler from Sect. (1.3.3), i.e., a NIW distribution with hyper-parameters obtained by updating the prior G 0 parameters with Eq. (1.9) and Eq. (1.10).

With these updated hyper-parameters, the integral Eq. (1.11) is equivalent to the following multivariate t-distribution :

t νn-d+1 x|µ n , (κ n + 1)Σ n κ n (ν n -d + 1)
.

This direct computation of the posterior membership distribution allows to collapse the clusters parameters, and to avoid their estimation during inference.

At each iteration, the update of one row-membership z i involves the computations of (K (m) + 1) row-cluster membership probabilities (K (m) existing clusters and the new cluster). For each cluster k, the membership probability computation involves the cluster hyper-parameter update (based on Eq. (1.9) and Eq. (1.10)) and one posterior predictive distribution density estimation.

In practice, instead of updating the cluster k posterior hyper parameters (µ k , κ k , Ψ k , ν k ) with the full cluster contents x k at each iteration and at each membership update, it is more efficient to cache these hyper-parameters, and only remove the observation x i from its previous cluster before the membership probabilities computation, and add it to its new cluster after sampling, which is simply O(d 2 ). The posterior membership density estimation has a complexity O(d 3 ) because of the multivariate t-distribution computation. This membership must be performed once for each observation and iteration, which gives a total complexity of

O(M n(d 2 + Kd 3 ))
This complexity is comparable to the one of the Gibbs Sampler, and highlights the interest of keeping a representation space dimension low. Both the Gibbs Sampler and Collapsed Gibbs Sampler are MCMC that produces samples that mimics the joint distribution p(x, z | χ). These samples must in turn be processed to obtain the final partition.

Inferring the final partition

These Gibbs Sampler outputs a chain of partitions, that correspond to samples from an approximation of the posterior distribution. These sampled must, in turn, be aggregated over the iterations, usually after a given number of burnin iterations. Given a chain of partitions (ẑ m ) M , the objective is to estimate

ẑ = arg min Z m d(ẑ m , z),
with d a dissimilarity measure between permutations and Z the set of all possible rowpartition with size n. This consensus partition estimation is an NP-complete problem [START_REF] Mirko | NP-hard problems in hierarchical-tree clustering[END_REF], that several works have addressed in the past (c.f. [START_REF] Petros | A review on consensus clustering methods[END_REF] for an extensive review).

This question is crucial for real-life applications, where the inference process can be sensitive to initialization and produce different solution. One objective of the consensus partition is to obtain a "better" partition, in the sense that it is more stable and relevant. The consensus partition also has the advantages to resolve the label switching problem.

We advocate to use the recent method from [START_REF] Thomas | Finding the mean in a partition distribution[END_REF], which proposes an efficient extension of the combinatorial optimization method from [START_REF] John | Inference of population structure under a Dirichlet process model[END_REF] that constructs the partition with minimal partition-distance [START_REF] Dmitry | Partition-distance via the assignment problem[END_REF] to the samples, without assumptions on the final number of clusters or on the clustering structure. The complete inference process is described in Algorithm 3. Add x i to its new cluster.

Compute the average partition ẑ on the last samples (c.f. averaging methods in Sect. 1.3.3).

The methods presented in this previous section are dedicated to univariate time series clusters. When several temporal variables are present in a dataset, additional difficulties appear, that must be addresses with specific methods.

Multivariate Time Series Model-Based Clustering

As stated in this section preamble (c.f., Sect. 1.1.1), the time series clustering is applied in many domains and contexts. In most of the use cases, the time series can be described by several variables simultaneously, i.e., the time series datasets are multivariate. For instance, in Health, the patient state is often observed with several sensors (e.g., temperature, hearth rate, blood pressure, toxin concentration). In finance, the analyst is often interested in understanding the correlation between the stock price of several items.

In our industrial applications, the ADAS numerical validation also generates an important number of variables (position, angle, activation systems, radar, engine, . . . ). This Multivariate Time Series (MTS) context is a multi-dimensional extension of the univariate case, with additional difficulties that we present in the next sections.

Notation

As a natural extension of the univariate case notation, let denote S = (s i,j ) n×p the matrix of time series observations, such that each cell s i,j = (s i,j (t)) T i is a time series indexed by T i . This notation highlights one key property of this dataset : the time series length T i is shared among simulation (i.e., matrix row) and not among variable (i.e., matrix column). From the ADAS perspective, it comes from the fact that the simulated variables observe a simulation with a common duration. An instance of MTS dataset is given in Fig. Let denote (X i,j ) n×p the dataset transformed in a lower dimension space, where each cell share the same dimension d. In line with this matrix notation, we denote (abusively and for brevity) x i,. the i-th row/simulation of X, and x .,j its j-th column/variable.

Multivariate Context Constraints

In the multivariate setting, the observation space dimension is the product of all the individual temporal variable dimension, which aggravates the separation effect of the dimensionality curse. Even after individual time series transformations (which corresponds to the analysis of dataset X where each x i,j is a d-dimensional vector), the representation space has an d × p dimension. In our application, p can be as large as several hundreds, which makes the transformation space dimension prohibitively important even for small values of d.

Multivariate clustering must therefore include one or several dimension reduction features to address this dimension increase problem.

In the presence of several variables, the methods must also take the variables relationships into account. For instance, the Gaussian MM model use specific hypothesis on the covariance matrices to model the variables correlation.

In the following we describe some existing multivariate clustering methods, before introducing the block-clustering framework.

Related Work on Multivariate Time Series Clustering

Existing MTS methods are the direct result of an existing univariate method extension.

Multivariate Similarity Measure

An important number of papers have been focused on the extension of DTW [START_REF] Mohammad | On the non-trivial generalization of dynamic time warping to the multi-dimensional case[END_REF][START_REF] Mohammad | Generalizing DTW to the multi-dimensional case requires an adaptive approach[END_REF]. In [START_REF] Mohammad | Generalizing DTW to the multi-dimensional case requires an adaptive approach[END_REF], the authors discriminates the possible DTW extensions in two categories :

-Independent strategies (also called Match-By-Dimension), that compute the DTW dissimilarity on each variable independently, and use the sum as global dissimilarity score.

-Dependent strategies (also called overall Matching), that sums up the multidimensional DTW dissimilarity step by step. One main advantage of this approach is that the time axis warping is shared among dimensions, which conserves the time dependency information in the dissimilarity evaluation.

In a recent work, [START_REF] Ahmed | Elastic Similarity Measures for Multivariate Time Series Classification[END_REF] uses this independent/dependent categorization as a basis for a broad similarity extension framework, and applies it to several existing measures (L2 Norm, Longest Common Subsequence, Move-Split-Merge, . . . ).

At the crossover between feature extraction and MTS similarity, several MTS similarity measures have been defined on the basis of Singular Value Decomposition (SVD). It includes :

-Weighted Sum SVD (WSSVD) [START_REF] Cyrus | Real-time Pattern Isolation and Recognition Over Immersive Sensor Data Streams[END_REF], which consists in computing the scalar product of the MTS matrices eigenvectors, weighted by the eigenvalues.

-PCA similarity factor (S P CA ) [SS05] considers the independent PCA projection of two MTS, keeping only the first d axes, and evaluates the difference between every PCA axes of the first projection with every PCA axes of the second, as the squared value of their angle cosine, weighted by the eigenvalues -Extended Frobenius norm (Eros) [START_REF] Kiyoung | A PCA-based similarity measure for multivariate time series[END_REF] is close to S P CA , but considers the absolute value of the angle cosines, and considers every axes of the PCA representations.

Multivariate Transformation

The feature extraction strategy (c.f., Sect 1.1.4), which consists in representing time series with discriminate characteristic, has also been an important research field in the last years.

In the similar way than the MTS similarity measure can be categorized in dependent or independent extensions, the transformation can lead to an independent transformation (each dimension is individually transformed), or dependent transformation (the overall MTS is transformed). Note that the strategy choice has not implications on the final dataset dimensions, as independent transformations can lead to vectorized representation, and dependent transformation may produce multi-dimensional representations.

Few works have addressed the unsupervised transformation, compared to the univariate case. For instance, [DSP05 ; WWW07] that perform single-point feature (mean, deviation, quantile, . . . ) extraction.

The supervised feature extraction methods have received more attention. The vector linear autoregression [START_REF] Hiro | Vector autoregression and causality : a theoretical overview and simulation study[END_REF] method uses the autocorrelation to represent the dynamical information in the MTS. For instance, [START_REF] Pierpaolo | Wavelet-based self-organizing maps for classifying multivariate time series[END_REF] uses a multivariate wavelet transform and uses this representation as basis for the training of a SOM.

The PCA projection has also been extensively used for MTS representation. In [START_REF] Zheng-Xin | Dimension reduction method for multivariate time series based on common principal component[END_REF], the authors use the common principal component analysis [START_REF] Hailin | Multivariate time series clustering based on common principal component analysis[END_REF] representation, which consists in estimating a covariance matrix on each MTS independently, perform the PCA on the basis of the covariance matrices average, and projects the MTS in this common representation space. An important number of neural-network representations have also been recently (an extensive review is available in [START_REF] Ismail | Deep learning for time series classification : a review[END_REF]). For instance, [II20] uses a deep embedding based on a recurrent auto-encoder, or [Kip+18 ; PRB19] that models the interaction between MTS with a latent interaction graph.

As extensions of the univariate u-shapelet [START_REF] Jesin | Clustering time series using unsupervised-shapelets[END_REF], multivariate u-shapelets methods have been developed [Spe+18 ; OW20], based on dependent or independent extensions of the DTW.

Another approach consists in extracting relevant subsequences with a multivariate segmentation (e.g. with a model-based joint segmentation [DTS07 ; Cha+13 ; HNB19]).

Multivariate Time Series Clustering Algorithm

As in the univariate case, MTS clustering algorithms can be based on multidimensional models, similarity measures, in the original space or in a lower-dimensional representation. The clustering algorithms can be categorized in the same taxonomy than the univariate methods :

-The most popular partition-based algorithm, k-means, has been applied in [START_REF] Tamraparni | Grouping multivariate time series : A case study[END_REF] on single-points features, and in [Lee+20 ; SRV21] with multivariate DTW in the original space. In [START_REF] Dino | Deep multivariate time series embedding clustering via attentive-gated autoencoder[END_REF], the k-means is applied on a non-linear representation based on auto-encoder representation.

-In [START_REF] Chandrakala | A density based method for multivariate time series clustering in kernel feature space[END_REF], the authors have extended the density-based algorithm DBSCAN to handle MTS.

-The Hierarchical clustering (used for instance in [START_REF] Jens | Lets Do the Time Warp Again : Non-linear time series matching as a tool for sequentially structured data in ecology[END_REF]) can be applied directly on any MTS dissimilarity matrix, and does not require additional extensions.

-Model-based clustering of MTS has been addressed by [Kir05 ; GGM14] with multivariate HMMs

Extending an univariate time series clustering method is not straightforward, for several reasons. First, when several variables are present, a new information emerges : the relationship between variables. This information can take several form (same row-partition, same distribution, . . . ).

Another difficulty is that, when the number of variable is important, the model dimension also becomes high, which can cause inference issues and amplify the curse of dimensionality effects.

By assuming the presence of a latent variable partition, such that variables in a given cluster share the same distribution, the Latent Block Model addresses both the model dimension problem and the variable relationship modeling.

Multivariate Clustering With Block Structure

As described in Sect. The block clustering methods extend this principle by inferring a column-partition and one or several row-partitions. The crossing of these row and column partitions divides the MTS dataset into sub-matrices, called blocks.

In the following we detail the existing work on model-based coclustering and multiclustering.

Coclustering

The coclustering, also called biclustering, simultaneous clustering, two-way clustering or two-mode clustering, consists in estimating simultaneously a column partition and a row partition, which creates a block representation of a dataset. The associated structure is represented on Fig. 1.13-b).

The following section presents the related works on the topic, focusing on the models used in this thesis (see [BL15 ; PGA15], for exhaustive reviews).

Originally developed by [START_REF] Good | Categorization of classification. Mathematics and Computer Science in Biology and Medicine[END_REF] and [START_REF] John | Direct clustering of a data matrix[END_REF]), coclustering methods are specifically designed to deal with multivariate datasets. [START_REF] Govaert | Classification croisée[END_REF] and CROINFO [START_REF] Gérard | Co-clustering : models, algorithms and applications[END_REF]). For instance the CROEUC algorithm (designed for quantitative dataset coclustering) minimizes the least square criteria :

k,l x∈X k,l (x -µ k,l ) 2 ,
where µ k,l is the mean of block (k, l). Some other dissimilarity-based methods express the coclustering as a matrix factorization [START_REF] Lee | Algorithms for non-negative matrix factorization Advances in Neural Information Processing[END_REF], [START_REF] Chris | Orthogonal nonnegative matrix t-factorizations for clustering[END_REF] or an optimal transport problem [START_REF] Charlotte | Co-clustering through optimal transport[END_REF]. The problem has also been addressed with NN-based optimization [START_REF] Dongkuan | Deep co-clustering[END_REF] and evolutionary algorithm [START_REF] Stefan | An EA framework for biclustering of gene expression data[END_REF].

Existing coclustering variants includes the Stochastic Block Model [NS01 ; CL15] (which performs a block clustering based on a similarity graph observation-observation) and the diagonal Latent Block Model [Roo95 ; Li05 ; LN11 ; LN17] which adds the constraints of a block-diagonal structure.

In the following we detail the model-based coclustering approach.

Model-Based Coclustering

Model Definition

The Latent Block Model (LBM), first introduced in [GN03], assumes the presence of block components, such that the contents of a given block follow independently the same component distribution. The number of row-clusters K and column-clusters L are specified by the user.

The method has been applied to binary [START_REF] Gérard | Block clustering with Bernoulli mixture models : Comparison of different approaches[END_REF], contingency [START_REF] Gérard | Clustering of contingency table and mixture model[END_REF], ordinal [START_REF] Julien | Model-based co-clustering for ordinal data[END_REF] and continuous data [START_REF] Mohamed | Model-based co-clustering for continuous data[END_REF].

In the following, in addition to z = (z i ) n that designates the observation partition (equivalently called row-partition), the vector w = (w j ) p designates the variable partition (equivalently called column-partition) such that w j = l means that the variable x .,j belongs to the column-cluster l. p l denotes the column-cluster l cardinality. In the multivariate continuous case and with respect to this notation and the notation detailed in Sect. 1.4.1, the model likelihood is given by :

p(X) = z∈Z w∈W p(X, z, w) = z∈Z w∈W p(z)p(w)p(X| z, w),
where Z and W respectively denote the sets of all possible row and column partitions. The row-membership distribution p(z) is defined as n i=1 p(z i ) = n i=1 π z i , with π = (π k ) K the row proportions, and p(w) = p j=1 p(w j ) = p j=1 ρ w j , with ρ = (ρ l ) L the column proportions. The conditional density of X given the block memberships, is

p(X| z, w) = k,l x∈X k,l p(x | θ k,l ), with x ∼ F (θ k,l ). The complete parameter set is denoted Θ LBM = (π, ρ, (θ k,l ) K×L ).
The corresponding graphical model is displayed in Fig. 1.14.

Inference

The LBM likelihood optimization is not straightforward, and requires specific algorithms [START_REF] Gérard | Co-clustering : models, algorithms and applications[END_REF]. As in the clustering case, the inference process is usually performed in an Expectation-Maximization fashion (c.f., Sect. 1.2.1). However, the EM algorithm would require to access the block membership posterior distribution p(z, w | X, Θ), which is not tractable [START_REF] Gérard | Co-clustering : models, algorithms and applications[END_REF]. Alternative approaches include the variational inference (VEM), Stochastic-Gibbs EM (SEM) , and the Classification-EM (CEM) algorithm.

In the following we detail the SEM, that has the advantage to reduce both the risk of getting stuck in a local optima and the dependence to initialization.

SE step

Considering fixed values of the parameter set Θ, the block membership posterior distribution p(z, w | X, Θ) is approximated by a Gibbs sampler that alternate the sampling of z given w, and w given z. For a given row membership z i update, the sampling is performed by drawing from the conditional distribution p(z i | w, x i,. ) :

p(z i = k | w, x i,. ) = p(z i ) p(x i,. | z i = k, w) p(x i,. ) = π k p(x i,. | z i = k, w) K k=1 π h p(x i,. | z i = h) , where p(x i,. | z i = k) = p j=1 F k,w j (x i,j
) is the distribution of the row x i,. in the blockclusters. The column membership sampling is performed symmetrically, according to the following conditional distribution

p(w j = l | z, x .,j ) = π k p(x .,j | z i = k, z) K k=1 π h p(x .,j | z i = h) , M step
Considering fixed values of the block membership (z, w), the parameter set Θ update is, multivariate Gaussian case, given by :

n k = n i=1 1 k (z i ), p l = p j=1 1 l (w j ), n k,l = x∈X k,l 1, π k = n k n , ρ l = p l p , µ k,l = 1 n k,l x∈X k,l x, Σ k,l = 1 n k,l x∈X k,l (x -µ k,l )(x -µ k,l ) T .
The algorithm iterates these two steps for a given number of iteration or until membership convergence. This inference process is described in Algorithm 4.

As in the BNP clustering case, the final block membership can be inferred with a consensus partition estimation [START_REF] Thomas | Finding the mean in a partition distribution[END_REF]. The BIC criterion has also been derivated [START_REF] Aurore | Sélection de modèle pour la classification croisée de données continues[END_REF] for continuous data :

BIC LBM = ln p(x | Θ) - K -1 2 ln(n) - L -1 2 ln(p) - KLC 2 ln(np),
with C is the number of parameter per block.

Functional Latent Block Model

As stated in Sect. 1.5.2, the LBM has been applied to various types of data. Recently, several applications to time series clustering have been developed. The first application [START_REF] Faicel | Model-based co-clustering of multivariate functional data[END_REF] extends the clustering method [START_REF] Faicel | Joint segmentation of multivariate time series with hidden process regression for human activity recognition[END_REF] to coclustering, and assumes that the contents of a block can be represented by a piecewise polynomial regression model. This method supposes that the time series can be adequately modeled by piecewise representation, and the resulting segmentation can be a valuable information for the expert in that case. However, this method has the drawback to work in the original time series space, which becomes computationally prohibitive in high-dimension cases.

The other existing functional LBM [SAJ18 ; Bou+18 ; Sch+19] share the same time series preprocessing step, based on a fPCA representation that allows to work in a lowerdimensional space. One important difference between these methods is the modeling of the fPCA coefficient : while, in [START_REF] Ben | Model-based co-clustering for functional data[END_REF], the time series are represented by the first fPCA coefficient, the extended methods in [Bou+18 ; Sch+19] apply a subspace projection representation that associates a specific combination of fPCA coefficients to each block representation.

In [START_REF] Charles | Co-Clustering of Multivariate Functional Data for the Analysis of Air Pollution in the South of France[END_REF], the authors apply the same principle to MTS (i.e., each cell x i,j is itself an MTS). Very recently, the method [START_REF] Alessandro | Co-clustering of time-dependent data via Shape Invariant Model[END_REF] has been developed, that models the block components in the original space with the Shape Invariant Model.

Specific model selection criterion have been developed for the Functional Latent Block Model, including an ICL-BIC extension [START_REF] Ben | Model-based co-clustering for functional data[END_REF] and the ICL criterion [Bou+18 ; Cas+21], based on the existing ICL developed for LBM [START_REF] Aurore | Sélection de modèle pour la classification croisée de données continues[END_REF].

Non-Parametric Latent Block Model

As the DPMM can be seen as an extension of the MM with additional BNP prior, the Non-Parametric Latent Block Model can be seen as an extension of the LBM.

For the univariate case, the BNP extension of the LBM has been developed and presented in [START_REF] Edward | Nonparametric bayesian biclustering[END_REF], and consists in specifying two separate prior on the proportions, and a prior on the block component distribution. In [START_REF] Edward | Nonparametric bayesian biclustering[END_REF], the authors use Pitman-Yor (PY) processes as proportion priors (equivalent to the Dirichlet prior with an additional discount parameter, and which tends to create clusters with uniform sizes). With Dirichlet priors (used for instance in [SB08 ; WF12]), the model is defined by :

x i,j | z i , w j , θ z i ,w j ∼ F (θ z i ,w j ), z i | π ∼ M ult(π), w j | ρ ∼ M ult(ρ), π | α ∼ SB(α), ρ | β ∼ SB(β), θ k,l ∼ G 0 ,
with SB the Stick-Breaking generation process. The associated graphical model is displayed in Fig. 1.15. As in the clustering case, the inference can be performed in different ways, including variational inference, split-merge, and symmetrical Gibbs-sampler. This last inference method can be seen as two separate DPMM Gibbs-sampler inference on the rows and on the columns, where the objective is to approximates the joint distribution of the column and row memberships. The update of a row-membership z i given the

     n k n -1 + α p(x i,. | w, X -i , z -i , G 0 ), existing row-cluster k, (1.12) α n -1 + α p(x i,. | w, G 0 ), new row-cluster, (1.13)
Where

p(x i,. | w, X -i , z -i , G 0 ) = p j=1 F (x i,j , θ k,w j )
is the product of the block component densities, conditional to the fixed column-memberships values. In eq. (1.13), the joint prior predictive distribution of x i,. is p(x i,. | w, G 0 ) = L (m) l=1 p(x i,l | G 0 )), with L (m) the current number of column-partition in w at iteration m, and x i,l = {x i,j : w j = l} the elements of row x i,. belonging to column-cluster l. With this notation, the joint prior predictive of x i,l is the density of a multivariate t-distribution similar to the one used in Eq. (1.11).

The extension of the NPLBM to the case of time series coclustering is the topic of the contribution detailed Chapt. 4, which also details an adapted symmetrical Collapsed Gibbs Sampler.

Parametric or BNP versions of the LBM share one major assumption : there is only one row-partition, and, by construction, this row-partition is the same for every columnclusters. The next collection of methods relaxes this assumption and consider that every column-cluster can be associated with a specific row-partition.

Multi-Clustering

In a multivariate dataset, every variable taken independently is susceptible to produce a specific and distinct row-partition. The Multi-Clustering consists in assuming a small number of row-partitions are shared by the variables, which defines a column-partition. The associated structure is represented on Fig. 1.13-c).

It is different from clustering subspaces, which consists in estimating partition associated to specific linear variable combinations, in the way that Multi-Clustering estimates a variable partition such that variables in a column-cluster share the same row-partition. In some way, the Multi-Clustering can be considered a hard subspace clustering method, i.e., with hard column-cluster membership.

Parametric Multi-Clustering

The model-based Multi-Clustering has been developed by [START_REF] Giuliano | Model-based methods to identify multiple cluster structures in a data set[END_REF], which assumes a multi-partition block structure adapted to continuous multivariate datasets, where each cell is an univariate continuous observation. This method regroups the variable with a forward/backward search on the basis of a model selection criterion (in the example, variants of AIC and BIC), and models each column-cluster content with a full-covariance GMM. In this setting, the variables inside a block are considered dependent. The method has been extended in [START_REF] Giuliano | Modelling the role of variables in model-based cluster analysis[END_REF], which relaxes the block independence assumption, and separates the variables in three group : classifying, redundant and non-classifying.

Recently, [START_REF] Matthieu | A tractable multi-partitions clustering[END_REF] proposed a the Multiple Partitions Model (MPM) that deals with heterogeneous data (including categorical), and which consider the variable independence inside the block (resulting, for the Gaussian case, in diagonal covariance matrices). This model can be described by :

x i,j | {v j = h, z h i = k} ∼ F (θ h k ), v j ∼ M ult(η), z h j ∼ M ult(π h ),
where π h are the row-cluster proportions corresponding to partition h, η are the columncluster proportions and v = (v j ) p are the column-cluster membership. Note that, in this model, the block component distribution is a p h × d-dimensional distribution and not ddimensional, as each variable is modeled independently. The associated graphical model is displayed in Fig. 1.16. The authors propose two methods to perform simultaneously that maximizes a penalized log-likelihood. However, the authors state that these methods may suffer from combinatorial issues when the number of column-cluster is high (more than 5), in which case they suggest to use a forward/backward search.

Bayesian Non-Parametric Multi-Clustering

In parallel to the parametric model-based methods, the BNP Multi-Clustering (BNPMC) has also been developed in two independent works, the Cross-Categorization model from [START_REF] Vikash | Cross-categorization : A method for discovering multiple overlapping clusterings[END_REF] and the Non-Parametric Multiple Clustering model from [START_REF] Yue | Variational inference for nonparametric multiple clustering[END_REF]. These two papers share the model definition, which consists in a hierarchical Dirichlet prior : first on the column partition, which allows to automatically infer the number of column- The model, with graphical model displayed in Fig. 1.17, can be described as :

x i,j | {v j = h, z h i = k, θ h k } ∼ F (θ h k ), θ h k ∼ G 0 , v j ∼ M ult(η), z h i ∼ M ult(π h ), η j (r) = r j j-1 j ′ =1 1 -r j ′ , r j i.i.d. ∼ Beta(1, γ), π h j (t h ) = t h j j-1 j ′ =1 1 -t h j ′ , t h j i.i.d. ∼ Beta(1, α).
The inference is performed with a variational inference approach in [START_REF] Vikash | Cross-categorization : A method for discovering multiple overlapping clusterings[END_REF], and a MCMC-based in [START_REF] Yue | Variational inference for nonparametric multiple clustering[END_REF].

The contributions detailed in Chapt. 3 and Chapt. 5 consists in two Multi-Clustering models : the first one addresses the problem of time series Multi-Clustering with a parametric conditional latent block model, the second contribution is a BNP model that combines Multi-Clustering and Coclustering, and automatically infers the model selection without assumption on the partition dimensions. This knowledge is expressed as assumptions on the information that the user seeks to use for the cluster creation, and that should be conserved and outlined in the clustering result (e.g., a specific feature, a shape, the presence of a shapelet, the dissimilarity in a lower dimensional space, . . . ). When several variables are present in the datasets, additional assumptions are needed to deal with the dimension increase and make use of the variable dependencies information.

Conclusion

In the following chapters, the proposed univariate time series clustering is based on a latent scenario hypothesis, and the multivariate approach is based on a block clustering structure.

Time Series Clustering with

Model-Based Dictionary Representation 2

Introduction

In this chapter, we present an univariate time series clustering method applied to time series produced by autonomous driving numerical simulations. The method is based on a driving pattern dictionary construction and consists in three steps : automatic segmentation of each time-series, regime dictionary construction, and clustering of produced categorical sequences. In this chapter's first section, we present the detailed simulation method and the time series structure. In the second part, we discuss the existing approaches and describe our contribution. In the third section, we present the results obtained on public datasets and on an industrial use case : the Autonomous Emergency Braking (AEB) system validation. Finally, we conclude on our method's capabilities and perspectives.

Regime-changing Time Series Clustering

AEB use case time series are the result of the chaining of distinct phases, also known as regimes. Provided the ability to detect those regimes, it is possible to use their estimated distribution (order, frequency, amplitude. . . ) to characterize the observations and discriminate them. During the last decades, several papers have been proposed to detect optimal regime change points. Those methods sum up to piece-wise polynomial regression models. The common strategy relies on optimizing an approximation error in different ways : sliding windows of increasing size as in [START_REF] Eamonn | Segmenting time series : A survey and novel approach[END_REF] and [START_REF] Erich | Online segmentation of time series based on polynomial least-squares approximations[END_REF], by dynamic programming as in [START_REF] Marc | Least-squares estimation of an unknown number of shifts in a time series[END_REF], Hidden Markov models in [START_REF] Ath | A hidden Markov model segmentation procedure for hydrological and environmental time series[END_REF] or by regression mixture models in [START_REF] Faicel | Time series modeling by a regression approach based on a latent process[END_REF]. We selected this last model for two reasons : on the one hand, the benefits of using a mixture model (confidence intervals, model selection strategy. . . ) and on the other hand, the particular performances of this model compared to hidden Markov model approaches and its computational efficiency compared to dynamic programming methods [START_REF] Faicel | Time series modeling by a regression approach based on a latent process[END_REF][START_REF] Faicel | A hidden process regression model for functional data description. application to curve discrimination[END_REF].

In mixRHLP from [START_REF] Samé | Model-based clustering and segmentation of time series with changes in regime[END_REF], the same author combines the piece-wise regression model in a finite mixture to construct a one-step model-based clustering method. The proposed approach aims at regrouping time series with common regimes cut-points. mixRHLP also assumes that the number of regimes is known. These assumptions are not suitable to our use case, in which time series belonging to the same class can exhibit varying regime duration and cut-points. Moreover, the number of segment is not known in advances, nor the polynomial regression order.

Our contribution is an attempt to adapt mixRHLP to our constraints. It consists in a three-steps workflow with the addition of an original strategy of segmentation model selection. In the first step, we apply Individual time-series segmentation with a polynomial regression mixture. In the second step we build a standard dictionary of regimes by clustering the extracted segments. Finally, The clustering of these sequences using Levenshtein distance in categorical sequence space produces the final result. Our method, called SDLHC for Segmentation, Dictionary construction, Levenshtein Hierarchical Clustering, has the following advantages :

-Clustering based on regime detection is intuitive and easily interpretable by experts.

-The method can be applied to a dataset of time-series with unequal lengths. Moreover, it is independent of the time-series synchronicity and the regime's moment of appearance synchronicity.

-The segmentation phase can be applied independently on each time-series, which makes the computation an embarrassingly parallel task. This step drastically reduces the data dimension.

-Our segmentation strategy optimizes automatically both the number of segments and polynomial regression on each segment, which allows to automatizes this step and makes the segmentation step a turn-key solution.

A three-step time-series clustering algorithm (SDLHC)

The method SDLHC is composed of three steps : segmentation, dictionary construction, and categorical sequence clustering. The first two steps are addressed with a mixture model approach.

Segmenting time-series with a mixture of polynomial regressions

In a first step, each time series is individually segmented with a piecewise polynomial regression model with hidden logistic process from [START_REF] Faicel | Time series modeling by a regression approach based on a latent process[END_REF]. This model is based on a polynomial regression model mixture, with time-dependent proportions following a hidden logistic process. Given a time-series s = (s t ) T and ϕ = (ϕ q (t) = t q ) Q a functional basis of size Q. A Polynomial Regression Model (PRM) of sequence s in the basis ϕ is defined by

x = Q q=1 β q ϕ q (t) + σ 2 ϵ, with (β q ) q∈1,...,Q ∈ R Q , σ ∈ R +
* and ϵ ∼ N (0, 1). These PRMs are the segmentation mixture model components.

Given a number of components K, the time-indexed vector (z t ) T designates the segment memberships. At any given time t, z t follows a Multinomial distribution with parameters π(t) = (π k (t)) K . The distribution of s t is given by

p(s t ) = K k=1 π k (t)f θ k (s t ),
and the sequence s log-likelihood by

l(s; θ) = T t=1 log K k=1 π k (t)f θ k (s t ) , ( 2.1)
with f θ k (s t ) the density associated to a PRM component. With this notation, the timevarying proportions π k (t) can be seen as the parameters of a Multinomial distribution followed by the clusters memberships at a given time t. These proportions vary according to a logistic process. More formally, for k ∈ {1, . . . , K} and t ∈ T ,

π k (t) = p(z t = k) = exp( S s=1 w k,s ϕ q (t)) K h=1 exp( S s=1 w h,s ϕ q (t)) , ( 2.2)
with w k = (w k,q ) Q the associated model parameters. In the following paragraphs, we denote by W the set of parameters (w k ) K . The complete set of parameters is finally θ = (W, β, σ). The log-likelihood (2.1) optimization requires a specific version of the Expectation Maximization (EM) algorithm, with steps E and M described as follows.

Expectation step (E) Given the parameters θ, the first step of the EM algorithm consists in optimizing the complete log-likelihood defined as :

E x,θ [l(x, z; θ)] = E x,θ n i=1 K k=1 I z i =k log (p(x i , z i = k; θ)) = n i=1 K k=1 τ i,k log (π k f θ k (x i )) .
The development of the equation (2.3) shows that this step is simplified to the estimation of τ i,k = p(z i = k|x i ; θ), the posterior distribution of (z t ) T conditionally to s. The Bayes theorem gives the following estimation of this quantity :

τ t,k = p(z t = k|x i ; θ) = p(z t = k, s t ; θ) p(s t ) = π k f θ k (s t ) K h=1 π h f θ h (s t )
.

Maximization step (M) At each iteration, the model parameters are updated during the Maximization step. In this phase, the following decomposition of the complete loglikelihood expectation is maximized :

E x,θ [l(x, z; θ)] = T t=1 K k=1 τ t,k log (π k f θ k ,t (s t )) = T t=1 K k=1 τ t,k logπ k + T t=1 K k=1 τ t,k logf θ k ,t (s t ) = Q 1 (π) + Q 2 ((θ k ) k∈{1,...,
β k T t=1 τ t,k (s t - R r=1 β k ϕ r (t)) 2 , ( 2.3)
σ2 k = 1 T t=1 τ t,k T t=1 τ t,k (s t -μk (t)) 2 , ( 2.4)
with μk (t) = S s=1 βk,s ϕ q (t) the expected value of s t in the regression model of component k.

The maximization of Q 1 , which consists in estimating the time-dependent segment membership probabilities (π k (t)) K , is performed with an Iterative Reweighted Least Squares (IRLS) algorithm [START_REF] Faicel | A regression model with a hidden logistic process for feature extraction from time series[END_REF].

The algorithm has a linear complexity in the time series length and in the number of EM iteration, but quadratic in the polynomial regression order and cubic in the segment number. If this number is important (>15), we suggest to use a top-down hierarchical segmentation (not presented here) to reduce the effects of this cubic complexity.

Adaptive model selection strategy

In the initial model [START_REF] Samé | Model-based clustering and segmentation of time series with changes in regime[END_REF], the regression polynomial basis is common to every component, while in our contribution each regression order is specific. Moreover, we do not make a priori assumptions on the segment's number, which is also estimated by our strategy.

To estimate both the segment's number and the polynomial regression order on each segment, we combine this model with a new top-down strategy. This strategy is iterative and consists, at each step, in identifying the 'worst' component, in terms of the partial likelihood defined as :

l k (x; θ) = 1 T t=1 π t,k T t=1 π t,k log (f θ k (s t )) , k ∈ 1, .., K.
This criterion can be seen as the component representation quality weighted by the conditional membership probabilities. By improving the component k old ∈ {1, . . . , K} that minimizes this score, two candidate models are created and compared. Splitting k old in two sub-components, while conserving the other components, produce the first candidate model. We denote by k 1 and k 2 these new clusters. We denote t m the weighted median of the sequence {1, . . . , T } with weights π k old , and consider this time as the optimal cut-point for splitting the component π k old . The new components membership probabilities associated are based on the former component membership probabilities.

The membership probabilities of component k 1 are defined by :

π k 1 = π t,k old , t ∈ {1, . . . , t m } ϵ , t ∈ {t m + 1, . . . , T } , ( 2.5)
with ϵ the threshold precision. The new cluster k 2 membership probabilities are obtained likewise, with inverted time indices. A regularization of the (π k ) K is necessary at this point to enforce the constraint k=1 π k,t = 1, ∀t ∈ {1, . . . , T }.

The second candidate is obtained by increasing the polynomial regression associated to k old by one. Two runs of EM are then launched, each of them considering one of the candidates as the initial state. After the convergence of both EM, the candidate optimizing the BIC is selected for the next iteration. This strategy is summarized in algorithm 5.

Algorithm 5 : Top down segmentation strategy. Fix the convergence threshold c > 0 Choose an initial state for the first EM run : After convergence of BIC criterion, the regime change cutpoint are estimated by postprocessing the time-dependent memberships probabilities. This segmentation method is applied independently to each time-series and transforms each one in a set of subsequences.

θ init old := ((w k , β k , σ 2 k ) k∈{1,...,

Dictionary construction

Expressing the extracted segment in a common basis is mandatory to compare and cluster the sequences. The objective is to encode the original time-series in the new dictionary, as represented in the left part of Figure 2.1. This common basis, or dictionary, is constructed with a clustering algorithm applied to the dataset composed of all segments from time series.

The sub-segments are first scaled, expressed on a common support, and regressed in a polynomial regression basis. Other informative descriptors can be added depending on the case, as the regime's duration, offset, or variance. These features are then clustered with a GMM to produce the dictionary. In section 2.2.1, we mentioned an implicit assumption based on the segmentation polynomial basis. In this section we make the additional implicit assumption that the GMM is adapted to the regimes density estimation and makes sense from the field expert point of view.

The EM algorithm is initialized with the k-means++ algorithm, which is a standard approach [START_REF] Johannes | Simple methods for initializing the em algorithm for gaussian mixture models[END_REF]. At the end of this step, the modes of the Gaussian mixture components are the reference regimes, entitled "patterns" in the following, with which to recode the original time series. The dictionary size is determined by the field experts, assisted by the BIC. The final dictionary is denoted R = (r u ) U , with U the number of patterns.

After this re-coding phase, data dimension is greatly reduced : for a time-series of size n, the dimension goes from R n to R K , with R the categorical space and K the number of regimes composing the sequences. The third and last step of SDLHC regroups these sequences to produce the final clustering result.

Categorical Sequences Clustering

We use the Levenshtein distance [START_REF] Vladimir | Binary codes capable of correcting deletions, insertions, and reversals[END_REF] combined with Ward's hierarchical clustering method to obtain the final clusters. Levenshtein distance between two categorical sequences a and b (respectively composed ofK a and K b segments) is defined as the minimum number of operations (insertion, deletion, substitution) needed to transform a into b. In this categorical space, Levenshtein Distance complexity is O(K a × K b ).

However, with the original Levenshtein distance, replacing a symbol with another has a fixed unit cost, independent from the target and replacement symbols. Therefore, it does not take into account the possible similarity between patterns.

For instance, in Fig. 2.4, some speed patterns appear similar (different phases of acceleration of cruising speed) and distance on categorical sequences should take this similarity into account.

Weighted Levenshtein Distance

In order to take this information into account, we propose to use a Weighted Levenshtein Distance (WLD). Considering a set of patterns R, the weight on the edition cost between r 1 , r 2 ∈ R is symmetric and defined by :

||r 1 -r 2 || p max ra,r b ∈U ||r a -r b || p , ( 2.6)
where ||.|| p is the p-norm on the pattern space. The choice of the p-norm influences moderately the final clustering. Based on our practical experience in the AEB use case, we advocate the use of ||.|| ∞ .

Variant representations

The following paragraphs describe other possible representations and dissimilarities, that can be use in combination or in replacement of WLD.

Removing duplicate successive regimes

The categorical sequence contains sometimes successive identical symbol (e.g., the sequence abaac). This effect can have several origins : it can be that the time series exhibits two close (but different) regimes, that are clustered together in the dictionary construction step. It can also be caused by an overestimation of the segment number in the independent segmentation step.

To mitigate these effects, we suggest to assume that the successive symbols do not add any information and to delete the redundant symbols (e.g., to transform abaac into abac).

Integrating additional features in the Levenshtein Distance

During the second phase of SDLHC, the dictionary has been constructed based on scaled and segments on same support, with the optional addition of the offset, variance and phase duration information. These characteristics can also be integrated in the weighting function of the Levenshtein Distance.

It is also possible to assume a hierarchy of these features, and consider a nested clustering (for instance, a partition based on pattern shape, then duration, then offset, . . . ).

Histogram-based representation

A common representation choice in dictionary approaches consists in considering the recoded symbol histograms [START_REF] Jessica | Finding structural similarity in time series data using bagof-patterns representation[END_REF][START_REF] Patrick | The BOSS is concerned with time series classification in the presence of noise[END_REF]. This representation can also be used with our approach. However, as the number of patterns may be sensitive to the over-estimation of the segment number or of dictionary size, we suggest to also take into account the proximity between patterns in the comparison of these histograms. The final clustering can be obtained with the same hierarchical clustering approach based on the ℓ 2 norm.

Based on the same logic, the frequency or an absence/presence binary value representations can be considered instead of the histograms.

Integrating temporal relevance information

In some situations, the user may be interested in comparing specific temporal sub-sequences of the time series. For instance, the last moments of the emergency braking scenarios are often more relevant. Instead of extracting the sub-sequences before SDLHC run, it is possible to integrate this prior information in this last symbol sequence clustering step. For instance, by weighting the ℓ 2 norm in the histogram dissimilarity or in adding a temporal weighting in the edit cost of the Levenshtein Distance computation.

Once the weighted Levenshtein Distance Matrix computed, Ward's hierarchical clustering method is applied to produce the final clusters.

In the following section we compare SDLHC (with the WLD) to other state-of-the-art methods.

Experiments

We present, in this section, the results of several experiments on public datasets and on a real-world use case AEB obtained from Renault's simulation system. The method described in this article was implemented in Scala for the segmentation step and R for the hierarchical step. Code and (public) datasets available at https ://tinyurl.com/sdlhc. The following baseline methods are selected :

-Three methods based on classical measures (Euclidean distance and DTW) associated with Partitional Around Medoid clustering approach. We have also tested the combination of DTW with center construction using the popular Dynamic Barycenter Averaging (DBA) method [START_REF] François | A global averaging method for dynamic time warping, with applications to clustering[END_REF]. These methods are denoted respectively ℓ 2 , DTW p , and DTW d in the following.

-The K -Shape method [PG15], a partitional clustering using the shape-based distance based on the cross-correlation measure.

-The SAX method, a dictionary-based methods from [START_REF] Jessica | A symbolic representation of time series, with implications for streaming algorithms[END_REF] that builds representations of the time series based on uniform time step segmentation. Based on this representation and associated distance, hierarchical clustering with Ward's criterion produce the clusters.

-In order to compare to the original method we aimed to extend, the results of mixRHLP are also reproduced here.

Whenever needed, we interpolate time-series to equal-length sequences. We used the R package TSclust's distance-based and SAX methods implementations and mixRHLP using flamingos R package. Some of these methods depend on parameters, usually estimated by optimizing a risk in a supervised framework. The ARIs obtained here are always the maximal ARI obtained when testing the method on a parameter grid, displayed in Table 2.1, reproducing the results that experts can obtain after fine-tuning. The comparison is based on the Adjusted Rand Index (ARI), a popular score in the clustering validation context. This criterion represents the proportion of correctly grouped and separated observations with respect to the observed classes.

Public datasets results

In order to validate SDLHC adequation to the regime-changing time series clustering problematic, we selected a subset of the UCR archive [START_REF] Hoang Anh | The UCR Time Series Classification Archive[END_REF] whose data exhibit regime structure. The ARI score obtained are shown in Figure 2.2. Although performant when applied to Renault's dataset (c.f. next subsection), we found that the weighted Levenshtein hierarchical clustering requires fine-tuning to adapt to the considered data characteristics. The test ran in this section therefore use the non-weighted Levenshtein distance.

The results confirm that the method perform well when addressing regime-changing time series. In these tests, the considered datasets contain equal-length time-series. However, SDLHC can also be applied, without data preprocessing, to unequal-length timeseries, which is the case in our application.

Real dataset results

In the following section, we evaluate the clustering performance of SDLHC on an industrial use case : the Autonomous Emergency Braking (AEB) system validation. In this case, a ground truth is available, and it is possible to compare clustering methods based on the similarity between the observed labels and the produced clusters. The clustering methods performances are, as in the previous section, measured by the ARI score. Renault's dataset is composed of 150 time series, with a duration varying from 13 to 52 seconds and length varying from 415 to 573 data points. An illustration of the segmentation step result is shown in Figure 2.2.

The obtained dictionary, composed of five driving patterns, is shown in Figure 2.4. Using this dictionary, Figure 2.3 shows the encoded sequence.

Two stationary patterns can be recognized (b and c), corresponding to cruise speed phases, as well as two accelerating (d and e) and one decelerating (a).

The scores are obtained in the same conditions than the previous tests on public datasets, displayed in Figure 2.5.

Two versions of SDLHC are tested : SDLHC -LEV and SDLHC -W LEV corresponding to the use of the standard and weighted Levenshtein distance in SDLHC's last step. ARI criterion confirms that the SDLHC -W LEV method slightly improves the score obtained by SDLHC -LEV . Among the distance-based methods, the K -Shapes method is the best performer without, however, reaching the ARI threshold of 0.45 regardless of the number of clusters. With high cluster numbers, SAX method nearly reaches the performance of SDLHC -LEV . This result seems logical given the proximity between the proposed workflow and the dictionary-based methods.

Conclusion

This first contribution SDLHC performs an univariate time series clustering based on the independent segmentation of time series, a common segment pattern dictionary construction, and a categorical sequences clustering.

This method makes several assumptions : a) the time series are composed of several phases and these phases can be extracted with a piecewise polynomial regression ; b) A small number of pattern can adequately represent all the phases from every time series and these patterns can be obtained with a GMM on the segments representation coefficients ; c) a relevant partition of driving behaviors can be obtained with a weighted levenshtein hierarchical clustering of the time series recoded as categorical sequences.

It is our strong belief that SDLHC can be applied to other domains complying with the latent scenario hypothesis (e.g., human activity recognition).

One drawback of this contribution, however, is that it only addresses univariate time series dataset, while the ADAS validation use cases are mainly composed of multivariate data. Because the direct extension of SDLHC is not straightforward (c.f. perspective discussions in Sect. 6.1), we considered a different strategy for these use cases. This strategy, the model-based block-clustering framework, is the topic of the next three chapters.

Multivariate Time Series And Functional Conditional Latent Block Model 3

This chapter details the first block clustering method contribution of this thesis, which consists in a model-based parametric multi-clustering model.

Introduction

Coclustering methods infers simultaneously one row partition and one column partition (c.f. Sect. 1.5.1). However, this class of method makes one strong structural assumption : in addition of the variable partition, there is only one observation partition. This assumption means that the observation partition is shared among every column cluster. When the true latent block structure is composed of several row-partitions, inferring a coclustering forces the creation of a unique consensus row-partition made of the crossing of the multi-clustering row-partitions. This effect leads to the over-estimation of the consensus partition cluster number, and to the inappropriate split of block components. This effect is illustrated on Fig. 3.1, with three variable clusters, each containing containing 3 row-clusters.

When the multi-clustering structure contains an important number of variable clusters, this over-estimation can severely hinder the interpretation.

In the following we present a new model, called Functional Conditional Latent Block Model (FunCLBM), a multivariate time series method that extends the Functional Latent Block Model with a dependency structure that allows each column-cluster to be associated with a specific row-partition.

As opposed to existing model-based multi-clustering from [START_REF] Giuliano | Model-based methods to identify multiple cluster structures in a data set[END_REF] or [START_REF] Matthieu | A tractable multi-partitions clustering[END_REF], that handle univariate cells elements and where each variable follow a different distribution, in CLBM the cells belonging to the same block follow independently the same multivariate distribution, as is the case in the LBM [GN03 ; GN08] or FunLBM [START_REF] Ben | Model-based co-clustering for functional data[END_REF]. This difference is outlined in Fig. 3.2. This choice of model combines the efficiency of the LBM parameterization sparsity with the multi-clustering framework.

The output of this method is a dataset structure that highlights the variables that are linked together and the multiple partitions, from which it is easy to discard groups of uninformative (i.e., the one that are associated with a one-component row-partition) or uninteresting groups of variables.

As a model-based method, this approach allows to detect outlier observation from several point of view, and to precisely understand which group of variable is associated to the outlier values.

Functional Conditional Latent Block Models

This section presents the FunCLBM model, as well as its inference and model selection strategies. The proposed approach relies on the projection of the time series in a specific space.

Representation with Principal Components of Interpolated Periodogram

As described in the Related Work section, there is a important number of existing time series representations (c.f., Sect. 1.1.4). In this chapter we chose to use an interpolated log-scaled Fourier periodogram representation, as it proves its ability to compactly represent the time series obtained in our driving simulation context.

Definition Each time series S = (s i,j ) n×p is first independently represented as a periodogram with a Fourier Transform. Given a time series s i,j with length T i , this transformation outputs a decomposition of the time series in a weighted combination of trigonometric functions, with f i,j the associated frequency basis. Each frequency basis has the same dimension than the original time series, T i , and the associated frequencies are different from one time series to the other (i.e., the time series are represented in different basis).

In order to compare the time series, it is mandatory to construct a common representation basis f , based on (f i,j ) n×p .

A possible solution consists in computing the sample average gap ∆f between two consecutive frequencies over all time series, and to build an "average" frequency basis f = ⟨0, ∆f , 2 ∆f , . . . , ( l -1) ∆f ⟩ that best represents the set of time series.

The periodogram values in this common basis can be estimated using linear or cubic interpolation techniques [START_REF] Jorge | Comparison of times series with unequal length in the frequency domain[END_REF] from B. After this step, as advocated in [START_REF] Jorge | Comparison of times series with unequal length in the frequency domain[END_REF], we take the logarithm transformation of the normalized periodograms, which eases the detection of the dependency structure and reduces the importance of the scale.

As a final transformation, the obtained log-normalized-interpolated periodograms are linearly projected with a PCA, which produces the dataset X, such that each cell x i,j contains the representation of time series S i,j . Every cell has the same length d, corresponding to the number of PCA axis kept. The number of PCA axis, and the periodogram dimension are both hyper-parameter oh these representations, that are studied in the next chapter.

This preprocessing method is a close variant of the fPCA used in [SAJ18 ; Bou+18], with a trigonometric basis of function and an additional log normalization transformation.

Interpretation As stated in this thesis introduction, it is important for the user to understand the underlying transformations that compose the clustering algorithms. Therefore, it is important to detail the characteristics of the produced representations, and the link between the similarity in the final representation space and the original time series similarity.

Because the time series are represented in a frequency basis, this transformation is suited to capture any seasonality information, but also the regime change similarity. Because the periodogram are normalized, the transformation reduces the importance of the time series scale.

These elements can be observed in practice in the use cases : some variables grouped together may have different shapes, but they always exhibit an important correlation and contains synchronous evolution and events. For instance, a variable that records a braking system activation (which is a temporal binary variable exhibiting a rectangular shape) is often seen grouped with the braking strength.

Model definition

Let denote K h the number of row-clusters associated to column-cluster h, 1 ≤ h ≤ H. We also denote by v = (v j ) p the column membership vector and η the column proportions. Given a column-cluster h, the associated row-clusters partition is denoted as z h = (z h i ) n×K h . We denote z h i the row-cluster membership of observation i in the column cluster h. The complete set of row-partitions is denoted as Z = (z h ) H . The row proportions are denoted π = (π h ) H , with π h = (π h k ) 0≤k≤K h . Finally, the joint model density can be described by :

p(X) = V p(v)p(x | v) = V j η v j H p(X h ),
where V is the set of all possible column-partition in H clusters, and p(X h ) is the density of the sub-matrix containing only the variable belonging to column-cluster h. This density can, in turn, be expressed as :

p(X h ) = Z h p(z h )p(X h | z h ) = Z h i η z h j k x∈X h k p(x | v j = h, z h i = k),
where Z h is the set of all possible row-partition in K h clusters, and p(x

| v j = h, z h i = k) = f (x, θ h k )
is the block-component density, in our case a d-dimensional multivariate Gaussian density with parameter θ h k = (µ h k , Σ h k ). The complete set of parameter θ=(η, π, (θ h k ) K h ×H ) is inferred with a dedicated SEM-Gibbs algorithm.

Inference with SEM-Gibbs algorithm

As detailed in the related work chapter (c.f. Sect 1.2), the SEM algorithm is popular practice in the model-based clustering framework. As is the case for the coclustering, the direct optimization of the likelihood is not straightforward, and requires specific strategies.

In the following we propose an SEM-Gibbs algorithm that alternates the Gibbs sampling of the proportions and the parameter inference. Starting from an initial parameter state θ 0 and an initial column partition w 0 , the algorithm alternates between these two steps :

1. SE step :

-Update the column-partition v given the row-partitions. For each column j, the update of the column membership v j consists in sampling from p(

v j = h | x .,j , θ, Z) ∝ η h p x .,j | v j = h, z h ; θ H r=1 η r p (x .,j | v j = r, z r ; θ)
, where p x .,j |z h ; θ is the density of column j in cluster-column h given the row-cluster memberships, i.e., p x .,j

|z h ; θ = n i=1 f (x i,j , θ h z h i ).
-Update the row-partitions Z given the column-partition. For each columncluster h and each row i, the row-membership

z h i update consists in sampling a new value from p(z h i = k | x i,. , v, θ) ∝ p x i,h | z h i = k, θ) = π h k p x i,h | z h i = k, θ H s=1 π s i p x i,h | z h i = k, θ
, where x i,h designates the content of row x i,. restricted to column-cluster h, and

p x i,h | z h i = k, θ = j:v j =h f (x i,j , θ h k )
is the density of the row i in block-cluster (k, h).

M

Step : given the sampled block partition, and denoting by X h k the observations belonging to block (k, h), the mixture proportions are updated by :

π h k = 1 n i 1 k (z h i ), η h = 1 p j 1 h (v j ), µ h k = 1 n h k i,j z h ik w jl v ij , Σ h k = 1 n h k i,j z h ik w jl v ij -µ h k v ij -µ h k T
The SEM algorithm alternates these two steps for a given number of iterations, and output a set of sample (v, Ẑ).

As in the mixture model and coclustering case (c.f. Sect. 1.2.1, the initialization choice is crucial to ensure the good behavior of the algorithm.

Several methods are often considered : populating components with a small random sample of the observations, shuffling the column and block partitions, or using another clustering algorithm to get a good initial starting point. In sect. 3.3, these different initializations are experimented.

Model Selection

With a good initialization choice, SEM-Gibbs may converge to a solution for a given clustering structure, i.e. a column cluster number H and a set of row clusters numbers K = (K h ) H . Several criteria have been developed to address the model selection problem. In this work, we propose a dedicated criterion based on the Integrated Classification Likelihood (ICL) [START_REF] Christophe | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. Initially developed for GMM Selection, extended by [START_REF] Aurore | Sélection de modèle pour la classification croisée de données continues[END_REF] to coclustering and in [START_REF] Charles | The functional latent block model for the co-clustering of electricity consumption curves[END_REF] to functional coclustering, we propose the following extension to functional Multi-clustering :

ICL(K, H) = log p(x, v, ŵ; θ) - H -1 2 log p - 1 2 h ((K h -1) log n) - h,k ν h k 2 log(np),
where ν h k is the component parameter number of block (k, h). This score penalizes the log-likelihood with a function of the number of parameters. The best model is the one maximizing this score.

In the coclustering case, finding the best structure can be done by an exhaustive grid search. This strategy cannot be applied in the Multi-Clustering case, where the number of possible model is huge, even for a small number of blocks.

With a fixed number of column-clusters H, the number of possible model with a maximum of K m row-clusters is the number of possible combinations of H elements taken from a collection of K m objects, with repetition (because several column cluster can share the same number of row-cluster -e.g., (2 × 2 × 2) is a valid model choice-), and without order (because the model likelihood is invariant to label switching -e.g., choosing the model dimensions (1 × 2 × 3) is equivalent to choosing (3 × 2 × 1)-). The exhaustive combination list for H = 2 and K m = 3 is given in Table 3 The number of possible combinations is the number of multisets Km h = Km+h-1 h . For a maximum number of column cluster H m and a maximum number of row-cluster K m , the total number of model is therefore given by :

Hm h=1 K m h = Hm h=1 K m + h -1 h .
This number grows quite fast, and becomes huge for low values of H m and K m (e.g., for H m = K m = 10, it amounts to 184755 combinations, and, for H m = K m = 20, 1.37E12 combinations). Estimations of this number are given in Table 3.2 for a range of low values.

As a consequence of this important number of models, it is not possible to test the combinations exhaustively, and alternative model selection strategies must be considered. Our proposal of heuristics consist in estimating the column-partition with a coclustering model selection, and then to estimate independently the best model for each columncluster.

In the next section, two of these approaches are tested on a simulated dataset, based on an LBM greedy search.

Experiments on Synthetic Data

In order to test the capabilities of FunCLBM, several experiments are first conducted on a simulated dataset. These experiments help us verify that the model is suited to the use case and that the SEM-Gibbs algorithm behaves as expected in a controlled environment.

The Scala source code is available at https ://tinyurl.com/FunCLBM, along with the data simulation script.

Simulated dataset

The first experiment is conducted on a dataset sampled from a known generative model. The objective is to check the behavior of FunCLBM, its initialization and model selection strategies. The dataset is generated by sampling around one of several "prototypes" denoted (ϕ h k ) and which represents the components modes in the original space. For each block (k, h), several time series are drawn following N ϕ kh (t + t s ), s 2 with s = 0.02 and t s a random shift ∼ N 0, s 2 . These modes are depicted in Fig. 3.3 according to the dataset structure.

In the experiments, the quality of the estimated block partition is compared to the known generative partition, based on the Adjusted Rand Index (ARI). This is a popular criterion choice in the clustering domain, which represents the proportion of correctly grouped and separated observations with respect to the observed classes. In our particular context, we compare the obtained partition based on three aspects : the column cluster partition, the rows cluster partitions (made of the binning of every row cluster partition per column), and the block partition. We generate a dataset of size 90x90, with column cluster of size (45, 15, 30) and row cluster sizes of respective sizes (20, 40, 30), (60, 30) and (40, 50).

Model Adequacy

As a preliminary test, we first verify the ability of FunCLBM to regroup time series after the pre-processing step and with the latent block structure hypothesis. To do so, we compare the model Log-likelihood produced after 100 launches of SEM-Gibbs inference to the corresponding ARI scores. The results, shown in Fig. 3.4, seem to indicate a correlation between ARI and likelihood, which tends to confirm the model suitability.

We verify this relationship with the Pearson's correlation coefficient between the two scores, and Kendall's correlation test. We use this latter test on scores ranks to avoid making assumptions on ARI or Log-likelihood normality and because of the presence of ex-aequos values (that can be produced if the "true state" is reached).

The results (with 95 % confidence level for the Kendall test) are displayed in Table 3 correlation for every partition dimension and by Kendall's correlation test p-value at 95% confidence level.

Initialization

The next experiments aim at evaluating the different initialization strategies. The clustering quality is again evaluated with the ARI criterion after 30 runs of SEM-Gibbs. We compare four strategies :

-Populate blocks with samples : the model is initialized by sampling a small number of variables for each column-cluster, and then a small number of rows for rowclusters.

-Random shuffle of the column partition, then of the row partition.

-One k-means algorithm run on the variables (i.e., on the transposed dataset), then one k-means to estimate each row-partition.

-Initialize w with the column-partition obtained after an LBM run, then the rowpartitions at random.

The results distributions are described in Fig. 3.5.

The figure shows important differences between row and column cluster results, as expected since the model does not treat rows and column symmetrically. The k-means algorithm has an unexpected behavior : while performing well on average, the results show a high dispersion for column cluster ARI. Its row cluster ARI however is slightly better than the other methods.

The FunLBM initialization strategy gives the best results on average, closely followed by the randomPartition strategy. For real-life applications, we suggest to use FunLBM when the datasets dimensions allows it, and otherwise the randomPartition strategy, which is cheaper alternative. Whether with the FunLBM or randomPartition initialization, an additional strategy can be to select the best result (w.r.t. model likelihood) among several runs. The following experiment highlights the interest of this approach, by comparing the ARI score increase with respect to the number of candidate runs. The results are displayed in Fig. 3.6.

Row clusters

Column clusters

The results shows that, whichever initialization strategy is applied, the concurrent run of several methods allows to stabilize the results, and finding the perfect structure is an easy task whenever the number of concurrent launches is higher than 4.

Model selection

The previous experiments compared the initialization methods for a given choice of model. In the following, we experiment two approaches of model selection, which is a challenging aspect of this method, and an important feature for the field expert. In [START_REF] Charles | The functional latent block model for the co-clustering of electricity consumption curves[END_REF], the coclustering model selection is performed with a grid search strategy, which requires inferring K × H models. As stated in Sect. 3.2.4, this exhaustive search is not feasible in the multi-clustering case, because of the prohibitively high number of model possibility.

To overcome this limitation, we propose and compare two strategies that both avoids a direct grid-search of every multi-clustering combinations. These two strategies are based on a different estimation of the column-cluster number H and then on a column-wise grid search.

In the first case, ŵ and H are estimated from an LBM exhaustive grid search, which implies an exhaustive search of K M × H M LBM models and then H × K M to estimate the best number of row-clusters per column-cluster. The total number of model to infer is bounded by 2K m H m (K m H m for the LBM grid search, then K m possible models for each of the selected column number bounded by H m ).

The second strategy is an iterative algorithm that chooses, at each iteration, the best functional Latent Block Model between the one with an added row cluster and the one with an added column cluster, based on the ICL criterion [START_REF] Aurore | Sélection de modèle pour la classification croisée de données continues[END_REF][START_REF] Charles | The functional latent block model for the co-clustering of electricity consumption curves[END_REF]. The search stops when the ICL stops increasing. In this LBM search, the number of model to infer is bounded by K m + H m , and the total model number is

K m H m + K m + H m .
Each methods are run 50 times, and the model selection is evaluated with the model dimension difference to the ground-truth model (in this case, a (3 × 2 × 2) model). These differences, displayed in Fig. 3.7, show that the LBM grid search performs well on that case, and is a valid strategy for low model dimensions candidates.

Application to autonomous driving system validation 3.4.1 Use case description

In this situation, the objective is to test the reactions of a car (called Ego) equipped with the control logic. Ego runs in a straight line and starts drifting laterally towards the road side or the other lane, simulating a sleeping driver. We expect the drifting detection system to trigger the control logic, which in turn puts the car back in its line center, as an emergency maneuver. The situation is depicted in Fig. 3.8. The simulated datasets contain the data from 56 simulations, each described by 20 signals. Some signals are duplicated in order to test FunCLBM ability to regroup them, and some uninformative ones are kept on purpose.

Results

The experiments on simulated datasets lead us to choose the following setup for the real case analysis : the initialization is always performed by sampling column and row cluster partitions, and the FunLBM grid search approach is applied for model selection.

Each combination is tested with 30 concurrent runs. In the second, and largest, column-cluster can be found the uninformative signals, that either give constant values (vehicle length, width, distance between wheels, road bend radius) or increasing linearly (distance to origin). Fig. 3.14 and Fig. 3.15 illustrates some of these signals. The third column-cluster regroups two other interesting features : the rectangular function indicating the activation of the control logic, and the changes in Ego's heading. While the first column-cluster was grouping position, this one gathers the control features. Fig. 3.13 shows the content of subcluster (3,3) which illustrates the relationship between them. Overall, every set of duplicate features have been correctly grouped together.

This conditional clustering partition shows, in conclusion, that the FunCLBM approach has correctly discriminated uninformative signals, while creating meaningful clusters of features (position and leverage). In each column-clusters, the observation clusters are relevantly regrouped and provide insights of the dataset content.

Conclusion

This chapter describes FunCLBM, a model-based method that addresses the problem of clustering multivariate time series in multi-views. This new model enables regrouping redundant signals, discriminating uninformative ones and provides the user with multiple clustering views of a multivariate time series dataset. The time series are transformed into interpolated log-periodogram before being projected into low-dimensional space.

Several initialization methods and model selection strategies are proposed and experimented on a simulated dataset. The experiments show the model adequacy and give insights on the most interesting implementation strategies. Finally, we apply the method to a real-case dataset from the autonomous driving system validation domain. In this application, FunCLBM has been able to simultaneously discriminate groups of signals and produce meaningful driving behavior clusters. These results shows the usefulness of the model and the effectiveness of the initialization and model selection strategy.

A major drawback of this method is the model selection aspect, that cannot be performed without the help of specific strategies. The proposed greedy search algorithms allow to explore a part of the model space but at the cost of additional structural assumptions.

For instance, the LBM grid search-based model selection assumes that the LBM with the ICL can adequately estimate the true number of cluster column (which, by assumption, have different row-partitions).

In addition, for higher numbers of clusters (e.g., H m , K m >= 20), these strategies still involves the inference of hundreds of models.

One solution to circumvent this model selection limitation is to use the BNP framework. By adding a BNP prior to the existing block structure, it is possible to estimate simultaneously the model dimensions and parameters.

Bayesian Nonparametric Coclustering 4

This chapter describes a Bayesian Non Parametric Functional Latent Block Model, which is probabilistic coclustering that extends the recent Functional Latent Block Model [START_REF] Ben | Model-based co-clustering for functional data[END_REF] in a Bayesian Non Parametric framework. This model also serves as basis component of the next chapter contribution.

Introduction

As stated in Sect. 1.5.2, LBM methods for time series coclustering, where each cell is a temporally-indexed vector, have only been introduced recently. The existing modelbased methods [CB17 ; SAJ18 ; Bou+18 ; Sch+19 ; Bou+20 ; Cas+21] are parametric model, and assume that the number of blocks is known a priori. This assumption is rarely true in practice, and these dimensions must be estimated with an additional model selection step. This selection is usually performed either with a grid-search or with a greedy strategy that hierarchically explore the coclustering solutions. These strategies present several drawbacks :

-In the grid-search selection context, the computation cost can be prohibitive as every combination of block number is tested, and the user is never certain that the true model is within the grid -The greedy optimization heuristic is sub-optimal, by picking iteratively local optima and assuming a hierarchical structure of the mixture components -Whether with the greedy optimization or grid-search, choosing the model selection criterion [START_REF] Gilles | Model selection for mixture modelsperspectives and strategies[END_REF]For+21] is not an easy task and influences the final results

The Dirichlet Process Mixture Model (DPMM), described in 1.3.2, is a Bayesian BNP model-based clustering approach that can infer the number of latent clusters. As a nonparametric model, its parameter set dimension may increase indefinitely with the dataset size. This property makes it extremely interesting for massive dataset exploration, especially when the user can augment the dataset content at will. This is precisely our situation in the driving validation context, in which it is possible to allocate additional resources to explore specific areas of the observation space.

Non-parametric approaches to LBM (NPLBM) have been studied in few works [MR07 ; WF12], but, to the best of our knowledge, never applied to multivariate time series coclustering use cases. This paper proposes the functional non-parametric Latent Block Model (FunNPLBM) method, the first non-parametric model-based method applied to time-series coclustering, which closes the gap between FunLBM and NPLBM. In addition, our contribution includes a practical use case illustrating the method's interest, a more compact definition of the NPLBM, a detailed study of the method's capacities based on benchmark and experiments, and Scala source code put at disposal for reproducibility.

In the following, Sect. In the following, we recall that the dataset X is an interpolated periodogram fPCA representation of the vectorized multivariate time series dataset as in the previous chapter (c.f. Sect. 3.2.1) .

Model Definition

In another context than the multivariate time series analysis, [START_REF] Edward | Nonparametric bayesian biclustering[END_REF] proposed a definition of the NPLBM. This work assumes Pitman-Yor Process (PYP) priors for the rowmemberships and column-membership. However, PYP distributions (as DP) are distributions over parameters and not on memberships. Using PYP, the model from [START_REF] Edward | Nonparametric bayesian biclustering[END_REF] implicitly defines different sets of parameter distributions that are not linked to the block distributions.

In the following, we propose a model definition of the NPLBM. Instead of generating each cell independently, this definition models the overall dataset X. Given fixed rowpartition and column-partition z, w and given the block parameter matrix (with dimension 1 ≤ K ≤ n and 1 ≤ L ≤ p, Θ = (θ k,l ) K×L , the element x i,j follows the distribution associated with its block membership. In our application, the block distribution F (theta z j ,w j ) are multivariate Gaussian.

The overall process is defined by the following expressions :

x i,j | {z, w, Θ} ∼ F (θ z i ,w j ), z i ∼ M ult(π), w j ∼ M ult(ρ), π ∼ SB(α), ρ ∼ SB(β), θ k,l ∼ G 0 .
With this definition, the generation of the matrix X is done by generating the proportions π and ρ, sampling the row-memberships z ∼ M ult(π) and column-memberships w ∼ M ult(ρ) separately, then sampling the block parameter matrix Θ (in practice, only the ones associated with observed block defined by given z and w) and finally drawing the cells value x i,j from F (θ i,j ). The likelihood of X is given by p(X | z, w, Θ) = i,j p(x i,j | θ z i ,w j ), and the joint prior distribution of the hidden variables by :

p(z, w, π, ρ, Θ | G 0 , α, β) = p(z | π) p(π | α) p(w | ρ) pρ | β) p(Θ | G 0 ).
In the next section, we describe the two-steps collapsed Gibbs sampling used to simulate draws from the posterior.

Inference

The objective is to obtain good estimates of z, w, based on the prior (α, β, G 0 ) and the observed dataset X, by sampling from the joint posterior distribution p(z, w | X, G 0 , α, β). Sampling directly from this joint distribution is not feasible, but can be approximated with a Gibbs sampler that iteratively update the values of the memberships z given w, X, G 0 , α, and then w given z, X, G 0 , β. In our case we use a collapsed version of the Gibbs sampler, which uses directly the predictive distributions closed form and therefore does not require sampling block parameters. At each iteration m, the sampler alternates the following two-steps :

1. Draw z (m+1) | w (m) , X, α, G 0 , 2. Draw w (m+1) | z (m+1) , X, β, G 0 .
During the first step, the row memberships update is performed sequentially : each row-cluster membership z i is updated with the other row-memberships z -i and columnpartition w (m) fixed, following p(

z i = k | z -i , X, w (m) , α, G 0 ) ∝      n k n -1 + α p(x i,. | w (m) , X -i , z -i , G 0 ), existing row-cluster k, ( 4.1)
α n -1 + α p(x i,. | w (m) , G 0 ), new row-cluster, (4.2)
where n k is the size of row-cluster k. We emphasize that the parameters Θ do not appear in these formulas, as they are integrated over in the predictive distributions.

In eq. (4.2), the joint predictive distribution of the row p(x i,. | w (m) , G 0 ) is the product of the joint predictive distributions in the block components :

p(x i,. | w (m) , G 0 ) = l p(x (m) i,l | G 0 ), with x (m)
i,l the elements of row i in column-cluster l at iteration m. The joint prior predictive distribution of x (m) i,l is obtained by integrating over the component's parameter :

p(x (m) i,l | G 0 ) = θ p(x (m) i,l | θ) p(θ | G 0 ) dθ.
Because G 0 is a prior conjugate to F , this integral is analytically tractable (cf. Sect. 4.2.3 for the detail in the multivariate Gaussian case).

The posterior predictive distribution in eq. ( 4.1) has a similar definition, with G 0 updated with the observations inside the blocks cf. Sect. 4.2.3).

The second step of the Gibbs-sampler is performed symmetrically on column clusters. Once the maximum number of iterations reached, the row and column final partitions are estimated with the mean of the partitions sampled after burn-in (c.f. Sect. 4.2.4 - §3). In the next subsection we detail how eq. (4.1) and (4.2) simplify with our choice of G 0 .

Multivariate Gaussian case

After the time series preprocessing step, each dataset cell x i,j is a d-dimensional numeric vector produced by the fPCA, that we model with a multivariate Gaussian density. As as conjugate prior, we choose G 0 to be the Normal Inverse Wishart (NIW) distribution with hyper-parameters (µ 0 , κ 0 , Ψ 0 , ν 0 ).

The next paragraphs details the joint predictive distributions closed forms that are necessary to compute the memberships probabilities (cf. proof and computation details in Appendix A.3).

Initialization strategy.

Initializing a Bayesian non-parametric MCMC inference algorithm is usually performed with single-component partition [NBA13 ; RM18]. However, [START_REF] David | Sampling from Dirichlet process mixture models with unknown concentration parameter : mixing issues in large data implementations[END_REF] shows that, for high-dimensional datasets with a high number of components, this strategy can give poor performances and recommends initializing the algorithm with a random partition with more components than the actual number of clusters. However, this number is unknown, and our objective is precisely to avoid its specification. A tempting workaround is to initialize the partitions with one component for each observation. However, this choice is computationally expensive because the membership update has linear complexity in the number of blocks. In practice, as a safety measure, we propose an heuristic, consisting in running the inference process twice. In a first run, the inference is initialized with a onecluster partition ; after this first run completion, the maximum block number sampled during the inference is kept and used as the initial number of components for the second run. If this maximal number is aberrantly high, it can be replaced by an upper quantile of the sampled cluster numbers.

Inferring the final partitions.

In the Bayesian Non-Parametric context, the MCMC-based inference process outputs samples of partitions that are drawn from an approximation of the posterior distribution. These sampled must, in turn, be aggregated over the iterations (usually after a given number of burnin iterations) with a consensus partition estimation (c.f. Sect. 1.3.3 for the clustering case).

In our application, the objective is to estimate the row and column partitions modes ẑ, ŵ. Because the row-partition and column-partition are independent, it is possible to resolve the consensus partition estimation problem independently on each dimension. With (ẑ m ) M and ( ŵm ) M the row and column-memberships sampled during the MCMC inference, the problem is equivalent to solve :

ẑ = arg min Z m d(ẑ m , z), ŵ = arg min W m d( ŵm , w).

Functional Bayesian Non-Parametric Latent Block Model

Consequently, the problem can be solved by applying separately the same method [START_REF] Thomas | Finding the mean in a partition distribution[END_REF] than in the clustering case. The complete inference is summarized in algorithm 6.

Algorithm 6 : FunNPLBM Inference input : Dataset X, n × p × d tensor α, β, G 0 (cf. specification strategies in Sect. 4.2.4 - §1) Iteration number M output : Estimated row-partition ẑ and column-partition ŵ

Initialize z (0) and w (0) (c.f. initialization methods in Sect. 4.2.4 - §2) for m ← 1 to M do for i ← 1 to n do Remove x i,. from the blocks composing its current row-cluster. Compute p(z i | z -i , X, w (m) , α, G 0 ) as defined by eq. (4.1) and (4.2) Sample z

(m+1) i Add x i,.
to its new row-cluster.

for j ← 1 to p do Remove x .,j from the blocks composing its current column-cluster. Compute p(w j | w -j , X, z (m+1) , β, G 0 ) as defined by eq. (

Add x .,j to its new column-cluster.

Compute the partitions averages ẑ and ŵ on the last samples (c.f. averaging methods in Sect. 4.2.4 - §3).

Algorithm Complexity

During the row membership update step, the inference process updates the memberships one at a time. With K (m) and L (m) respectively the current number of row-cluster and column-cluster at iteration m, the update of row-membership z i involves the computations of (K (m) + 1) membership probabilities ((K (m) existing clusters and the new cluster), and, for each membership probability, L (m) block parameters posterior parameters updates (based on Sect. 4.2.3) and predictive distribution density estimations (one for each block in the row-cluster).

In practice, instead of updating the block posterior parameters (µ k,l , κ k,l , Ψ k,l , ν k,l ) with the full block contents X k,l at each iteration and at each membership update, it is more efficient to cache these parameters, and only remove x i,. from its previous block-clusters before the membership probabilities computation, and add it to its new block-clusters after sampling.

With this caching, updating the block posterior distribution parameters of every rowclusters has a complexity reduced from O(npd 2 ) (dominated by the covariance matrix estimation cost), to two operations (removing and adding) with complexity O(pd 2 ).

Because of the matrix determinant computation cubic cost, the complexity of computing the joint predictive distributions for one of the L (m) column clusters is O(d 3 ) , and the membership probabilities computations in eq. ( 4.1) and eq. ( 4.2) have a O(L (m) d 3 ) complexity. This computation must be repeated K (m) times, once for each row-cluster. Overall, sampling one row membership has a complexity O(pd

2 + K (m) L (m) d 3 ).
This row-membership update is performed n times for the row memberships, and the symmetrical operation p times for the column memberships, which sums up to a complexity of O(npd 2 + (n + p)KLd 3 ).

With K = max M K (m) and L = max M L (m) , the global complexity is O M npd 2 + (n + p) K Ld 3 .
The complexity is therefore linear in the datasets dimensions, the number of block components and the iteration number. This expression also highlights the interest to keep the projection space dimension d as low as possible.

Experimental setup

Benchmark and experiments are conducted on a dataset sampled from a known generative model. The dataset is generated by sampling from the distributions N f k,l (t), s 2 where f k,l is a given function and s = 0.02. In the following benchmark and experiments, we work on a dataset of dimension 140 × 140, with unbalanced row cluster sizes (20, 30, 40, 30, 20) and column cluster sizes (40, 20, 30, 20, 30), which amounts to 19600 time series.

The quality of the estimated block partition is compared to the known generative partition, based on several scores : the Rand Index (RI), Adjusted Rand Index (ARI) and the Normalized Mutual Information (NMI). The RI is a popular criterion choice in the clustering domain, which represents the proportion of correctly grouped and separated observations with respect to the observed classes. The ARI is a corrected-for-chance version of the RI that takes into account the probability of getting good RI at random. The NMI is an entropy-based criterion from the information theory literature that estimates the quantity of knowledge a partition gives on another. In addition, we also compare the number of co-clusters estimated by the different algorithms to the true number of co-clusters (5 × 5 = 25).

Benchmark and experiments were run on a 64-bit Ubuntu 18.04 LTS with 32GiB RAM and 12 processors Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz. The studied datasets, data generation scripts and the Scala code are available at the following github repository https ://tinyurl.com/funNPLBM, along with the data simulation method.

Baselines and compared methods

As detailed in the introduction, coclustering methods have been the subject of numerous works and applied to various data types. However, most of the existing models and applications concern unidimensional datasets (i.e., d = 1). The coclustering on datasets containing multidimensional cells (d > 1) has been considered only recently in the literature. Apart from the model-based method from [START_REF] Faicel | Model-based co-clustering of multivariate functional data[END_REF] which does not include a dimension reduction aspect and cannot be applied to large datasets, LBM [START_REF] Charles | The functional latent block model for the co-clustering of electricity consumption curves[END_REF] is the only existing model-based parametric method dealing with this use case.

In addition to FunLBM and FunNPLBM (which we denote LBM and N P LBM for brevity in the following), we consider two coclustering that infer the row-partition and column-partition independently : a bi-dimensional Gaussian Mixture Model (B-GM M ) and a bi-dimensional DPMM (B-DP M ).

In the parametric methods cases, the selection step is performed with a grid-search and the ICL selection criterion [START_REF] Christophe | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF][START_REF] Christine | Estimation and selection for the latent block model on categorical data[END_REF]. This model selection tries every model possibility with a maximum of seven row clusters and seven column clusters, which represents respectively 2 × 7 models for the B-GM M , and 7 × 7 for the LBM . We denote B-GM M M S and LBM M S these parametric methods with model selection. For comparison, we also recorded the performances when directly considering the true block number (the methods are then simply denoted B-GM M and LBM ). In each case, we also tested the effect of keeping the best result among several runs, corresponding to the same number of models tested in the model selection (respectively denoted B-GM M 14 and LBM 49 ). In the non-parametric case, the model selection is natively performed. The concentration hyper-parameters are set to α = 0.1, β = 0.1 and G 0 specified as described in 4.2.4. The baselines are summarized in Table 4.1. For each method, the results are the average scores and median block number over ten runs. These performances are displayed in Table 4.2. The method B-GM M 14 and LBM 49 both give perfect scores for this dataset. Their scores show that when the true number of blocks is known and given enough re-trys, finding the true block partition with the parametric approaches is likely. On the contrary, the methods B-GM M and LBM have the worst performances, which proves that using re-trys is mandatory when using random initialization. This phenomenon comes presumably from a known sensibility of parametric model-based methods to initialization [START_REF] Christophe | Choosing starting values for the EM algorithm for getting the highest likelihood in multivariate Gaussian mixture models[END_REF]. However, a high number of re-trys increases the computation time, as illustrated by comparing the computation times of LBM and LBM 49 . We also note that the performances associated to B-GM M are higher than the ones of LBM , which suggests that the LBM convergence might be harder to reach than the one of B-GM M .

Performing a model selection also increases the computation time, as is illustrated by the computation times of B-GM M M S and LBM M S . The scores of these methods illustrate their capacities to find the true block partitions. However, when the grid-search spectrum is wide, when it is necessary to launch several re-trys, and when the dataset dimensions are high (or a combination of these three conditions), the associated computation cost becomes prohibitively high. In comparison, the N P LBM method only runs once for an equivalent result and in less time. We acknowledge that these scores and runtimes depend on appropriate hyper-parameter specifications, which we discuss in the next sections.

Hyperparameters specification study

The complete set of the method's hyperparameters is composed of the base NIW distribution G 0 , the concentration parameters (α, β), the number of iterations M and the preprocessing parameters : the Fourier basis dimension and the number of PCA axes.

In the unsupervised analysis context, hyperparameters specification remains an active research topic because there is no label to support hyperparameter inference. Consequently, it is not possible to give definitive good choices of values for the following hyper-parameters, which depend on the dataset contents and must be hand-tuned by the experts. These specifications, however, can be based on the knowledge of each hyperparameters impact on NPLBM's behavior, which we illustrate in the following. In each case, we compare NPLBM's performances when one hyperparameter varies while keeping the others equal to given default values : α = β = 0.1, M = 10, a Fourier expression basis of dimension 30 and 3 PCA axes.

Concentration parameters and number of iterations.

In the NPLBM setting (as in the DPM) the prior distribution of the number of components is an increasing function of the concentration parameters : the higher the concentration parameters the higher the probability of producing high numbers of components, without data knowledge. Because the whole method is symmetric on the dataset rows and columns and because the experiment dataset dimensions n and p are equal, we consider the case α = β.

As shown in Fig. 4.1, the concentration parameters effects are negligible for this dataset and only have an impact for extreme values. When the concentration is extremely high (> 10 12 ), the number of components is highly overestimated. When its value is on the contrary extremely small (< 10 -10 ), only one-cluster partition is inferred. This small impact of α comes presumably from the high separation of the components in the time series high-dimensional observation space. This separability is simulated for this experiment but is consistent with what we observe in practice. This separation also explains the small number of iterations needed for convergence (here, less than 4 for 10 -6 ≤ α ≤ 10 10 ) and the high stability of the MCMC chain.

For a given dataset, if the values of α and β seem to have a strong influence on the inference results, we advise to add a Gamma hyper-prior assumption for α and β and estimate their values during the inference algorithm. This strategy is a common practice in the DPM setting [START_REF] Mike | Hyperparameter estimation in Dirichlet process mixture models[END_REF]. However, it implies to study specification strategies for the Gamma distribution parameters. The Fourier Basis dimension and the PCA axes numbers influence both the performance of NPLBM and the trade-off between sparsity and quality of time series representation.

Experimental setup

Fig. 4.2 shows the effects of under-estimating or over-estimating the number of PCA axes. Low numbers of PCA axes (<3) are associated with poor scores and few block components due to poor representations of the time series, which are too close in the projection space. If the number of axes is too high (>5), the high dimensionality exaggerates the time series separation and the component number is overestimated, which explains the sharp decrease of the ARI and NMI.

In our use cases, we observed that three PCA axes lead to the most interesting results. In Fig. 4.3, we observe that, with a fixed number of 3 PCA axes, high polynomial basis dimensions (> 50) are correlated with poor scores, presumably because of lower time series representation quality (reflected by the low variance explained score). On the contrary, when this basis dimension is low (<10), the 3-axes PCA adequately represents the information and the variance explained is high (>= 0.97). However, in this case, the Fourier basis dimension is too low to adequately represent the time series, which explains the poor scores. The studied datasets, data generation scripts, and the Scala code used for the benchmark and the experiments are available at the following github repository https ://tinyurl.com/funNPLBM, along with the data simulation method. In the next section, a real-case situation is studied and illustrates the method's interest for Advanced Driving-Assistance Systems validation.

ADAS validation

Validating an ADAS is a complicated task, that is partly addressed with on-track tests. This practice has been gradually replaced with numerical simulations. The first motivation is the reduced testing cost compared to real physical tests which requires specific track infrastructure, equipment management, and significant human intervention. One digital simulation is estimated 10,000 times cheaper than its physical counterpart. The savings achieved through the use of digital simulation add up to millions of euros.

The second motivation comes from the fact that numerical validation remove uncertainties due to physical conditions (e.g. sensor error, initial conditions, weather conditions, ...). As a consequence, we are able to test the driving software independently from uncontrolled external factors as .

Another major disadvantage of physical tests is the prohibitive number of test sample to check for high reliability of the system. A validation objective may be the assessment of vehicle incident odds (e.g. < 10 -8 casualties per hour). With a classical sampling method, estimating such probability would require real driving sessions of billions of kilometers.

Even if such a large amount of real-life data were available, as is the case in some data science applications, there would be no guarantees of the data quality or value. In our case, this value lies in the specific driving situation in which to test the system behavior. Indeed, some situations are rarely observed in reality such as emergency procedures, like the one considered in this numerical application.

Use case description

This section illustrates the use of the coclustering approach to the Emergency Lane Keeping (ELK) assistance system validation. In a straight lane scenario, the vehicle under test (called ego) is drifting towards an oncoming car on the other lane (c.f. Fig. 4.4). The ELK system is expected to detect the drifting, the oncoming car, and to put the vehicle back to its lane center with an emergency maneuver.

For our dataset, we consider n = 400 simulations that were generated by simulating the ADAS according to a particular design of experiments (varying ego speed, the drift angle, . . . ). The objective of this analysis is to find characteristic operating modes of the intelligent system described as p = 22 time series and to discriminate relevant groups of correlated variables.

Results

The final coclustering is the block partitions mean (c.f. averaging methods in Sect. 4.2.4) of 10 samples obtained after a burn-in of 10 iterations and is composed of 6 row-clusters and 13 column-clusters. With color indicating block membership, Fig. 4.5 shows the global coclustering structure and Fig. 4.6 an extract of the block contents. The row-clustering is also insightful : each row-cluster correspond to well-separated driving behaviors. This separation is best seen in Fig. 4.6 (top-right) that shows the following driving behaviors : 1) ego is drifting left and the ELK system fails (light green) ;

2) the symmetric behavior on the right (dark green) ; 3) the ELK system corrects the car trajectory (light orange).

Finally, the three other row-clusters (regrouping the remaining 5% of the observations) are composed of outliers simulations, with driving behavior displayed on Fig. 4.6

(bottom-right). In this situation, the oncoming car is correctly detected, ego heading is changed accordingly, but still the collision cannot be prevented. In conclusion, this co- clustering partition shows that the approach has correctly discriminated uninformative signals while creating meaningful clusters of features and clusters of simulations. From this information, we can visualize the variety of driving behaviors that compose the datasets. Moreover, we can understand them from the variable perspectives which was the original objective of the application. The next step is to link the driving behavior to the control logic parameters and, if need be, optimize them to reach the performance objectives.

Conclusion

This chapter describes FunNPLBM, a Bayesian non-parametric based method which addresses the problem of coclustering multivariate time series without prior specification of the model dimension. This work offers a novel definition of the NPLBM model, the description of an adapted collapsed Gibbs sampling, and overall the first BNP coclustering method dedicated to multivariate time series data analysis. This model enables regrouping redundant signals and discriminating uninformative ones. Moreover, the methodology provides a two-dimensional overview of a multivariate time series dataset to the users.

The specification of the different hyperparameters (concentrations, NIW, pre-processing) is discussed and experimented on a simulated dataset, which shows the model performance and gives insights on the most interesting hyper-parameter specification strategies.

Finally, the method is applied to a real-case dataset from the advanced driving-assistance system validation domain. In this application, FunNPLBM proved its ability to simultaneously discriminate groups of signals and to produce meaningful driving behavior clusters. These results show the usefulness of the model and of the model selection strategy in concrete modern challenging use cases.

The FunNPLBM approach was applied here to an autonomous driving context, however we are confident that it can be used in many other domains. A first interesting perspective would be the extension of a variational inference algorithm, a strategy that has been proven interesting in [START_REF] David | Variational inference for Dirichlet process mixtures[END_REF]. Another extension would be to define similar BNP concept to datasets of higher dimension and develop tensor coclustering methods for higher-dimensional functional problems.

Finally, the combination of FunNPLBM and CLBM from Chapt. 3 brings us to the next contribution titled Bayesian Non Parametric Multi-Clustering, that models each columncluster with a specific FunNPLBM component.

Multivariate Time Series Multi-Coclustering 5

This chapter introduces a block structure model that combines two layer : a Multi-Clustering structure that groups together the variables with common row-partition, and a coclustering layer that associate one specific NPLBM structure to each columncluster.

Introduction

In a multivariate context, the possible relations between variables introduces additional difficulties (c.f. Sect. 1.4.3). In real-life use cases, it is likely that each variable is associated to a different partition of the observation space (called row-partition).

In the following, we assume the presence of a variable partition (called redundant column-partition), such that variables in a redundant column-cluster share the same row-partition. In our use cases, this assumption seems rational, as we know that many variables are related (e.g., position and acceleration, ADAS activation and braking). The Multi-Clustering [START_REF] Juhua | Subspace multi-clustering : a review[END_REF] designates a set of methods that infers this multiple-partitions representation of a dataset, as illustrated in Fig. 5.1 -b).

When in a parametric model-based context [MV19 ; Van20], these approaches requires to know the partitions sizes a priori (rarely true in practice), or to perform a model selection step. However, the combinatorial nature of Multi-Clustering makes the model selection a complex task (c.f. Sect. 3.2.4). Heuristic approaches could be considered at the cost of assumptions on model structure (e.g. with greedy strategies). The nonparametric approach circumvents this issue by natively performing a model selection during inference.

A non-parametric Bayesian Multi-Clustering model has been developed by [START_REF] Yue | Variational inference for nonparametric multiple clustering[END_REF] for continuous datasets. In this model, the rows belonging to a block follow an independent multivariate distribution (c.f. Fig. 5.1-b). It must be emphasized that, in this work, the block component distributions models the rows and not the cells as in the coclustering case. The reason for this is that each variable follows a different distribution, which means that the overall model to estimate has a p × d dimension, and is invariant to the number of column-cluster.

This Multi-Clustering method groups variables according to their row partition but cannot regroup variables with similar distributions. Moreover, estimating the block distributions parameters becomes problematic when the column cluster sizes are important (e.g., in our case, high dimensional covariance matrices estimation). We propose to model each redundant column cluster with a coclustering structure to deal with this limitation. In each coclustering structure, the column-partition is called the correlated variable partition, and regroups variables with common distribution..

The model-based coclustering [GN13] (Fig. 5.1 -a)) infers one row-partition and one column partition. This model assumes that all variables share the same row partition and that there exists groups of variables with common distribution. We propose to model each redundant variable cluster with a coclustering structure, which assumes that some variables sharing the same row-partition can also share the same distribution. This assumption seems also natural, because some of the simulated variables are physically correlated (e.g., car speed and wheel rotation speed). In this paper, we propose the Functional Non-Parametric Multi Coclustering (FunMCC) model (c.f. Fig. 5.1-c) that combines the Multi-Clustering and coclustering approaches and produces the best of both worlds : the Multi-Clustering layer infers the redundant variable partition and extracts the true sparse number of row partition while the coclustering layer reduces the parameter set dimension and regroup correlated variables.

To the best of our knowledge, there is only one comparable work by [START_REF] Tomoki | Multiple co-clustering based on nonparametric mixture models with heterogeneous marginal distributions[END_REF] in the bayesian non-parametric setting, which does not apply to multivariate time series.

Multiple

Coclustering of multivariate time series

FunMCC model definition

In the following, the dataset X is obtained with interpolated periodogram fPCA representation (c.f. Sect. 3.2.1) as in the previous chapters.

As in Chapt. 3, v denotes the Multi-Clustering partition, i.e., the redundant partition, and Z the n×H row-partitions indicator matrix, with H the number of redundant clusters. As in Chapt. 4, the column-partitions inside each coclustering structure (i.e., the correlated partitions) are denoted (w h ) H (with an additional superscript denoting the coclustering structure index). With this notation, the model is formally defined by :

x i,j | {v j = h, w h j = l, z h i = k, θ h k,l } ∼ N (θ h k,l ), θ h k,l ∼ G 0 , v j ∼ M ult(η), w h j ∼ M ult(ρ h ), z h i ∼ M ult(π h ), η j (r) = r j j-1 j ′ =1 1 -r j ′ , r j i.i.d. ∼ Beta(1, γ) , γ ∼ Gamma(a γ , b γ ), ρ h j (s h ) = s h j j-1 j ′ =1 1 -s h j ′ , s h j i.i.d. ∼ Beta(1, β h ) , β h ∼ Gamma(a β , b β ), π h j (t h ) = t h j j-1 j ′ =1 1 -t h j ′ , t h j i.i.d. ∼ Beta(1, α h ), α h ∼ Gamma(a α , b α ),
where each memberships proportions vector η, (ρ h ) h and (π h ) h follows a stick-breaking construction scheme [START_REF] Jayaram | A constructive definition of Dirichlet priors[END_REF]. Each concentration parameter follows a specific Gamma hyperprior.

The role of this additional hyper-prior is to grant the ability to infer and adapt the concentration parameter value to each coclustering structure, which are likely to contain different numbers of block components.

The block distributions are multivariate normal and G 0 is the associated conjugate Normal-inverse-Wishart prior. The complete hyperparameter set is noted χ.

It can be noted that, with this definition, constraining the correlated partitions to contain only one column-cluster makes the model equivalent to a Bayesian Non Parametric extension of the CLBM method described in Sect. 3. In addition, constraining the redundant partition to contain only one column-cluster makes the model equivalent to the FunNPLBM method described in Sect. 4. For these reasons, the FunMCC

Inference

We propose a stochastic iterative inference based on the following three-steps Gibbs sampler : 1) during the Multi-Clustering layer inference step, the redundant column partition v is updated given Z ; 2) in the Coclustering layer inference step, the coclustering partitions Z and (w h ) H are updated given v 3) Each concentration parameter is updated.

First step : Multi-Clustering layer

In the first step, the memberships (v j ) p are sampled sequentially according to :

p(v j | v -j , Z, x .,j , χ) ∝        p h p -1 + γ p(x .,j | z h , χ), existing cluster h, (5.1) γ n -1 + γ p(x .,j | χ), new cluster, (5.2)
with v -j = {v i : i ̸ = j}. In Eq. (5.1), the joint prior predictive distribution p(x .,j | z h , χ) reduces to a product of multivariate t-student densities [START_REF] Andrew | Bayesian data analysis[END_REF] thanks to an appropriate choice of conjugate prior G 0 (c.f. Sec. 4.2.3). In Eq. (5.2), p(x .,j | χ) is estimated with an univariate Dirichlet Process Mixture (once per variable, prior to the inference of the FunMCC), which also produces the row-partition ẑj associated with variable j and to a possible new cluster.

Second step : Coclustering layer

In the second step, an NPLBM [START_REF] Edward | Nonparametric bayesian biclustering[END_REF] is inferred on each sub-dataset (X h ) h defined by the partition v with a Gibbs algorithm that simulates p(z h , w h | v, X h , χ). The inference is performed in two symmetrical steps : first, z h is updated with fixed partition w h , then w h is updated with fixed z h . Because these steps are symmetrical, we only detail the column membership update in the following. For each variable j in sub-dataset X h , w h j is updated by sampling from p(w

h j | w h -j , z h , X, χ) ∝                p h l p h -1 + β K h k=1 p(x h k,j | x h k,l , χ), existing cluster l, ( 5.3
)

β p h -1 + β K h k=1 p(x h k,j | χ), new cluster, (5.4)
with w h -j and z h the column and row partition indicator vectors without variable j, p h l the size of the correlated column cluster l, x h k,j = {x i,j : v j = h, z h i = k}, and p(x h k,j | x h k,l , χ), p(x h k,j | χ), respectively, the joint posterior and prior predictive distributions in block (k, l). After the row and correlated partition updates, the concentration parameters α h and β h are updated.

Third step : Concentration parameters update

Following [START_REF] Mike | Hyperparameter estimation in Dirichlet process mixture models[END_REF], the concentrations parameters γ, α and β are updated based on a Gamma prior. This update process is the same for every concentration parameter with their respective Gamma prior parameters. We also emphasize that the coclustering structures share the same hyperprior values : the concentration parameters α = (α h ) H follow the prior Gamma(a α , b α ) and the concentration parameters β = (β h ) H follow the same prior Gamma(a β , b β ). In the following we give the detail for γ's update.

Multiple Coclustering of multivariate time series

Given H the current number of coclustering structures, γ the current value of the columnpartition concentration partition, p the number of variables, and (a γ , b γ ) the Gamma prior parameters, the update process is : The same process is applied to (β h ) H and (α h ) H , with corresponding priors, number of elements and number of clusters. In the following we detail how the appropriate choice of conjugate prior G 0 prior enables the tractable computation of the prior and posterior predictive distributions used in Sect. 5.2.2 and Sect. 5.2.2.

(i) Sample y ∼ Beta(γ + 1, p) (ii

Implementation

The implementation details are similar to the ones of the FunNPLBM discussed in Sect. 4.2.4 :

-We also use a conjugate NIW prior with a zero-valued mean µ 0 , a precision matrix Ψ 0 inferred from the dataset, and minimal values of κ 0 and ν 0 -The initial state is a one-cluster partition, and a second inference may be performed with block numbers higher than the one sampled in the first run.

The final partitions inference step is different than for the FunNPLBM method.

Inferring the final partitions.

In the Bayesian Non-Parametric context, the MCMC-based inference process outputs samples of partitions that are drawn from an approximation of the posterior distribution. These sampled must, in turn, be aggregated over the iterations (usually after a given number of burnin iterations).

In the FunNPLBM model, the row and column final partitions can be independently obtained by applying a consensus partition estimation [START_REF] Thomas | Finding the mean in a partition distribution[END_REF] separately on each dimension.

In our cases, the objective is to obtain the final sets of partitions modes v, (ẑ h ) H , and ŵh ) H . Although it is straightforward to obtain the consensus redundant partition v with [START_REF] Thomas | Finding the mean in a partition distribution[END_REF], the consensus on the other dimension partitions are much complex to estimate. This is caused by the label switching effect that makes the model likelihood invariant against component indices re-ordering.

In other words, if one wants to find the consensus correlated partition corresponding to redundant cluster h in v, one must first find, in every sampled FunMCC results, the FunNPLBM structure that match with the one associated with the redundant cluster h.

This problem can be seen as a consensus partition estimation with a hierarchical partitioning constraint.

Here are several workarounds solutions :

-A first method consists simply in keeping the last sample, which assumes that the MCMC has converged.

-It is also possible to select the sample that maximizes the model likelihood. This method seems the most natural, but is not a consensus estimation.

-A decoupled solution consists in first running the SEM-Gibbs to estimate the consensus redundant partition, then infer separately the FunNPBLM structure with fixed redundant clusters, and finally to estimate the consensus coclustering solutions of each FunNPLBM separately.

-As a coupled extension of the previous solution, a consensus partition estimation can be performed after each step (Multi Clustering and Coclustering layer) at each iteration. In that case, each step is based on the consensus partition(s) estimated in the previous.

Including this final partition estimation step, the complete inference process is summarized in Algorithm 7.

Algorithm Complexity

Before the main FunMCC inference, one DPM is inferred on each variable x .,j independently, in order to estimate the expression p(x .,j | χ) used in Eq.(5.2) and to obtain the row-partition ẑj associated with variable j "alone", i.e., the partition associated to a newly discovered cluster containing only the variable j. Each DPM has a complexity O(M nd 2 ), with M the iterations number (with same value than for the main inference), n the observation number, and d the time series representation space dimension. This inference is performed on each variable, therefore, which amounts to a O(M npd 2 ) complexity. These estimations are only performed once, then can be cached.

After this preparation step, the main FunMCC inference process begins, alternating the two phases : the Multi-Clustering layer update and the coclustering layer updates.

During the Multi-Clustering layer inference, the redundant cluster memberships are updated sequentially following 5.2.2, with a fixed value of Z. The update of one membership v j involves the computation of H + 1 membership probabilities (with H the number of existing redundant clusters plus one new cluster). For a given existing redundant cluster h, computing the redundant membership probability of column x .,j involves, in turn, to compute one joint prior predictive for each of the K blocks (one for each row-cluster) that composes the column-cluster h. These updates, detailed in Sect. 4.2.3, consists in combining the n observations by groups into K means, and K covariance matrices, with an overall O(nd 2 ) complexity (dominated by the covariance matrix estimation cost).

As opposed to the FunNPLBM case (cf. Sect. 4.2.4), it is not possible to cache the updated hyper-parameters values, because they are only computed once at each iteration (for this particular membership distribution computation) given a fixed value of Z that changes during the coclustering update set.

After having updated the prior hyper-parameters, the membership probability is finally obtained by computing p(x .,j | z h , χ) the joint predictive distribution of variable j. This involves K h computations (with K h the current number of row-clusters in redundant cluster h), each with complexity d 3 . This operation is performed once for each of the H existing cluster, and for each of the p variables, which amounts to a complexity

O(p H(nd 2 + Kd 3 )), (5.6)
where H = max m H (m) and K = max m,h K (m),h are the maximum of the cluster number values sampled during inference.

The coclustering layer inference consists in updating H FunNPLBM structures, therefore with complexity (derived in Sect. 4.2.4) O H npd 2 + (n + p) K Ld 3 , with L = max m,h L (m),h .

The multi-clustering and coclustering layer updates are repeated for M iterations. The final complexity is, after simplification :

O M H npd 2 + (n + p) K Ld 3 ,
with M the number of iterations. This complexity is linear in the datasets dimensions, the number of block components and the iteration number.

It may be noted that the coclustering structures can be updated in parallel, which would cancel the linear complexity in H. It is also possible to parallelize the independent DPM that are inferred before the main FunMCC inference to estimate the distributions p(x .,j | χ).

Experiments

In the following section we highlight the interest of using FunMCC by comparing the impact of uninformatives or misleading variables on several block clustering baselines. The experiments dataset is composed of observations generated from a ground truth row-partition, that we seek to estimate. These observations are described by several groups of variables sharing common distribution. This ground-truth dataset is therefore a coclustering structure.

This ground-truth information is "corrupted" by adding a proportion of other variables : uninformatives (generated from a one-cluster partition) or misleading (generated from another row-partition than the ground truth). These proportions are respectively denoted p u and p m . In the ground-truth coclustering structure, 30 variables and 150 rows are distributed in 5 row-clusters and 3 column-clusters. The uninformatives variables are simulated with one-partition variables, and the misleading structure with another coclustering structure with 4 row-clusters and 3 column clusters. The ground-truth and misleading row-partitions are randomly and independently generated. With the proportions p u and p m varying from 0 to 200%, this setup generates 7 datasets, each one corresponding to the ground truth and given uninformative or/and misleading proportions.

The objective is to find the ground-truth without supervision, and to discriminate and separate the uninformative and misleading information.

The elements belonging to a given block (i.e., to one particular redundant cluster h, and one block (k, l) in the coclustering structure h) are generated by sampling from the distributions N f k,l (t), s 2 where f k,l is a given prototype function and s = 0.02. An instance of such dataset is displayed in Fig. 5.2, with equal proportions of ground-truth, uninformative and misleading variables. This figure also represents three prototypes used to generate the associated block contents.

Observa ons

Variables Ground-truth Uninforma ve Misleading We compare FunMCC's performances to baselines non-parametric model-based approaches, that match to our specifications. We consider a Dirichlet Process Mixture model (DPM) that infers a row-partition without a column-partition and represents the default strategy without variable discrimination. The second baseline is a non-parametric coclustering method (CC) [START_REF] Edward | Nonparametric bayesian biclustering[END_REF] displayed in Fig. 5.1, that adds a column-partition.

As alternative Multi-Clustering solutions, we consider three methods derivated from our proposal FunMCC :

-Decoupled Multi-Clustering based on DPM (DPDM) consisting of one DPM on the columns, then one on the rows.

-Decoupled Multi-Clustering based on Coclustering (DCC) that first infers the redundant partition with a single global coclustering, then infer a coclustering structure in each of the obtained redundant cluster.

The results, displayed in Table . 5.1, show a gradual performance increase ranging from DPM to FunMCC. The DPM, which only creates one row-partition, is the most impacted by the presence of uninformative and misleading variables. The DPM performances are gradually worsened when the p m (the misleading proportion) increases. On the contrary, even a small proportion p u (the uninformative proportion) reduces DPM's performances dramatically. This is presumably induced by the non-parametric setting that selects the most likely number of row-clusters, and favours one-cluster partitions when adding uninformative signals. On the contrary, the presence of misleading information does not reduce the number of row-clusters, but corrupts their contents, resulting in lower performances with an approximatively correct number of row-clusters.

The CC method is not impacted by the presence of uninformative signal, but tends to overestimate the row-cluster number in presence of misleading variables. This overestimated K is associated to the row-partition that results from crossing the ground-truth row-partition (5 row-clusters) with the misleading row-partition (4 row-clusters), resulting in a 8-clusters row-partition, as approximated by CC. The Decoupled methods (DDPM and DCC) have similar performances than CC, with a nearly null impact of p u , but with a greater impact of p m . This effect is linked to a low-quality redundant column-clustering, due to the decoupled aspect, that in turn impacts the row-partition quality. With the decoupled methods, the Multi-Clustering is the result of two steps : a column-partition inference, then one row-partition inference in each column-cluster. However, the row-partitions information, obtained in this second step, cannot be used in the column-cluster construction.

On the contrary, the pure non-parametric Multi-Clustering uses the row-partitions information to construct the redundant column-partition, which explains a better columnpartition quality, and, as a consequence, a stronger stability against the presence of misleading variables. We observe a light performance increase with the addition of a coclustering structure layer (our proposal FunMCC), that regroups variables with similar distributions. Our intuition is that the coclustering layer acts as a shrinkage method that reduces the parameters number and the global model complexity while keeping the same overall density representation capacity. This phenomenon presumably increases the Multi-Clustering layer quality, in addition to enhancing the model interpretability.

Applications

Emergency lane Keeping Assist

After experimenting FunMCC on a synthetic dataset, we apply it on a real ADAS validation dataset. In a straight lane scenario, the vehicle under test is drifting towards an oncoming car on the other lane. The ELK system is expected to detect the drifting, the oncoming car and put it back to its lane center with an emergency maneuver. The simulation generates a dataset of n = 500 described by p = 150 temporal variables, totalling 75000 time series. We emphasize that these simulations are generated with a simulation black-box that faithfully recreates the real-life driving conditions, and are not produced by the generative model proposed in this article.

The objective is to infer a synthetic structure that highlights interesting driving patterns and discriminate relevant groups of sensors, in order to isolate and understand the ADAS activation contexts.

One row = one observa on

One column = one variable A first interesting results is that the non-informative variables variables are grouped together in the red and orange column-clusters. This group contains normally uninformative variable such as the car length, the time from start, but also use case specific (e.g., road radius in a straight lane scenario, headlights activation while in broad daylight, . . . ).

Several other column-clusters are particularly interesting :

-The leftmost green column-cluster contains orientation and car direction variables (e.g., the car orientation, the wheel angle, the steering wheel angle, . . . )

-The next ochre column-cluster on the figure contains position variables : lateral position, index and type of the current lane, distance to the road center, . . .

-The next small blue column-cluster contains ADAS activation variables, which are categorical time series. These variables can also indicates the activation "mode" (for instance, activation on the left or right).

The column-clusters composition exhibits one outstanding feature of this approach : the column-cluster acts as meta-variables, that regroup the meaning of several others. This process acts simultaneously as a variable selection and dimension reduction step, and helps the interpretation.

Another interesting information that appears on Fig. 5.3 is the hierarchical structure between the row partitions. Because the simulations are sorted by row-partitions, and based on the order of the row-partitions, the successive sorting creates a hierarchical clustering. The specification of the row-partitions order must be specified by the user. This problematic is specific to the Multi-Clustering case, and is discussed in the perspectives (Sect. 6.2.2).

The Fig. 5.4 illustrates the content of a column-cluster, and shows the car trajectories (surrounded by a red rectangle in Fig. 5.3) distributed in the row-clusters, where several scenarios are represented : the sub-graph I shows a late ADAS activation scenario, subgraph II regroups cases when ELK did not activate at all. In sub-graph III, the systems correctly change the car trajectory and prevent the collision. The last graph IV shows a very small cluster of anormal driving behaviors detected and grouped by FunMCC.

Finally, the understanding of the Multi-Coclustering result structure can be enhanced by representing the variables in their own space and in their column-partition. This Based on these results, the next steps for the ADAS system developer consists in relating the partitions to the input simulation parameters and assessing the ADAS compliance with its specifications.

Conclusion

This chapter describes a new Bayesian non-parametric based method designed for the exploration of multivariate time series datasets produced by driving simulations. This solution infers a multi-level dependency structure that highlights the relationships between sensors and discriminates driving simulations patterns. This probabilistic method can natively be used for anomaly detection based on probabilistic predictive intervals. The next step in the ADAS validation pipeline is to use the multivariate time-series multicoclustering results as labels for a supervised analysis of the simulation parameters, in order to understand and explain ADAS failure and success conditions.

Conclusion

Conclusion And Perspectives 6

Dictionary-based time series clustering

In the context of unsupervised classification of regime-changing time-series, this chapter proposes a dictionary-based method that consists in three steps : automatic segmentation of each time-series, regime dictionary construction, and clustering of produced categorical sequences. SDLHC shows good results when applied to time-series complying to the regime construction assumption, and is competitive with other state-of-the-art methods in this case. The ability to handle unequal-length time-series, a-synchronized time-series, and time-series exhibiting asynchronous regimes are its best assets.

One weak point of this method is the complexity associated to the third step, which is O(n 3 ). This complexity makes the method quite slow when the number of time series becomes important. In that case, alternatives methods should be considered, e.g., using a partition-based clustering on histograms, at the cost of losing the time-dependency information.

One other weak point is that this method can only be applied to univariate time series, which limits its application to specific use cases. In the ADAS validation context, this means that the user is limited to the exploration of one simulated variable. The multivariate extension can be addressed in several ways.

Joint Multivariate Segmentation

The easiest extension to multivariate clustering consists in assuming that, for a given MTS observation, every dimension share the same segment cut-points (i.e., the same univariate hidden logistic process). Under this assumption, it is possible to perform a multivariate joint segmentation, similar to [START_REF] Faicel | Time series modeling by a regression approach based on a latent process[END_REF] with multidimensional segments. This method has been described in [START_REF] Faicel | Joint segmentation of multivariate time series with hidden process regression for human activity recognition[END_REF] by the same authors.

In this case, the dictionary construction step can be extended in a multivariate GMM that regroups multivariate segments, and the produced driving patterns are also multivariate. The MTS are, in turn, recoded as univariate symbol sequences, which leaves the final step (univariate symbol sequence clustering) unchanged. The overall transformation and clustering can be considered a dependent approach, as every dimensions are considered jointly.

The joint multivariate segmentation can be suited for dataset composed of a small set of correlated variables (e.g., position and speed variables). However, this method is not suited to deal with a large number of variables, when the synchronous cutpoint hypothesis becomes more and more unrealistic.

Independent Univariate Segmentation

A short-cut to the previously described limitation is to independently segment every dimensions. In this case, three alternatives can be considered.

Cut-points Intersection

A first solution consists in taking the union of the cutpoints cut-points of every dimensions, and segments the MTS in smaller sub-sequence.

After this step, the dictionary can perform, as in the previous variant, a multidimensional segment clustering that produces multivariate patterns. The MTS are recoded in univariate symbol sequences, and the final hierarchical clustering remains unchanged.

However, intersecting the cut-points of every dimensions can lead to an excessive subsampling of the segments, and to an important loss of information.

For instance, given a p-dimensional MTS with length T , and assuming that the pdimensions have, on average, K distinct cut-points, the union of these cutpoints generates a pK cutpoints. As a consequence, the average length of the d-dimensional subsequences is T pK . With plausible values (T = 1000, p = 100, k = 10), the sub-sequences becomes singleton sub-sequences, which cancels the dimension reduction effects of the segmentation step, and makes the dictionary construction step (a GMM on a R d = R 100 space) unfeasible.

Symbol Sequences Intersection

Instead of taking the union of the cut-points, which may greatly over-estimate the number of segments, it is possible to apply the union operation on the uni-dimensional symbol sequences.

In this variant, after the univariate segmentation of each dimension, the dictionary reconstruction step is applied independently on each dimension. As a results, each MTS is recoded as a multivariate symbol sequence.

After this phase, each unique combination of univariate state is associated to a new multivariate symbol. For instance, given an MTS T 1 (t) represented with a two-dimensional symbol sequence ((a,a,a,b,b,b),(d,d,e,e,f,f )), the resulting representation would be ( A,A,B,C,D,D This method implies to construct a new dictionary that associates a unique symbol to each multidimensional symbol state. As a result, each MTS can be represented with a univariate symbol sequence, and the final clustering step remains unchanged.

The drawback of this method is the same as the previous cut-points union method : when the number of variable is high, the length of the sub-sequences associated to each multivariate state decreases, and the number of possible multivariate state increases sharply.

The problem is the same if the dictionary construction is performed on the overall set of segments from each dimension rather than dimension-wise (with D replaced by the size of the overall dictionary).

Independent Symbol Sequence Comparison

Instead of combining the multi-dimensional symbol sequences to obtain univariate strings, this variant proposes to keep a multi-dimensional string representation and to perform the last clustering step based on this multi-dimensional representation.

In that case, the dictionary method remains unchanged and deals with univariate segments, either applied dimension per dimension or to every segments from every dimensions), as in Sect. 6.1.2. As a result, each MTS is recoded as a tuple of symbol sequences.

Based on this representation, the MTS are compared by summing their dimension-wise difference string dissimilarity, e.g., the Levenshtein distance. This setting implies to compute the dimension-wise sub-sequences dissimilarities and the Levenshtein distance matrix, which adds a linear complexity in d to the step complexity. This solution can, up to the final hierarchical clustering, be parallelized or distributed per dimension, which makes it scalable and interesting for industrial applications.

A major drawback of this solution is that each dimension has the same weight, which makes the overall method sensible to uninformative or misleading variables. Moreover, summing the distances in this high-dimensional space produces strongly separated observations.

In order to solve this problem, a perspective would be to consider weighted version of this multi-dimensional Levenshtein Distance.

Coclustering extension

Another approach to extend the method to the multivariate case is to consider a coclustering extension. Using the fact that the block component are uni-dimensional distribution (shared by every variables inside a variable cluster), our first idea was to develop a model-based version of SDLHC that could be extended to the multivariate case with the coclustering framework.

In order to obtain this model-based extension, it is necessary to replace the last step of SDLHC (the symbol sequence clustering), with a model-based clustering method. After this replacement, the objective was to collapse the entire method and to perform a joint inference of the segmentation, dictionary construction and symbol sequence clustering in a single inference process.

However, the development of a model-based categorical sequence clustering has proven difficult, partly because of the complexity of estimating a consensus categorical sequence. This limitation led us to modify our strategy, and adopt a specific representation for the multivariate case.

Multi-Coclustering

The Multi-Coclustering method can be extended in several ways. As a model-based method, it can natively be used for anomaly detection (e.g., by considering a componentwise likelihood ratio). It can also serve as a supervised classification method, by simply constraining the redundant, correlated and row partition values. In that case, the inference would lead to a multiple coclustering density estimation, that would enable a multiple classification of a given observation. In a semi-supervised application, the method can be used to infer missing value in one coclustering structure based on the others.

FunNPLBM can be useful in every domains that deal with datasets containing multiple row-partitions and correlated temporal variables. For instance in industrial contexts for predictive maintenance based on multiple sensors, in health for ECG and biological signals data analysis, in finance for stock trade data analysis.

Another natural extension would be the online multiple coclustering applied to the realtime detection of driving states and the online anomaly detection.

Method Scalability

One drawback of the inference method is that, as is the case for the DPM, the memberships (of every dimensions) are updated one by one given the other observation and other memberships. Because of this sequential inference, it is not straightforward to parallelize the inference at the individual level. This can cause some inference issues when the datasets dimensions are too important. However, following recent advances in the field [START_REF] Khadidja | Dirichlet process mixture models made scalable and effective by means of massive distribution[END_REF], it is possible to distribute the computation when the dataset is important. The DPMM distribution consists in running independent DPMM on sub-datasets, and then regroup the obtained clusters. In the next works we suggest to develop an extension of the work of [START_REF] Khadidja | Dirichlet process mixture models made scalable and effective by means of massive distribution[END_REF] by distributing the FunNPLBM, and the FunMCC.

Representing the variables relationships

The MCC method outputs several row-partitions, but also a nested column-partition : based on the observation clustering, then on the distributions. An interesting perspective for the user is to visualize the relationships between the variables, which can help the results interpretation.

A Partition of variables

This relationship can be represented with a partition-based clustering, as illustrated in Fig. 5.5. In this example, an MDS projection is used based on the combination of two dissimilarities : the row-partition dissimilarity, and the time series distance. This representation can help interpreting the contents of the column-clusters (e.g., associate a column-cluster to Ego's direction, or lateral position).

However, the weights of this combination have been set manually, which is not a straightforward task. In the following we consider developing alternative representations of the variable.

One interesting candidate can be the dandelion plot, that allows a multi-level categorical grouping based on hierarchical dissimilarities. It is also possible to simply ignore the partition similarity information and construct the row-partitions on the only basis of the distributions.

Choosing the row-partition ordering

In the coclustering framework, the user is interested in visualizing the inferred block structure, which is obtained by re-ordering the dataset matrix by group. For instance, this re-ordering produces the grid visualization shown in Fig. 1.13 b), or in Fig. 3.2 a).

With MCC (and more generally with the Multi-Clustering methods), the re-ordering of the rows is performed hierarchically : based on z 1 , then z 2 . . . until z H . Because it is hierarchical, this re-ordering depends on the redundant cluster order (e.g., swapping z 1 and z 2 will modify it). However, these different visualisations are not equivalent to the user, because some row-partitions are more interesting than the others. In addition, several partitions are included into others, which means that they are the result of a hierarchical subdivision, similar to to a hierarchical clustering dendrogram. Choosing this ordering manually is not straightforward (especially when the number of redundant clusters H is large) and, in the current implementation of FunMCC, depends on expert knowledge. An interesting extension would consist in inferring automatically this ordering. For instance, assuming that some row-partitions are nested into each others, a tree representation might be able to represent the dependencies between row-partitions. A score that measures the mutual partition information (e.g., NMI, ARI, . . . ) or the partition inclusion seems an interesting candidate for this construction. This hierarchical representation adds an interesting layer of information : the dependency between row-partitions.

Time Series Representation Alternatives

In this thesis contributions, the time series are transformed with a PCA on the coefficient of an Interpolated Log-Periodogram representations. Alternative time series representations can be considered, with a strong incentive on keeping a final representation dimension low (cf. inference complexity description in Sec. 5.2.4).

First, it is possible to replace the trigonometric function basis by another polynomial basis (e.g., Legendre, ..). It is also possible to represent the time series with a wavelet basis function decomposition, and to apply a PCA on the resulting coefficients. This method, presented in [START_REF] Jo | Feature extraction across individual time series observations with spikes using wavelet principal component analysis[END_REF], has the advantage to represent the time series in both frequency and time domain, and could be an interesting alternative.

It is also possible to avoid the functional basis decomposition approach. If the user is interested in grouping time series on the basis of a given similarity measure, the spectral dimension reduction is an interesting candidate. This approach consists in computing the Laplacian of the time series similarity matrix and use its eigenvectors as time series features. With this approach, any distance can be considered (e.g., DTW, Mahalanobis, . . . ). Other ways to make use of the dissimilarity matrix for time series representation include the Multi-Dimensional Scaling, which is similar to performing a PCA on the dissimilarity matrix. The results is a linear transformation in a smaller space, that conserves the similarity information (i.e., two elements close in the final space are also close in the original space). This transformation can also be performed with auto-encoder, or t-SNE. However, this set of approaches may be computationally expansive, as they require to compute np(np-1) 2 similarities before the block-clustering inference.

Inspired from [START_REF] Rohit | Using dynamic time warping distances as features for improved time series classification[END_REF], it is also possible to represent a time series with its dissimilarity to other reference elements. By using a small number of elements (for instance a ratio specified by the user), this method is equivalent to sub-sampling the dissimilarity matrix, and results in a lower dimensional feature basis representation. Again, a dimension reduction step may be useful if this feature space remains excessively large. This approach seems promising, but requires to correctly choose the reference time series.

Alternative Bayesian Non Parametric Prior

The DPMM framework used throughout this thesis is an handy solution to mixture model indeference and is one of the most popular BNP mixture model. However, the DPMM has several drawbacks :

-The DP prior tends to create partitions with few large clusters and several small ones, which is not always desirable.

-Even though the DPMM consistently estimates the mixture density [START_REF] Subhashis | Posterior convergence rates of Dirichlet mixtures at smooth densities[END_REF] and the mixing distributions [START_REF] Xuanlong | Convergence of latent mixing measures in finite and infinite mixture models[END_REF], it has been shown in [START_REF] Jeffrey | A simple example of Dirichlet process mixture inconsistency for the number of components[END_REF] that the method does not consistently converges to the true number of latent components.

The first limitation can be addressed with the introduction of a Pitman-Yor Process prior in place of the DP. This process uses a discount parameter to leverage the creation of equal-size clusters. However, it has also been shown in [START_REF] Jeffrey | Inconsistency of Pitman-Yor process mixtures for the number of components[END_REF] that this method still inconsistently estimates the number of components.

A recent proposal [START_REF] Jeffrey | An elementary derivation of the Chinese restaurant process from Sethuramans stick-breaking process[END_REF] introduces the Mixture of Finite Mixtures model, that is consistent for the number of components. This model is similar to the DPMM, but adds an explicit definition of a prior directly on the component number K. This prior distribution is only a consequence of the stick-breaking construction with the DP prior (cf. Eq. 1.6 in Sect. 1.3.2). The same paper presents the associated inference methods, that are similar to the existing DPMM inference methods.

Based on this thesis contributions, extending the Mixture of Finite Mixtures to a coclustering, multi-clustering and multi-coclustering framework seems straightforward, and would allow a better control of the number of components.

Conclusion

In the broad and wide domain of time series clustering, this thesis presented several contributions specifically designed for the exploration of time series datasets generated for ADAS validation.

The first contribution addresses univariate time series datasets, and assumes the presence of latent scenario. This method uses the mixture model framework for time series segmentation and driving pattern detection, and the hierarchical clustering heuristic for the final clustering of time series re-coded as symbol sequences.

The next three contributions are block clustering methods that make full use of the model-based coclustering and multi-clustering methods. These methods address the clustering of MTS based on a Fourier representation of the time series. The use of the BNP framework allow to integrate a model selection step while keeping the anomaly detection and prediction interval features.

These methods are experimented on simulated datasets, but also applied on real-life datasets generated for ADAS validation. These applications proved the contributions interest for industrial purposes. The methods source code is open and made available at github repositories for reproducibility.

This thesis, staged in convention between the LIPN at Université Sorbonne Paris Nord and Groupe Renault, opens several perspectives of extension, both for the dictionarybased and the block-clustering-based methods. 1 + κ 0 κ 0 + 1 (yµ 0 ) T Ψ -1 0 (yµ 0 )

Glossary

-(ν 0 +1)/2
, which is the density of a multivariate t-distribution t ν 0 -d+1 y|µ 0 , (κ 0 +1)Ψ 0 κ 0 (ν 0 -d+1) .

A.2 Posterior Predictive Distribution

With the same notation, and denoting x = (x i ) n a set of observations, the posterior distribution follows a NIW distribution [Ben+19 ; Mur07] :

p(µ, Σ | x, µ 0 , κ 0 , Σ 0 , ν 0 ) = N IW (µ, Σ | µ n , κ n , Σ n , ν n ),
with the updated hyper-parameter values obtained by : µ n = κ 0 µ 0 + nx κ 0 + n , κ n = κ 0 + n, ν n = ν 0 + n,

Ψ n = Ψ 0 + i x i x T i + κ 0 µ 0 µ T 0 -κ n µ n µ T n .
With

Ψ ′ = κ n µ n vec(κ n µ n µ T n + Ψ n ) and ν ′ = κ n ν n + d + 2
, the posterior predictive distribution can be written as : , and

p(y | Ψ ′ , ν ′ ) =
ν ′ + 1 = κ 1 ν 1 + d + 2 = κ n ν n + d + 2 + 1 1
. As in the prior predictive case, injecting the updated hyper-parameters values into Eq. A.3 gives the following expression :

p(y | x, Ψ, ν) = π -d/2 • κ d/2 n κ d/2 1 • Γ d (ν 1 /2) Γ d (ν n /2) • |Ψ n | νn/2 |Ψ 1 | ν 1 /2 ,
which is the density of the multivariate t-distribution t νn-d+1 y|µ n , (κn+1)Ψn κn(νn-d+1)

A.3 Joint Predictive Distribution

In the Bayesian Coclustering or Multi-Coclustering, the inference algorithm uses the (prior or posterior) predictive distribution of a row, i.e., a joint predictive distribution. This distribution also has a closed-form when the prior is conjugate. In the joint prior predictive case, the distribution of a vector of observations y = (y i ) m is : with T (y) = i T (y i ). After simplifying the integral (c.f. Eq. A.1), the posterior can be written as : The same derivations yields the following joint posterior predictive distribution expressions :

p(y | Ψ, ν) =
p(y | x, Ψ, ν) = i h(y i ) f (Ψ ′ , ν ′ ) f (Ψ ′ + T (y), ν ′ + m1) ,
that reduces to the following closed-form :

π -md 2 • κ d/2 n κ d/2 m • Γ d (ν m /2) Γ d (ν n /2) • |Ψ n | νn/2 |Ψ m | νm/2 .
The prior and posterior predictive distributions closed forms are summarized in Tab. A.1.

Clustering Multi-Blocs et Visualisation Analytique de Données Séquentielles Massives Issues de Simulation du Véhicule Autonome

Le développement de systèmes d'aide à la conduite demeure un défi technique pour les constructeurs automobiles. La validation de ces systèmes nécessite de les éprouver dans un nombre considérable de contextes de conduites. Pour ce faire, le Groupe Renault a recourt à la simulation massive, qui permet de reproduire précisément la complexité des conditions physiques de conduite et produit une grande quantité de séries temporelles multivariées. Le rôle de l'expert métier est alors d'explorer ces données pour déterminer précisément les capacités du système d'aide à la conduite étudié et, si nécessaire, de raffiner sa conception.

Un jeu de données simulées peut contenir jusqu'à plusieurs centaines de milliers de simulations, décrites par plusieurs centaines de variables. Face à de telles dimensions, le travail du développeur requiert une grande expertise, ainsi qu'un travail d'investigation minutieux et chronophage.

L'objectif de cette thèse est de produire des algorithmes et outils d'aide à l'exploration, la structuration et l'analyse de ces jeux de données issues de simulation massive. Basées sur des hypothèses de construction argumentées, nous proposons trois approches probabilistes de classification non supervisée de séries temporelles, adaptées à des jeux de données temporelles univariées et multivariées.

La première contribution de cette thèse est construite sur une hypothèse de construction par scénario, et réalise une classification de séries temporelles univariées par dictionnaire. Partant d'un modèle existant de segmentation polynomiale par morceaux, cette méthode crée un dictionnaire de régimes de conduite pour recoder les séries temporelles en phases similaires. Les contributions suivantes se concentrent sur la classification de jeux de données multivariés par bloc, qui regroupent les variables en se basant sur leur distribution et leur partitionnement.

La suite des contributions se concentre sur l'analyse de jeux de données multivariés, grâce à des modélisations multi-blocs qui regroupent les différentes variables en fonction de leur distribution.

La pertinence de ces contributions est illustrée par des applications issues des cas d'usage industriels Renault : la validation d'un système de freinage d'urgence, d'un système d'aide au maintien dans la voie, et d'un système d'évitement d'urgence.

  Systèmes d'aide à la conduite (ADAS) . . . . . . . . . . . . . . . . . . . . . . .

Figure . 1 .:

 1 Figure .1. : Système de classification des ADAS (source : SAE).

Figure . 2 .

 2 Figure .2. : Représentation schématique du système de simulation.

  3.

Figure . 3 .:Figure . 4 .:

 34 Figure .3. : Illustration du contenu d'un jeu de données simulées. Celui-ci est représenté par une matrice, contenant une série temporelles dans chaque cellule.

Figure . 5 .:

 5 Figure .5. : Structure des modèles proposés : a) Le Coclustering infère une partition colonne et une partition ligne ; b) le Multi-Clustering infère une une partition colonne et une partition ligne dans chaque groupe de variable ; c) le Multi-Coclustering combine les deux aspects et ajoute une couche de structures de Coclustering dans le Multi-Clustering.

Figure . 6 .:

 6 Figure .6. : Chronologie des publications.

Figure . 8 .

 8 Figure .8. : Simulation System Representation.

Figure . 9 .:Figure . 10 .

 910 Figure .9. : Illustration of the simulation dataset contents. The dataset can be represented as a matrix, where each cell is a time series. Every time series in a given row share the same length.

Chapter 2 :

 2 Time Series Clustering with Model-Based Dictionary Representation.This chapter details the first contribution of this thesis, an univariate time series clustering method consisting of three steps : a) a model-based piecewise polynomial regression used for time series segmentation ; b) a model-based dictionary construction that regroups similar driving patterns ; c) a character chains clustering hierarchical clustering.

Figure 1 . 1 .

 11 Figure 1.1. : DTW measure alignment between two rectangular time series (generated with [Gio09]).

  Figure 1.2. :The Symbolic Aggregate ApproXimation (SAX) method recodes the time series with a dictionary built with the Piecewise Aggregate Approximation (PAA) representation (here, 9 segments described with a four-symbol dictionary).

Figure 1 . 3 .:

 13 Figure 1.3. : On the left the original time series, on the right its decomposition in a threedimensional trigonometric basis.

Figure 1 . 4 .

 14 Figure 1.4. : Non-linear time series representation based on reconstruction error optimization with a convolutional autoencoder.

Figure 1 . 5 .

 15 Figure 1.5. : Deep Self Organizing Map [For+19] architecture. The clustering is based on the simultaneous optimization of the reconstruction and SOM loss functions.

Figure 1 .

 1 Figure 1.6. : Dendrogram representation produced with hierarchical clustering algorithm.

Figure 1 . 7 .

 17 Figure 1.7. : Voronoï diagram associated with a k-means solution.

Figure 1 . 8 .:

 18 Figure 1.8. : Illustration of a Gaussian Mixture Model assumption : there are several latent gaussian components that generate the dataset. The ellipses are the Gaussian components isoprobability contours.

  and F (θ k ) = p(x|z = k). With this generative model definition, sampling an observation x is performed by first drawing a component membership z ∼ M ult(π) then drawing from F (θ z ). The associated graphical model is represented in Fig. 1.9.

Figure 1 . 9 .

 19 Figure 1.9. : Mixture Model graphical model

Algorithm 2 :

 2 (-/S/C) Expectation Maximization Algorithm for Mixture Model Inference Input : cluster number K, iteration number M Output : MM parameters Θ, partition ẑ Initialize partition z (0) Estimate initial mixture parameters

Figure 1 .

 1 Figure 1.10. : a) Illustration of G with four components, when the component parameters are bivariate normal distributions modes, i.e, elements of R 2 . b) Stick representation of the component proportions.

Figure 1 .

 1 Figure 1.11. : Dirichlet Process Mixture Model graphical model

. 5 )Figure 1 .

 51 Figure 1.12. : a) Illustration of the first components of G, when the component parameters are bivariate normal distributions modes, i.e, elements of R 2 . b) Stick representation of the component proportions. This illustration illustrates the link between the Stick-Breaking construction and the fast proportion decrease.

  the density function associated with component k. It is interesting to note the presence of the same prior membership expressions n k n-1+α and α n-1+α than in the CRP definition.

Algorithm 3 :

 3 DPMM Inference input : α, G 0 , Iteration number M output : Estimated row-partition ẑ Initialize z (0) for m ← 1 to M do for i ← 1 to n do Remove x i from its cluster. Compute p(z i | z -i , X, α, G 0 ) as defined by eq. (1.7) and (1.8) Sample z (m+1) i

  9.

  Figure 1.13. : Comparisons of the multivariate clustering, coclustering, and multi-clustering structures.

Figure 1 .

 1 Figure 1.14. : Latent Block Model graphical model

Figure 1 .

 1 Figure 1.15. : Non-Parametric Latent Block Model graphical model

Figure 1 .

 1 Figure 1.16. : Multiple Partitions Model graphical model

Figure 1 .

 1 Figure 1.17. : Non-Parametric Multi-Clustering graphical model

Figure 2 .

 2 Figure 2.1. : SDLHC step 2 : From time-series to categorical sequences.

Figure 2 . 2 .:

 22 Figure 2.2. : Segmentation result sample.

Figure 2 . 3 .:

 23 Figure 2.3. : Segment sequence encoded using the dictionary. Two stationary patterns can be recognized (b and c), corresponding to cruise speed phases, as well as two accelerating (d and e) and one decelerating (a).

Figure 2 . 4 .:

 24 Figure 2.4. : Dictionary produced in the AEB use case.

Figure 2 . 5 .:

 25 Figure 2.5. : ARIs scores of various clustering approaches as a function of the number of clusters.

Figure 3 . 1 .:

 31 Figure 3.1. : The true multi-clustering partition and a block partition that can be produced with a coclustering method. The coclustering inadequation produces more row-clusters than necessary and several blocks are wrongly splitted.

Figure 3 . 2 .:

 32 Figure 3.2. : Differences between the LBM, MPM and CLBM. In CLMB c), every cells belonging to the same block follow the same block component distribution, as is the case in the LBM case a)

Figure 3 . 3 .:Figure 3 . 4 .

 3334 Figure 3.3. : Prototypes used as block mode for the simulations

Figure 3 . 5 . 3 . 3 Figure 3 .

 35333 Figure 3.5. : Distribution (median and quantiles 0.1,0.9) of Row, Cluster and Block partition ARI obtained after 30 SEM-Gibbs runs with different initialization methods.

#FrequencyFigure 3 . 7 .:

 37 Figure 3.7. : Results of 50 runs of each selection model strategy, in terms of differences to the generative model structure.

Figure 3 . 8 .:

 38 Figure 3.8. : Use case illustration : Ego drifts from its lane and cross the white line on the side of the road, before being put back in the runway center

Figure 3 . 9 .:

 39 Figure 3.9. : Final structure obtained on real case dataset

Figure 3 . 3 . 4

 334 Figure 3.10. : Ego lateral position in Block Cluster (5,1)

Figure 3 .Figure 3 . 3 . 4

 3334 Figure 3.11. : Ego lateral position in Block Cluster (2,1)

Figure 3 .Figure 3 .Figure 3 .

 333 Figure 3.13. : Control logic activation and changes in Ego's heading in Block Cluster (3,3)

Figure 4 . 1 .:

 41 Figure 4.1. : Scores versus Alpha and Iterations -with # Block in log scale

Figure 4 . 2 .:

 42 Figure 4.2. : Scores versus number of PCA axes.

Figure 4 . 3 .:

 43 Figure 4.3. : Scores versus logPeriodogram dimension.

  The time series are represented in a common log-periodogram with 40 coefficients and then projected on the 3 dominant PCA axes. The concentration parameters are both set to 10 -2 and the NIW parameters to the default values discussed in Sect. 4.2.4.

Figure 4 . 4 .:

 44 Figure 4.4. : Use case illustration : ego drifts from its lane, crosses the center line and heads toward an oncoming vehicle. The system detects the target and change ego's direction.

Figure 4 . 5 .:

 45 Figure 4.5. : Resulting coclustering on Emergency Lane Keeping (ELK) dataset. The result consists of 6 row-clusters and 13 column-clusters

Figure 4 .

 4 Figure 4.6. : Top-left : two highly negatively correlated direction change signals ; top-right : ego lateral position in the 3 biggest observation clusters ; bottom-left : 2 correlated speed variables ; bottom-right : 3 outlier driving patterns in the 3 smallest observation clusters.

Figure 5 . 1 .

 51 Figure 5.1. : Coclustering, Multi-Clustering and Multi-Coclustering, with k, l, h the row-cluster, correlated cluster and redundant cluster indices (respectively). Color and pattern designate block membership.

)

  With probability proportional to (a γ + H + 1), draw γ's new value from Gamma(a γ + H, b γlog(y)) (iii) With probability proportional to p(b γlog(y)), draw γ's new value from Gamma(a γ + H -1, b rlog(y)) (5.5)

Figure 5 . 2 .:

 52 Figure 5.2. : Illustration of the ground truth coclustering structure, the uninformative variables and the misleanding variables. Color and pattern indicates block memberships. The time series profiles displayed on the left represent a sample of the prototypes used for block generation.

Figure 5 . 3 .

 53 Figure 5.3. : Multiple Coclustering structure overlaid on the ELK dataset (one cell = one time series). The red rectangle outlines the trajectory variables depicted in Fig. 5.4.

Figure 5 . 4 .:

 54 Figure 5.4. : Car trajectories distributed in the clusters, corresponding to the content of the red rectangle outline in Fig. 5.3.

Figure 5 . 5 .:

 55 Figure 5.5. : Redundant variable partition projected with MDS

  ), with A = (a, d), B = (a, e), C = (b, e), D = (b, f ).

1. 3 .

 3 On the left the original time series, on the right its decomposition in a threedimensional trigonometric basis. . . . . . . . . . . . . . . . . . . . . . . . . . 1.4. Non-linear time series representation based on reconstruction error optimization with a convolutional autoencoder. . . . . . . . . . . . . . . . . . . . . 1.5. Deep Self Organizing Map [For+19] architecture. The clustering is based on the simultaneous optimization of the reconstruction and SOM loss functions. 1.6. Dendrogram representation produced with hierarchical clustering algorithm. 1.7. Voronoï diagram associated with a k-means solution. . . . . . . . . . . . . . 1.8. Illustration of a Gaussian Mixture Model assumption : there are several latent gaussian components that generate the dataset. The ellipses are the Gaussian components isoprobability contours. . . . . . . . . . . . . . . . . . 1.9. Mixture Model graphical model . . . . . . . . . . . . . . . . . . . . . . . . . 1.10. a) Illustration of G with four components, when the component parameters are bivariate normal distributions modes, i.e, elements of R 2 . b) Stick representation of the component proportions. . . . . . . . . . . . . . . . . . 1.11. Dirichlet Process Mixture Model graphical model . . . . . . . . . . . . . . . . 1.12. a) Illustration of the first components of G, when the component parameters are bivariate normal distributions modes, i.e, elements of R 2 . b) Stick representation of the component proportions. This illustration illustrates the link between the Stick-Breaking construction and the fast proportion decrease.

4. 4 .Lemma 2 . 2 . 1 . 2 π d/ 2 Γ( ν 0 -d+1 2 )

 4221222 Use case illustration : ego drifts from its lane, crosses the center line and heads toward an oncoming vehicle. The system detects the target and change ego's direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5. Resulting coclustering on Emergency Lane Keeping (ELK) dataset. The result consists of 6 row-clusters and 13 column-clusters . . . . . . . . . . . . . . . . 4.6. Top-left : two highly negatively correlated direction change signals ; topright : ego lateral position in the 3 biggest observation clusters ; bottom-left : 2 correlated speed variables ; bottom-right : 3 outlier driving patterns in the 3 smallest observation clusters. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1. Coclustering, Multi-Clustering and Multi-Coclustering, with k, l, h the rowcluster, correlated cluster and redundant cluster indices (respectively). Color and pattern designate block membership. . . . . . . . . . . . . . . . . . . . . 5.2. Illustration of the ground truth coclustering structure, the uninformative variables and the misleanding variables. Color and pattern indicates block memberships. The time series profiles displayed on the left represent a sample of the prototypes used for block generation. . . . . . . . . . . . . . . 5.3. Multiple Coclustering structure overlaid on the ELK dataset (one cell = one time series). The red rectangle outlines the trajectory variables depicted in Fig. 5.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4. Car trajectories distributed in the clusters, corresponding to the content of the red rectangle outline in Fig. 5.3. . . . . . . . . . . . . . . . . . . . . . . . 5.5. Redundant variable partition projected with MDS . . . . . . . . . . . . . .. The probability density function of the Normal Inverse Wishart distributionN IW (µ, Σ | ξ), with ξ = (µ 0 , κ 0 , Σ 0 , ν 0 ) can be written asp(µ, Σ | Ψ, ν) = f (Ψ, ν) exp η T Ψν T g(η) . µ 0 ) T Σ -1 (µµ 0 )κ 0 µ 0 vec(κ 0 µ 0 µ T 0 + Ψ 0 ) -κ 0 ν 0 + d + 2 T g(η)   = f (Ψ, ν) exp η T Ψν T g(η) with f (Ψ, ν) = κ κ 0 µ 0 vec(κ 0 µ 0 µ T 0 + Ψ 0 ) and ν = κ 0 ν 0 + d +TheoremThe prior predictive distribution of a random variable y drawn according to y ∼ N (µ, Σ), (µ, Σ) ∼ N IW (ξ), ξ = (µ 0 , κ 0 , Σ 0 , ν 0 ),is the multivariate t-distribution y | ξ ∼ t ν 0 -d+1 y|µ 0 , (κ 0 + 1)Ψ 0 κ 0 (ν 0d + 1).Démonstration. According to Lemmas 1 and 2, we find thatp(y | Ψ, ν) = µ,Σ p(y | µ, Σ)p(µ, Σ | Ψ, ν)d(µ, Σ) = η h(y) exp η T T (y) -1 T g(η) f (Ψ, ξ) exp η T Ψν T g(η) dη = η h(y) exp η T T (y) -1 T g(η) f (Ψ, ξ) exp η T Ψν T g(η) dη = h(y)f (Ψ, ν) η exp η T (Ψ + T (y)) -(ν + 1) T g(η) dη.and, using the matrix determinant lemma (i.e. |A + uv T | = (1 + v T A -1 u)|A|) gives |Ψ 1 | = 1 + κ 0 κ 0 + 1 (yµ T 0 )Ψ -1 0 (yµ T 0 ) |Ψ 0 |. Substituting these expressions into Eq. (A.2) gives the density p(y | Ψ, ν) = Γ( |Ψ 0 | 1/2

  | µ, Σ)p(µ, Σ | ψ ′ , ν ′ )d (µ, Σ) = η h(y) exp η T T (y) -1 T g(η) f (Ψ ′ , ν ′ ) exp η T Ψ ′ν ′T g(η) dη = h(y)f (Ψ ′ , ν ′ ) η exp η T (Ψ ′ + T (y)) -(ν ′ + 1) T g(η) dη.After integral simplification (c.f. Eq. A.1), the posterior can be written as :p(y | x, Ψ, ν) = h(y) f (Ψ ′ , ν ′ ) f (Ψ ′ + T (y), ν ′ + 1) . (A.3)These expressions are similar to the ones of the prior distribution, with updated hyperparameter values, withΨ ′ + T (y) = κ 1 µ 1 vec(κ 1 µ 1 µ T 1 + Ψ 1 ) = κ n µ n vec(κ n µ n µ T n + Ψ n ) + y vec(yy T )

1 • 2 |Ψ 1 | ν 1 / 2 1 • 2 |Ψ 1 | ν 1 ••

 12121211 , µ, Σ | Ψ, ν)d(µ, Σ) = µ,Σ p(y | µ, Σ)p(µ, Σ | ψ, ν)d (µ, Σ) = µ,Σ i p(y i | µ, Σ)p(µ, Σ | ψ, ν)d (µ, Σ) = η i h(y i ) exp η T T (y i ) -1g(η) f (Ψ, ξ) exp η T Ψν T g(η) dη = i h(y i )f (Ψ, ξ) η exp η T (Ψ + T (y)) -(ν + m1) T g(η) dη, Γ d (ν 1 /2) Γ d (ν 0 /2) • |Ψ 0 | ν 0 /posterior predictive p(y | x, Ψ, ν) π -d/2 • κ d/Γ d (ν 1 /2) Γ d (νn/2) • |Ψn| νn/Γ d (νm/2) Γ d (ν 0 /2) • |Ψ 0 | ν 0 /2|Ψm| νm/2 joint posterior predictive p(y | x, Ψ, ν) π Γ d (νm/2) Γ d (νn/2) • |Ψn| νn/2 |Ψm| νm/2

•

  p(y | Ψ, ν) = i h(y i ) f (Ψ, ν) f (Ψ + T (y), ν + m1) , (A.4) with Ψ + T (y) = κ m µ m vec(κ m µ m µ T m + Ψ m ) = κ 0 µ 0 vec(κ 0 µ 0 µ T 0 + Ψ 0 ) + i y i vec( i y i y T i )the same hyper-parameters update expressions than in the posterior predictive case detailed in the previous section Sect. A.2.Injecting the values of h, f and the prior hyper-parameters in Eq. A.4 gives the closedform of the joint prior predictive distribution :p(y | Ψ, ν) = π Γ d (ν m /2) Γ d (ν 0 /2) • |Ψ 0 | ν 0 /2 |Ψ m | νm/2 .

  3 Milliards de dollars en 2019, et prévoyait une hausse jusqu'à 142 Milliards de dollars en 2027.

			DRAFT-Stand alone
		SAE J3016 TM LEVELS OF DRIVING AUTOMATION TM
		You are driving whenever these driver support features	You are not driving when these automated driving
	What does the human in the	are engaged -even if your feet are off the pedals and you are not steering	features are engaged -even if you are seated in "the driver's seat"
	driver's seat have to do?	You must constantly supervise these support features; you must steer, brake or accelerate as needed to	When the feature requests,	These automated driving features will not require you to take
			maintain safety		you must drive	over driving
					These are automated driving features
		These features	These features	These features	These features can drive the vehicle	This feature
		are limited	provide	provide	under limited conditions and will	can drive the
		to providing	steering	steering	not operate unless all required	vehicle under
		warnings and	OR brake/	AND brake/	conditions are met	all conditions
		momentary	acceleration	acceleration
			support to	support to
			the driver	the driver
		• automatic	• lane centering	• lane centering	• traffic jam	• local driverless	• same as
	Example Features	emergency braking warning • lane departure • blind spot warning	OR • adaptive cruise control	AND • adaptive cruise control at the same time	chauffeur	taxi • pedals/ steering wheel may or installed may not be	level 4, but feature conditions can drive everywhere in all
	Avec pour objectif d'équipper des millions de voitures à travers le monde, la sécurité des
	utilisateurs (conducteur ou piéton) est la première préoccupation lors du développement

Table des matières

 des

Table des matières

 des

  -a.

	Méthode	Chapitre	Lien
	SDLHC	Chapitre 2 https ://tinyurl.com/sdlhc
	FunCLBM Chapitre 3 https ://tinyurl.com/FunCLBM
	FunNPLBM Chapitre 4 https ://tinyurl.com/funNPLBM
	FunMCC	Chapitre 5 À paraître prochainement
	Cette thèse propose des extensions du Coclustering, avec trois méthode probabilistes :

Table .1. : Liens vers les dépôts Github associés aux contributions

  , M., Azzag, H., & Giraldi, L. (2021). Multivariate Time Series Multi-Coclustering. Application to Advanced Driving Assistance System Validation. In ESANN 2021-29th European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning. ESANN.

	Plusieurs publications sont, au moment de la rédaction de cette thèse, toujours en cours
	de révision :
	-Goffinet, É., Lebbah, M., Azzag, H., Giraldi, L., & Coutant. A. (2021). Multiva-
	riate Time Series Multi-Coclustering. Application to Advanced Driving Assistance

Table des matières

 des System Validation. Soumis à Applied Intelligence : Special Issue on Multi-view Learning. Springer .

	-Goffinet, E., Lebbah, M., Azzag, H., Giraldi, L. & Coutant. A. (2021). Functional
	Non-Parametric Latent Block Model : a Multivariate Time Series Clustering Ap-
	proach for Autonomous Driving Validation. Soumis à Computational Statistics &
	Data Analysis.

Table des matières

 des

Table des matières

 des

Table . 2

 . . : Links to the Github repositories associated with the contributions

Chapter 3 : Multivariate Time Series And Functional Conditional Latent Block Mo- del

  This chapter introduces the first contribution dedicated to multivariate time series with multi-block clustering. This method consists in a parametric model-based Multi-Clustering model that extends an existing time series coclustering work by allowing to infer several row-partitions instead of one.

Chapter 4 : Bayesian Non-Parametric Coclustering

This chapter presents a Bayesian Non-Parametric multivariate time series coclustering. Based on an existing coclustering model, this method adds a symmetric Dirichlet Process prior on the block components. With this extension, this new model natively performs model selection.

Table des matières

 des Notations

	Table .1. : Notations summary
		Dataset
	s = (s i ) n	Univariate time series da-	S = (s i,j ) n×p Multivariate time series
	= (s i (t j ) T i ) n	taset		dataset
	n	Observation number	p	Variable number
	(T i ) n	Time series length
	x ∈ R n×d	Dataset after transforma-	X ∈ R n×p×d Dataset after transforma-
		tion of s		tion of S
	d	Representation space di-
		mension
	z	Observation memberships	K	Component number
	F	Distribution family	(θ k ) K	Components parameter
	π = (π k ) K	Mixture proportions	(µ, Σ)	Multivariate Gaussian pa-
				rameters
	M	Iteration number	Θ	Parameter set
	(n k ) K	Cluster Size
	α	Concentration parameter	K obs	Observed	component
				number
	G 0	Base Distribution	µ

X -i,. /X .,-j X without row i/column j d(., .) Dissimilarity measure x i,. /x .,j row i / column j of X Mixture Model and Dirichlet Process Mixture Model 0 , κ 0 , Ψ 0 , ν 0 NIW prior parameters χ Hyper-parameter set µ k , κ k , Ψ k , ν k Cluster k posterior distribution parameters Dictionary-based time series clustering s Time series

  Sample membership w j Update Θ given (z, w). model selection can be performed in several ways (see [Lom12 ; Ker+15 ; Bra14] reviews), usually with an exhaustive search associated to a model selection criterion.

	The

Algorithm 4 : SEM algorithm for Latent Block Model inference Input : row-cluster number K, column-cluster L, iteration number M Output : LBM parameters Θ, block partition (ẑ, ŵ) Initialize row-partition z Initialize column-partition w for m ← 1 to M do for i ← 1 to n do Estimate p(z i = k | w, x i,. ) given Θ Sample membership z i for j ← 1 to p do Estimate p(w j = l | z, x .,j ) given Θ

  Time series clustering is a large and open problem that has been addressed with an important number of approaches, based on dissimilarity measures, feature extraction methods, clustering algorithms, . . . The approach choice is made by the user, based on prior knowledge of the dataset, and leads to the creation of natural groups of observations, w.r.t. the clustering objective.

  K} ).with τ t,k = p(z t = k|s t , θ) the posterior segment membership distribution estimated in (2.3) during the expectation step, and f θ k ,t the density associated to cluster k regression model at time t. This optimization can therefore be achieved by the separate maximization of Q 1 and Q 2 . The optimization of Q 2 with respect to the parameters θ k = (β k , σ k ) consists in the estimation of a polynomial regression model on the points (s t ) T weighted by the posterior membership distribution (τ t,k )

T , which gives the following expressions :

βk = arg min

Table 2 .1. :

 2 Parameters grid for ARI evaluation.

	Method	Parameters	Range
	SAX	Number of segments	(5,10,20,30,40,50)
		Number of gaussian bins	(2,3,5,7,10,20,30,40,50)
	SDLHC	Dictionary size	(2,. . . ,12)
	MIXRHLP	Number of segments	(1,. . . ,10)
		Polynomial regression order	(1,. . . ,3)

Table 2 .

 2 2. : Adjusted Rand Index on the UCR archive datasets.

	Name	ℓ 2	DT W p DT W d KS SAX M ixRHLP SDLHC
	CBF	0.28	0.66	0.68	0.63 0.46	0.47	0.71
	OliveOil 0.46	0.53	0.40	0.50 0.00	0.40	0.55
	Trace	0.32	0.40	0.66	0.57 0.32	0.41	0.94

Table 3 .1. :

 3 Every model combinations for H = 2 and K m = 3, where each row-cluster in each column component.

	H m	K m	2	3	4	5	6	7	8	9	10
	2		5	9	14	20	27	35	44	54	65
	3		9	19	34	55	83	119	164	219	285
	4		14 34	69	125	209	329	494	714	1000
	5		20 55	125	251	461	791	1286	2001	3002
	6		27 83	209	461	923	1715	3002	5004	8007
	7		35 119 329	791 1715 3431	6434 11439 19447
	8		44 164 494 1286 3002 6434 12869 24309 43757
	9		54 219 714 2001 5004 11439 24309 48619 92377
	10	65 285 1000 3002 8007 19447 43757 92377 184755

Table 3 . 2

 32

. : Number of possible models with respect to K m and H m

Table 3 .

 3

	.3.

  4.2.1 describes the model FunNPLBM and the proposed inference algorithm. Benchmark and experiments are studied in Sect. 4.3 and, finally, a real-case application on an industrial dataset is presented in Sect. 4.4. The next sections, Sect. 4.2.1 and Sect. 4.2.2 give the formal model definition and the associated stochastic inference process. Finally, Sect. 4.2.3 and Sect. 4.2.4 dive deeper into the algorithm implementation details : multivariate Gaussian observation model, initialization strategy, final inference, hyperparameter specification, and complexity.

	4.2 Functional Bayesian Non-Parametric Latent Block
	Model
	This section introduces our proposal of functional Bayesian non-parametric latent block
	model, FunNPLBM.

Table 4 .

 4

	1. : Baselines summary
	Method	Description	Grid	# runs
	B-GM M	Decoupled clustering	True cluster numbers	2
	B-GM M 14 -best among several runs	True cluster numbers	14
	B-GM M M S	-with model selection	{1, . . . , 7} × {1}, {1} × {1, . . . , 7}	14
	LBM	Coclustering	True block number	1
	LBM 49	-best among several runs	True block number	49
	LBM M S	-with model selection	{1, . . . , 7} × {1, . . . , 7}	49
	B-DP M	Decoupled DPM	-	2
	N P LBM		-	1

Table 4 .

 4

	2. : Performance comparison of different coclustering methods
	Method	ARI	RI	NMI	K runtime (s)
	B-GM M	0.825 0.982 0.946 25	18.6
	B-GM M M S 0.942 0.994 0.9803 30	85.9
	B-GM M 14 0.979 0.997 0.994 25	89.3
	LBM	0.823 0.913 0.887 25	17.5
	LBM M S	0.940 0.994 0.979 25	494.4
	LBM 49	1	1	1	25	603.2
	B-DP M M 0.670 0.958 0.906 16	25.3
	NPLBM	1	1	1	25	40

  The Symbolic Aggregate ApproXimation (SAX) method recodes the time series with a dictionary built with the Piecewise Aggregate Approximation (PAA) representation (here, 9 segments described with a four-symbol dictionary). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	1.2.
	ADAS	Advanced Driver Assistance	AEB	Autonomous Emergency Bra-
		System		king
	AIC	Akaike Information Criterion	ARI	Adjusted Rand Index
	BIC	Bayesian Information Crite-	BNP	Bayesian Non Parametric
		rion
	CLBM	Conditional LBM	CRP	Chinese Restaurant Process
	DP	Dirichlet Process	DPM	Dirichlet Process Mixture
	DTW	Dynamic Time Warping	DWT	Discrete WaveletTransform
	Ego	Vehicle under test	ELK	Emergency Lane Keeping
	-/C/S EM -/Classification/Stochastic Ex-	fPCA	Functional Principal Com-
		pectation Maximization		ponent Analysis
	GMM	Gaussian Mixture Model	HMM	Hidden Markov Model
	ICL	Integrated Completed Likeli-	LBM	latent Block Model
		hood
	LKA	Lane Keeping Assist	MCC	Multi-Coclustering
	MCMC	Markov Chain Monte-Carlo	MPM	Multiple Partitions Model
	MTS	Multivariate Time Series	NIW	Normal Inverse Wishart
	NMI	Normalized Mutual Informa-	NN	Neural Network
		tion
	NPLBM	Non-Parametric LBM	PAA	Piecewise Aggregate Approxi-
				mation
	PCA	Principal Component Analysis PRM	Polynomial Regression Model
	PYP	Pitman-Yor Process	SAE	Society of Automotive Engi-
				neers
	SB	Stick-Breaking	SDLHC Segmentation-Dictionary-
				Levenshtein-Hierarchical-
				Clustering
	SOM	Self Organizing Map	SVD	Singular Value Decomposition
	SAX	Symbolic Aggregate approXi-	RHLP	Regression with Hidden Logis-
		mation		tic Process
	RI	Rand Index	SNE	Stochastic Neighbor Embed-
				ding
	UMAP	Uniform Manifold Approxima-	WLD	Weighted Levenshtein Dis-
		tion and Projection		tance

  1.13. Comparisons of the multivariate clustering, coclustering, and multi-clustering structures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.14. Latent Block Model graphical model . . . . . . . . . . . . . . . . . . . . . . . 1.15. Non-Parametric Latent Block Model graphical model . . . . . . . . . . . . . . 1.16. Multiple Partitions Model graphical model . . . . . . . . . . . . . . . . . . . 1.17. Non-Parametric Multi-Clustering graphical model . . . . . . . . . . . . . . . 2.1. SDLHC step 2 : From time-series to categorical sequences. . . . . . . . . . 2.2. Segmentation result sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3. Segment sequence encoded using the dictionary. Two stationary patterns can be recognized (b and c), corresponding to cruise speed phases, as well as two accelerating (d and e) and one decelerating (a). . . . . . . . . . . . . 2.4. Dictionary produced in the AEB use case. . . . . . . . . . . . . . . . . . . . . 88 2.5. ARIs scores of various clustering approaches as a function of the number of clusters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 3.1. The true multi-clustering partition and a block partition that can be produced with a coclustering method. The coclustering inadequation produces more row-clusters than necessary and several blocks are wrongly splitted. . . 92 3.2. Differences between the LBM, MPM and CLBM. In CLMB c), every cells belonging to the same block follow the same block component distribution, as is the case in the LBM case a) . . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.3. Prototypes used as block mode for the simulations . . . . . . . . . . . . . . . 3.4. ARI versus Log-Likelihood in 100 launches of SEM-Gibbs on the simulated dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5. Distribution (median and quantiles 0.1,0.9) of Row, Cluster and Block partition ARI obtained after 30 SEM-Gibbs runs with different initialization methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6. Best ARI obtained among several SEM-Gibbs runs (median, quantile 0.9 ; with variable number of concurrents) . . . . . . . . . . . . . . . . . . . . . . 3.7. Results of 50 runs of each selection model strategy, in terms of differences to the generative model structure. . . . . . . . . . . . . . . . . . . . . . . . . 3.8. Use case illustration : Ego drifts from its lane and cross the white line on the side of the road, before being put back in the runway center . . . . . . . . . 3.9. Final structure obtained on real case dataset . . . . . . . . . . . . . . . . . . 3.10. Ego lateral position in Block Cluster (5,1) . . . . . . . . . . . . . . . . . . . . 3.11. Ego lateral position in Block Cluster (2,1) . . . . . . . . . . . . . . . . . . . . 3.12. Ego lateral position in Block Cluster (3,1) . . . . . . . . . . . . . . . . . . . . 3.13. Control logic activation and changes in Ego's heading in Block Cluster (3,3) . 3.14. Uninformative signals in Block Cluster (1,1) : linearly increasing feature (vehicle's width, length, headlights activation..) . . . . . . . . . . . . . . . . 3.15. Uninformative signals in Block Cluster (1,1) : linearly increasing feature (distance to origin) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. Scores versus Alpha and Iterations -with # Block in log scale . . . . . . . . . 4.2. Scores versus number of PCA axes. . . . . . . . . . . . . . . . . . . . . . . . 4.3. Scores versus logPeriodogram dimension. . . . . . . . . . . . . . . . . . . . .

	Table des figures

Table A .

 A 1. : Summary of predictive distributions closed forms.

Table des matières

Chapitre 1

1.1 Time Series Clustering

1.2 Mixture Model Framework

E step estimating the Expectation of the memberships posterior probabilities p(z | x)given a fixed value of Θ ;

1.3 Bayesian Non-Parametric Modeling

1.4 Multivariate Time Series Model-Based Clustering

1.5 Multivariate Clustering With Block Structure

1.6 Conclusion

2.2 A three-step time-series clustering algorithm (SDLHC)

2.3 Experiments

3.2 Functional Conditional Latent Block Models

3.3 Experiments on Synthetic Data

3.4 Application to autonomous driving system validation

3.5 Conclusion

4.2 Functional Bayesian Non-Parametric Latent Block Model

4.3 Experimental setup

4.4 ADAS validation

Chapitre 4

4.5 Conclusion

5.2 Multiple Coclustering of multivariate time series

Chapitre 5

5.3 Experiments

5.4 Applications

6.1 Dictionary-based time series clustering

6.2 Multi-Coclustering

6.3 Conclusion

Bibliographie

Acknowledgement

Joint Prior Predictive Distribution

The joint prior predictive distribution, used in Eq. (4.2) to compute p(x i,. | w (m) , G 0 ) = l p(x i,l | G 0 ), has the following expression :

where n i,l is the size of x i,l and where the updated hyper-parameters values are obtained with :

Functional Bayesian Non-Parametric Latent Block Model

Joint Posterior Predictive Distribution Given X k,l , the observations in block (k, l), the block parameters posterior distribution is formally defined by

) the block posterior distribution parameters, updated from the prior :

With these parameters, the joint posterior predictive distribution needed in eq. ( 4.1) has the following expressions :

Implementation

In this section we give more details on the inference implementation : hyper-parameters specification, initializatin strategy, inference complexity G 0 hyper-parameters specification.

The clustered objects are PCA coefficients, which are centered. Therefore we set µ 0 to be the d-dimensional zero vector. The precision matrix Ψ 0 specification is a bit trickier and depends on assumptions on the dataset. For non-parametric autoregressive models, [START_REF] Schuurman | A comparison of inversewishart prior specifications for covariance matrices in multilevel autoregressive models[END_REF] compares several prior specifications for Ψ 0 and concludes that the dataset precision obtained with maximum likelihood estimation is a good standard, which we keep in our application. κ 0 and ν 0 , which represent the user's confidence in µ 0 and Ψ 0 , are set to their lowest value, as we want them as uninformative as possible. We discuss the specification of the other hyperparameters in the experiments (c.f., Sect. 4.3).

4). for j ← 1 to p h do

Sample w h j based on the same process as (5.3) and (5.4).

Update concentration parameters γ, α, β based on equation (5.5) Aggregate the sampled partitions (c.f. Sect. 5.2.3 - §3) to obtain the final consensus partitions v, ŵ and Ẑ.

-A non-parametric Multi-clustering (MC), equivalent to FunMCC without the coclustering layer.

The performances comparison is based on two popular scores : the Adjusted Rand Index (ARI), and the Normalized Mutual Information (NMI). The ARI is a corrected-for-chance version of the Rand Index that evaluates the proportions of correctly grouped and separated elements in a partition. The NMI is an entropy-based criterion from the information theory literature estimating the quantity of knowledge a partition gives on another.

In addition, we also compare the estimated number of clusters, K, as a complementary information. Because the selected baselines are stochastic inference processes, we always consider the scores mean obtained over ten launches (or the median in the case of the cluster number). When several row-partitions are produced (i.e., Multi-Clustering methods), the closest row-partition is considered, with respect to the ARI.

Démonstration. First, observe that y T Σ -1 y = Tr(Σ -1 yy T ) = vec(Σ -1 ) T vec(yy T ). As a consequence, using the usual definition of the density function yields

where

and vec is the vectorization operator.

The integral part is obtained by considering the Normal Inverse Wishart density with hyper-parameter (Ψ + T (y), ν + 1) :

This last expression yields the following equality

,

where

, and ν +1 =

. Introducing these parameters into the last expression of p(y | ξ) gives the result of the lemma.

On the other hand, |Ψ 0 | ν 0 /2 /|Ψ 1 | ν 1 /2 reduces because

κ 0 µ 0 + y κ 0 + 1

) κ 2 0 µ 0 µ T 0 + κ 0 µ 0 y T + κ 0 yµ T 0 + yy T = Ψ 0 + κ 0 κ 0 + 1 (yµ 0 )(yµ 0 ) T = Ψ 0 1 + κ 0 κ 0 + 1 (yµ 0 ) T Ψ -1 0 (yµ 0 )

A.1 Prior Predictive Distribution