Keywords: Estimation en ligne, Géométrie Riemannienne, Gradient d'information, traitement des signaux et Online estimation, Riemannian geometry, Information gradient, Signal and image processing

Titre : Une approche basée sur la géométrie Riemannienne pour l'estimation en ligne des distributions elliptiques

Résumé : Durant les dix dernières années, les méthodes basées sur la géométrie Riemannienne et la géométrie de l'information ont eu un impact important sur le traitement des signaux et des images, la science des données, et l'intelligence artificielle. L'objectif de cette thèse est de proposer de nouveaux algorithmes, basés sur la géométrie Riemannienne et la géométrie de l'information, pour l'estimation en ligne des lois dites à contours elliptiques, et de leurs mélanges. En général, l'estimation en ligne est réalisée à travers la minimisation d'une divergence statistique, la divergence de Kullback-Leibler, grâce à l'application d'une méthode de gradient Riemannien stochastique. Afin d'implémenter cette méthode, l'espace des paramètres de la famille de lois elliptiques (ou de mélanges de lois elliptiques) doit être équipé d'une métrique Riemannienne, de préférence la métrique d'information de Fisher. Malheureusement, pour les lois elliptiques, cette métrique est souvent inconnue, ou n'a pas d'expression analytique exploitable. Pour répondre à cette difficulté, nous avons introduit une alternative à la métrique d'information de Fisher, que nous avons appelée métrique d'information par composantes. En utilisant cette métrique, nous avons développé la méthode du gradient d'information par composantes. La méthode du gradient d'information par composantes est une méthode en ligne, avec un faible coût calculatoire, qui lui permet de prendre en compte les jeux de donnés massifs ou de grandes dimensions. De plus, cette méthode a deux variantes, une à pas d'optimisation décroissants, et l'autre à pas d'optimisation adaptatifs. Cette seconde variante permet d'éviter le choix manuel (habituellement très long et pénible) des pas d'optimisation, et d'atteindre une vitesse de convergence qui s'approche d'une vitesse exponentielle. Nous avons appliqué la méthode du gradient d'information par composantes à l'estimation de deux familles de lois elliptiques, les lois Gaussiennes généralisée multivariées, et les lois de Student multivariées, ce qui nous a permis de mettre en évidence à la fois son faible coût calculatoire et sa vitesse de convergence optimale. Finalement, nous avons réalisé des applications concrètes, en traitement des images et vision par ordinateur, à la conversion de couleurs et à la classification de textures. Pour les images de haute résolution (avec plus de 2 millions de pixels), notre méthode du gradient d'information par composantes n'a besoin que d'une centaine de secondes pour effectuer le travail, avec des résultats nettement meilleurs qu'avec les autres méthodes.

List of Figures

List of Tables component-wise information gradient of function with respect to θ η (n) step-size

List of notations

Bibliography

Résumé en Français

Cette thèse s'appuie sur trois articles de revue [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF][START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF][START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF], dont le premier est déjà paru dans le journal Entropy, et les deux autres sont en cours de révision, dans les journaux Signal Processing et IEEE Transactions on Signal Processing, respectivement. Sa contribution principale est de proposer une nouvelle famille de méthodes pour l'estimation des lois à contours elliptiques et de leurs mélanges, que nous avons appelées méthodes CIG (en anglais, component-wise information gradient).

Les lois à contours elliptiques

Rappelons qu'une loi à contours elliptiques est une loi de probabilité sur R m , qui a une densité de probabilité de la forme

p(x|θ) = c(β) det(Σ) g β (δ) (1)
où θ = (µ, Σ, β), avec µ ∈ R m le paramètre de position, Σ une matrice définie positive, dite matrice de dispersion, et β > 0 le paramètre de forme. Dans [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF], c(β) est une constante de normalisation, g β une certaine fonction, dite fonction génératrice de la loi, et δ = (x -µ) † Σ -1 (x -µ).

Par ailleurs, un mélange de lois à contours elliptiques est une loi de probabilité sur R m , avec une densité de probabilité de la forme

f (x|θ) = K k=1 w k p(x|θ k ) (2)
où K est le nombre de composantes de mélange, et chaque composante de mélange p(x|θ k ) est de la forme [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF]. Les paramètres w k sont les poids du mélange, et vérifient w k ∈ (0, 1) avec K k=1 w k = 1. Le paramètre "général" θ est donné par θ = (w k , θ k ; k = 1, . . . , K),

où θ k = (µ k , Σ k , β k ).
Parmi les lois à contours elliptiques, nous avons souvent considéré les cas particuliers des lois gaussiennes généralisées multivariées, et des lois t de Student multivariées. Ces deux sous-familles de lois ont prouvé leur utilité en traitement des images, vision par ordinateur, traitement radar, et imagerie biomédicale [START_REF] Gupta | Generalized Gaussian scale mixtures: A model for wavelet coefficients of natural images[END_REF][START_REF] Verdoolaege | Multiscale colour texture retrieval using the geodesic distance between multivariate generalized gaussian models[END_REF][START_REF] Elguebaly | Bayesian learning of generalized Gaussian mixture models on biomedical images[END_REF][START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF].

Également, les mélanges de lois gaussiennes généralisées multivariées, ou de lois T de Student multivariées, sont reconnues en tant que généralisations des mélanges de lois gaussiennes multivariées, aux données à queues lourdes ou corrompues par des valeurs aberrantes [START_REF] Mario | Unsupervised learning of finite mixture models[END_REF][START_REF] Holzmann | Identifiability of finite mixtures of elliptical distributions[END_REF]. Ces mélanges ont été appliqués aux traitements des images et des vidéos, et plus généralement en analyse des données [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF][START_REF] Tan | Multivariate statistical models for image denoising in the wavelet domain[END_REF][START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF].

Contexte de la thèse

Le problème de l'estimation des lois à contours elliptiques est un problème classique, pour lequel il existe déjà plusieurs méthodes de référence. En revanche, l'estimation des mélanges de lois à contours elliptiques est un problème plus récent dans la littérature, et bien plus difficile.

En règle générale, les méthodes existantes sont basées sur l'estimation par maximum de vraisemblance (la méthode des moments, applicable à l'estimation des lois gaussiennes généralisées multivariées, est une exception à cette règle [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF]). Bien sûr, l'estimation par maximum de vraisemblance garantit la consistance (convergence des estimateurs vers les vrais paramètres) et l'efficacité asymptotique (la performance asymptotique optimale). Cependant, les estimateurs du maximum de vraisemblance, pour les lois à contours elliptiques et leurs mélanges, ne peuvent pas être calculés sous forme explicite. Plutôt, des méthodes itératives spéciales doivent être développées pour les obtenir.

Pour les lois à contours elliptiques, la méthode itérative la plus connue, pour calculer l'estimateur du maximum de vraisemblance, est la méthode du point fixe [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF][START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF]. Dans [START_REF] Sra | Conic geometric optimization on the manifold of positive definite matrices[END_REF], cette méthode est comparée à des méthodes d'optimisation riemannienne, comme le gradient conjugué ou BFGS riemannien.

Pour les mélanges de lois à contours elliptiques, l'approche standard pour calculer l'estimateur du maximum de vraisemblance reste la méthode espérance-maximisation (e.m.) [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Hasannasab | Alternatives to the EM algorithm for ML estimation of location, scatter matrix, and degree of freedom of the Student t distribution[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF]. Des méthodes plus sophistiquées n'ont été proposées que très récemment, comme la méthode d'optimisation proximale alternée de [START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF].

Dans le cadre de la thèse, nous avons commencé par recenser les méthodes existantes, pour réaliser l'estimation par maximum de vraisemblance des lois à contours elliptiques et de leurs mélanges. Nous avons classé ces méthodes en méthodes hors ligne/en ligne, et aussi en méthodes euclidiennes/riemanniennes.

Les méthodes hors ligne utilisent toutes les données disponibles, à chacune de leurs itérations. En revanche, les méthodes en ligne n'utilisent qu'une seule donnée, ou minibatch, à chaque itération. Les méthodes hors ligne sont mieux adaptées pour les jeux de données de taille modérée, alors que les méthodes en ligne sont mieux adaptées aux jeux de données de haute dimension ou de grande taille.

Par ailleurs, la différence entre les méthodes euclidiennes et riemanniennes est la suivante. Les paramètres des lois à contours elliptiques et de leurs mélanges obéissent à certaines contraintes non linéaires (par exemple, les matrices de dispersion doivent être définies positives). Les méthodes euclidiennes imposent ces contraintes à travers des techniques ad hoc (par exemple, rajouter une constante positive à la matrice de dispersion, afin qu'elle reste définie positive [START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF]), alors que les méthodes riemanniennes font appel à la géométrie intrinsèque des paramètres.

En réalité, les lois à contours elliptiques sont mieux adaptées à la modélisation des jeux de données de taille modérée, alors que les mélanges de lois à contours elliptiques répondent mieux aux jeux de données de grande taille (puisqu'on peut augmenter le nombre de composantes de mélange). Il n'est donc pas surprenant que les méthodes existantes pour l'estimation des lois à contours elliptiques soient toutes des méthodes hors ligne. Pour l'estimation des mélanges de lois à contours elliptiques, il est beaucoup plus intéressant de développer des méthodes en ligne. Or, on remarque que de telles méthodes n'ont été proposées que très récemment [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF][START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF]. De plus, il n'existe aucune méthode en ligne riemannienne, dans la littérature à présent.

C'est dans ce contexte que nous avons introduit, à travers cette thèse, des méthodes d'estimation en ligne riemanniennes, pour les lois à contours elliptiques et leurs mélanges.

Contribution de la thèse

La contribution principale de cette thèse est d'avoir introduit une nouvelle famille de méthodes, appelées méthodes CIG (pour component-wise information gradient), pour l'estimation des lois à contours elliptiques et de leurs mélanges. Le cas des lois à contours elliptiques a été traité en premier, dans [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF], et celui des mélanges un peu plus tard, dans [START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF].

Les méthodes CIG sont des méthodes d'estimation riemanniennes, basées sur une nouvelle géométrie pour les lois à contours elliptiques, également développée dans le cadre de la thèse. Cette géométrie est fondée sur une métrique riemannienne, que nous avons proposée, sous le nom de CIM (pour component-wise information metric).

La motivation pour l'introduction de la CIM, et ensuite des méthodes CIG, se trouve dans notre premier article de revue [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF]. Cet article a étudié une méthode d'estimation en ligne riemannienne, appelée méthode du gradient d'information riemannien.

La méthode du gradient d'information riemannien est une extension au cadre riemannien d'une méthode classique assez connue, qui est la méthode du gradient naturel [START_REF] Amari | Information geometry and its applications[END_REF]. C'est une méthode en ligne riemannienne, facile à utiliser, et avec d'excellentes propriétés de convergence.

En effet, la méthode du gradient d'information riemannien ne demande aucun calibrage de la part de l'utilisateur (précisément, la sélection des pas d'optimisation se fait de manière automatique). En même temps, elle est rapidement convergente et asymptotiquement efficace.

Cela étant dit, cette méthode a un champ d'application limité, comme elle ne peut être appliquée que si la métrique d'information de Fisher est connue sous forme explicite. Par exemple, elle ne permet d'estimer les lois à contours elliptiques que dans le cas où le paramètre de forme est déjà connu (dans ce cas, la métrique d'information de Fisher a été donnée sous forme explicite dans [START_REF] Berkane | Geodesic estimation in elliptical distributions[END_REF]).

Dans [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF], notre objectif était de trouver une méthode capable de remplacer la méthode du gradient d'information riemannien, pour l'estimation des lois à contours elliptiques "complètes", avec un paramètre de forme inconnu. Au lieu de chercher l'expression de la métrique d'information de Fisher, qui devenait de plus en plus compliquée, nous avons choisi de la remplacer par une autre métrique, plus simple et plus facile à calculer. C'est ainsi que nous avons introduit la CIM, qui est une approximation diagonale par blocs de la métrique d'information de Fisher.

Intuitivement, les méthodes CIG sont liées à la CIM de la même manière que la méthode du gradient d'information riemannien est liée à la métrique d'information de Fisher. Cela permet de retenir au moins une partie des propriétés de convergence de la méthode du gradient d'information riemannien, tout en élargissant son champ d'application aux lois à contours elliptiques avec un paramètre de forme inconnu.

La méthode CIG pour l'estimation des lois à contours elliptiques a deux versions, une version hors ligne et une version en ligne. La version en ligne est particulièrement intéressante, dans les applications aux jeux de données réels de haute dimension ou de grande taille. Par exemple, elle est capable de réaliser une transformation de couleurs, entre deux images qui font plus de 10 6 pixel, en 21 secondes, seulement. Pour cette même transformation, les méthodes existantes, dans l'état de l'art, mettent jusqu'à 2 heures, pour une performance comparable. En fait, les méthodes CIG en ligne ont même donné d'excellents résultats sur des images de haute définition, bien plus grandes.

Nous avons effectué des expériences numériques détaillées, pour comparer les méthodes CIG aux méthodes déjà existantes pour l'estimation des lois à contours elliptiques.

Ces expériences ont montré que la méthode CIG en ligne est similaire à la méthode des moments (de [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF]), par sa faible complexité calculatoire et son temps de calcul réduit, mais aussi similaire à la méthode du point fixe (de [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF]), grâce à sa performance supérieure. La méthode CIG en ligne réussit donc à réunir les meilleures qualités des méthodes existantes, à la fois en termes de temps de calcul et de performance.

Dans notre dernier article [START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF], nous avons étendu les méthodes CIG, de l'estimation des lois à contours elliptiques à l'estimation de leurs mélanges. Nous n'avons considéré que la méthode CIG en ligne, comme les lois de mélange sont habituellement utilisées pour modéliser des jeux de données plus compliqués, de grande taille. En revanche, nous avons encore développé deux versions de la méthode CIG en ligne, la première avec des pas d'optimisation décroissants, et la deuxième avec des pas d'optimisation adaptatifs. Nous avons appelé ces deux versions CIG-DS et CIG-AS (pour decreasing step-size et adaptive step-size, respectivement).

La méthode CIG-AS implémente une sélection adaptative des pas d'optimisation. De plus, on peut montrer qu'elle a un taux de convergence très rapide, puisqu'elle réalise ce qu'on appelle un taux de convergence linéaire (ce qui signifie que l'erreur d'optimisation décroît avec une vitesse exponentielle [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]). La méthode CIG-AS est très intéressante, d'un point de vue pratique aussi bien que théorique.

Cependant, dans les expériences avec des données simulées ou réelles, nous avons implémenté une méthode hybride, entre CIG-DS et CIG-AS. En effet, CIG-DS est préférable dans la phase préliminaire de l'optimisation, alors que CIG-AS est préférable dans la phase finale, à l'approche d'un point cible. De manière impressionnante, cette méthode hybride avait une meilleure performance que les méthodes existantes dans l'état de l'art, telles que la méthode du gradient stochastique euclidien et la méthode e.m. [START_REF] Hasannasab | Alternatives to the EM algorithm for ML estimation of location, scatter matrix, and degree of freedom of the Student t distribution[END_REF][START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF]. De plus, cette méthode hybride avait un temps de calcul deux à trois fois plus faible que les méthodes de l'état de l'art.

En résumé, les méthodes CIG ont plusieurs avantages sur les méthodes existantes pour l'estimation des lois à contours elliptiques et de leurs mélanges.

• les méthodes CIG en ligne n'utilisent qu'une seule donnée ou mini-batch (de taille constante), pour chaque itération. Cela leur permet de prendre en compte des jeux de données de haute dimension ou de grande taille, avec un temps de calcul assez réduit. A l'inverse, la majorité des méthodes existantes sont des méthodes hors ligne (point fixe ou e.m. [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF]), qui doivent utiliser des ressources en temps et mémoire beaucoup plus importantes, pour traiter des jeux de données plus grands.

• la méthode CIG-AS implémente une sélection adaptative, complètement automatique, des pas d'optimisation. Cela permet d'éviter une phase de calibrage avec sélection manuelle des pas, habituellement longue et pénible. La majorité des méthodes d'estimation en ligne sont très sensibles aux pas d'optimisation, et doivent faire appel à une phase de calibrage (voir la discussion dans [START_REF] Vaswani | Painless stochastic gradient: Interpolation, line-search, and convergence rates[END_REF]).

• les méthodes CIG en ligne garantissent une convergence rapide et une précision élevée des estimateurs, même en comparaison avec des méthodes hors ligne. Les autres méthodes en ligne décrites dans l'état de l'art (notamment [START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF][START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF]) restent difficiles à mettre en place lorsque les paramètres de forme des lois à contours elliptiques sont inconnus. Dans ce cas, ou bien elles manquent de précision, ou bien elles ont besoin d'utiliser des mini-batch d'une taille croissante.

Organisation de la thèse

La thèse contient sept chapitres, dont le contenu s'organise de la façon suivante.

• le chapitre 1 est un chapitre introductif. Il rappelle la définition des lois à contours elliptiques et de leurs mélanges, et présente de façon générale les différentes approches à leur estimation. Ensuite, il décrit la contribution et l'organisation de la thèse.

• le chapitre 2 contient une revue des travaux existants dans la littérature, sur la géométrie des lois à contours elliptiques, ainsi que sur l'estimation de ces lois et de leurs mélanges • le chapitre 4 introduit la méthode du gradient d'information riemannien, basée sur notre premier article en revue [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF]. Cette méthode est applicable au cas particulier des lois à contours elliptiques avec un paramètre de forme connu. Ses propriétés de convergence sont énoncées dans les Propositions 1 à 5. Celles-ci montrent que la méthode du gradient riemannien d'information converge rapidement, et qu'elle est asymptotiquement efficace, tout en ayant une sélection automatique des pas d'optimisation.

• le chapitre 5 introduit la méthode CIG, pour l'estimation des lois à contours elliptiques complètes, issue de notre article en revue [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF]. La section 5.2 décrit la méthode CIG hors ligne (Algorithme 1), et donne des conditions suffisantes pour sa convergence, dans la proposition 6. La section 5.3 décrit la méthode CIG en ligne (Algorithme 2), ainsi que ces propriétés théoriques (taux de convergence, normalité asymptotique), dans les propositions 7 à 9.

• le chapitre 6 introduit la méthode CIG en ligne pour l'estimation des mélanges de lois à contours elliptiques, issue de notre article en revue [START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF]. La section 6.

ECD and their mixtures

Elliptically-contoured distributions (ECD) are a far-reaching extension of multivariate Gaussian distributions. They go by this name because, when an ECD has a probability density function, its contours (level surfaces) are multi-dimensional ellipsoids.

The common centre of these ellipsoids is determined by the location parameter (or expectation) µ ∈ R m , while the squares of the principal axes are proportional to the eigenvalues of the scatter matrix Σ ∈ P m (here, P m denotes the space of m × m symmetric positivedefinite matrices). An additional shape parameter β ∈ R + (For some distributions it is also called degree of freedom or scale parameter) determines the factor for this proportionality.

Let X be an m-dimensional random vector following an ECD model. Its probability density function (if it exists) takes on the form

p(x|θ) = c(β) 1 det(Σ) g β (δ) (1.1)
where θ = (µ, Σ, β), the constant c(β) is a normalizing factor, and δ = (x -µ) † Σ -1 (x -µ).

Here, † denotes the transpose of the matrices. The density generator function g β depends on the specific sub-family of ECD distributions. ECD were originally introduced in [START_REF] Kelker | Distribution theory of spherical distributions and a location-scale parameter generalization[END_REF] and investigated in [START_REF] Cambanis | On the theory of elliptically contoured distributions[END_REF][START_REF] Fang | Generalized multivariate analysis[END_REF][START_REF] Frahm | Elliptical copulas: applicability and limitations[END_REF]. They inherit a lot of the nice properties of multivariate Gaussian distributions, and also make up for some of their deficiencies. Compared with multivariate Gaussian distribution, ECD are more flexible and powerful (e.g. they allow for heavy tails), and contain many widelyused sub-families of statistical distributions (elliptical Gamma, Pearson type II, elliptical multivariate logistic distributions, ...).

In terms of applications, the most popular sub-families of ECD are Multivariate Generalized Gaussian Distributions (MGGD), and Multivariate Student t-Distributions [START_REF] Gómez | A multivariate generalization of the power exponential family of distributions[END_REF][START_REF] Kotz | Multivariate t-distributions and their applications[END_REF][START_REF] Gómez Sánchez-Manzano | A matrix variate generalization of the power exponential family of distributions[END_REF]. These are location-scale distributions, and are further parameterised by a shape parameter, or a degrees of freedom parameter. MGGD are used in image processing, as models for wavelet and curvelet coefficients, and as models for three-channel colour vectors, in image denoising, context-based image retrieval, image thresholding, texture classification, and image quality assessment [START_REF] Bazi | Image thresholding based on the EM algorithm and the generalized gaussian distribution[END_REF][START_REF] Boubchir | Multivariate statistical modeling of images with the curvelet transform[END_REF][START_REF] Cho | Image denoising based on wavelet shrinkage using neighbor and level dependency[END_REF][START_REF] Gupta | Generalized Gaussian scale mixtures: A model for wavelet coefficients of natural images[END_REF][START_REF] Scharcanski | A wavelet-based approach for analyzing industrial stochastic textures with applications[END_REF][START_REF] Verdoolaege | Multiscale colour texture retrieval using the geodesic distance between multivariate generalized gaussian models[END_REF]. MGGD are also used in video coding and denoising, radar signal processing, and biomedical signal processing [START_REF] Mukund | Robust Gaussian and non-Gaussian matched subspace detection[END_REF][START_REF] Elguebaly | Bayesian learning of generalized Gaussian mixture models on biomedical images[END_REF]. Some applications of Student t-distributions are presented in [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF], involving image denoising. In radar imaging, the Student t-distribution is related to the so-called G 0 model within the family of spherically invariant random vectors (SIRVs), which is largely exploited in the context of SAR or PolSAR imaging, for tasks such as despeckling, classification, segmentation or detection [START_REF] Bombrun | Hierarchical segmentation of polarimetric SAR images using heterogeneous clutter models[END_REF][START_REF] Chen | A novel statistical texture feature for SAR building damage assessment in different polarization modes[END_REF][START_REF] Fernández-Michelli | Unsupervised polarimetric SAR image classification using G 0 p mixture model[END_REF][START_REF] Frontera-Pons | A class of robust estimates for detection in hyperspectral images using elliptical distributions background[END_REF].

Mixtures of ECD generalise mixtures of multivariate Gaussians (we will use the abbreviation MECD). An MECD is a convex combination of a finite number of ECD. Its density function is given by

f (x|θ) = K k=1 w k p(x|θ k ) (1.2)
Where K is the number of mixture components (or order of the mixture). The parameters w k are called mixture weights and must satisfy w k ∈ (0, 1) and K k=1 w k = 1. The function p(x|θ k ) is the k-th component of the mixture, with its parameters θ k = (µ k , Σ k , β k ), given as in equation (1.1). We will use θ to denote the ordered set of all the parameters of this mixture model θ

= (w k , µ k , Σ k , β k ; k = 1, • • • , K).
MECD are capable of modelling general, real-world probability distributions. For example, they have been recognized as generalizations of Gaussian mixture models, for data with heavy tails or outliers [START_REF] Mario | Unsupervised learning of finite mixture models[END_REF][START_REF] Holzmann | Identifiability of finite mixtures of elliptical distributions[END_REF]. Mixtures of MGGD have been applied to object or action recognition in videos [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF]. Moreover, mixtures of Laplace distributions have been applied to image denoising [START_REF] Tan | Multivariate statistical models for image denoising in the wavelet domain[END_REF]. Mixtures of Multivariate Student T-distributions have been used for robust modelling, clustering and classification of data with outliers [START_REF] Jeffrey | Model-based clustering, classification, and discriminant analysis via mixtures of multivariate t-distributions[END_REF][START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF][START_REF] Lin | Capturing patterns via parsimonious t mixture models[END_REF][START_REF] Peel | Robust mixture modelling using the t distribution[END_REF]. Based on MECD model, a new clustering algorithm for PolSAR images segmentation is introduced in [START_REF] Roizman | A new clustering algorithm for PolSAR images segmentation[END_REF] Because ECD and MECD have been successful in real-world signal and image processing applications, much attention has been devoted to developing effective methods for estimating their parameters. The estimation of ECD is an important problem that has received extensive attention. On the other hand, the estimation of MECD is currently an emerging problem, which is significantly harder than the estimation of ECD.

The estimation problem

There already exist several well-established methods for ECD estimation, while a few methods for MECD estimation have been recently proposed in the literature. Practically all of these methods are based on maximum-likelihood estimation (MLE), but one notable exception is the method of moments (MM), which can sometimes be used for ECD estimation.

The method of moments was introduced for MGGD estimation, in [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF]. As its name indicates, it is based on the idea of matching first, second, and fourth order moments of an unknown ECD with their empirical estimates.

The main advantage of this method is its low computational cost. However, it is much less accurate than maximum-likelihood estimation, and cannot be applied at all, to ECD which do not have higher-order moments (e.g., higher order moments are infinite), such as Student T-distributions [START_REF] Jeffrey C Fuhrer | Estimating the linearquadratic inventory model maximum likelihood versus generalized method of moments[END_REF][START_REF] Casella | Statistical inference[END_REF].

Maximum-likelihood estimation focuses on maximising the log-likelihood function1 . The maximum-likelihood estimate, based on a dataset X = {x 1 , • • • , x T }, drawn from an unknown ECD density p(x|θ * ), is a global maximiser of the log-likelihood function

L(θ; X ) = T t=1 p (θ; x t ) (1.3)
where p (θ; x) = log p(x|θ) for the ECD density p(x|θ) in (1.1), and the maximum is taken over all θ = (µ, Σ, β).

Similarly, if the dataset X = {x 1 , • • • , x T } is drawn from an unknown MECD density f (x|θ * ), then the maximum-likelihood estimate is a global maximiser of the log-likelihood function

L(θ; X ) = T t=1 f (θ; x t ) (1.4)
where f (θ; x) = log f (x|θ) for the MECD density f (x|θ) in (1.2), and the maximum is taken over all

θ = (w k , µ k , Σ k , β k ; k = 1, • • • , K) recall K is the number of mixture components, in (1.2).
The advantage of maximum-likelihood estimation is that the maximum-likelihood estimate is known to converge to the true parameter θ * or θ * , and to be asymptotically efficient (roughly speaking [START_REF] Ibragimov | Statistical estimation-asymptotic theory[END_REF][START_REF] Aad | Asymptotic statistics[END_REF], it enjoys the best possible asymptotic performance). On the other hand, it is not immediately clear how to compute maximum-likelihood estimates, and special methods must be developed for this purpose.

For ECD and MECD estimation, these methods can be divided into offline and online, and also into Euclidean and Riemannian methods. An iterative method for computing maximum-likelihood estimates uses available data to update (at each new iteration) an approximation θ (n) , and obtain a new and improved θ (n+1) .

The method is called offline or batch if it uses all of the available data for each update. It is called online or stochastic if it uses only one datapoint or one mini-batch of datapoints, for each update. Offline methods are suitable for dealing with small or moderate-sized datasets, while online methods are suitable for high-dimensional or large-scale datasets.

Furthermore, the distinction between Euclidean and Riemannian methods is the following. The parameters of ECD and MECD are subject to certain non-linear constraints (such as scatter matrices being positive-definite). Euclidean methods enforce these constraints using ad hoc techniques (such as adding a positive number to the scatter matrix, to make sure it is positive-definite [START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF]), while Riemannian methods appeal to the intrinsic Riemannian geometry of the parameters.

Usually, ECD are better at modelling moderate-sized datasets, while MECD are better with large-scale datasets. This is the reason why most currently existing methods for ECD estimation are offline methods, perhaps the most popular among them being the fixedpoint method (FP) [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF][START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF]. For now, existing methods for MECD estimation are mostly based on expectation-maximization (EM), which is an offline method [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Hasannasab | Alternatives to the EM algorithm for ML estimation of location, scatter matrix, and degree of freedom of the Student t distribution[END_REF]. Very recently, a few works have considered online methods [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF][START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF]. To our knowledge, there still do not exist any general Riemannian online methods for MECD estimation.

Within this context, the present thesis will introduce new Riemannian and online estimation methods both for ECD and MECD, designed to overcome the major difficulties which we have observed with currently existing methods.

Contribution of this thesis

The contribution made in the present thesis is based on our three journal papers [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF][START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF][START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF]. Above all, it consists in a new family of methods, for the estimation of ECD and MECD, called CIG methods. The abbreviation CIG stands for component-wise information gradient, and will be explained below.

The present section will discuss some of the main features of CIG methods. These features can be summarised as follows.

• the CIG methods include the CIG online estimation method. For each iteration, this only requires one datapoint or one mini-batch of datapoints. This makes it able to process high-dimensional or large-scale datasets, within short amounts of time.

• the CIG methods include the CIG adaptive step-size method. This implements a fully automatic selection of optimisation step-sizes (learning rates). This makes it easy to use, since it avoids the cumbersome task of manual selection of step-sizes.

• the CIG methods are able to guarantee fast rates of convergence, and high estimation accuracies, even in comparison to offline methods, which process the entire dataset at each iteration.

The importance of these features can be appreciated, in applications to real data. For example, in image editing, the CIG online method can be used to perform a colour transformation on a pair of images, with over 10 6 pixels, in only 21 seconds. For this same colour transformation, other state-of-the-art methods require two hours to do the same job, with a comparable performance. In fact, the CIG online method will be seen to produce excellent results, in applications with much larger, full HD images. Our motivation, for introducing the CIG methods, is the fact that the Fisher information metric (FIM) of ECD and MECD models can be very difficult to compute. This metric does have a tractable, closed-form expression, but only in the case of ECD models with a known and fixed shape parameter.

In this particular case, the closed-form expression of the FIM was given in [START_REF] Berkane | Geodesic estimation in elliptical distributions[END_REF]. In our journal paper [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF], we used it to introduce the Riemannian information gradient method. This is an online method, with automatic selection of step-sizes, which has, in a certain sense, the "best" convergence properties : fast rate of convergence and asymptotic efficiency.

The Riemannian information gradient method cannot be applied to complete ECD models (with unknown location and shape parameters), and cannot be applied to MECD models. This is because the FIM is not known, in a tractable closed form, for these models. In fact, for MGGD, there does exist an analytic expression of the FIM, but it is overly complicated [START_REF] Verdoolaege | On the geometry of multivariate generalized Gaussian models[END_REF].

In the journal paper [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF], our aim was to somehow extend the Riemannian information gradient method to complete ECD. Instead of pursuing increasingly complicated expressions of the FIM, we decided to replace it with another metric, which we called the component-wise information metric (CIM). Concretely, the CIM is a block-diagonal approximation of the FIM, and can be computed in a tractable closed form, even for complete ECD.

Roughly speaking, the CIG method for ECD estimation is based on the CIM, in the same way that the Riemannian information gradient method is based on the FIM. The CIG method for ECD estimation has two versions, CIG offline and CIG online.

The CIG method shows comparable, and sometimes superior performance, to that of state-of-the-art methods, such as the MM (method of moments) and FP (fixed-point) methods [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF]. However, the main advantage lies with its CIG online version, which requires significantly shorter amounts of time, in order to converge. Precisely, this CIG online version is similar to the more elementary MM method, for its short convergence time, and similar to the more sophisticated FP method, for its improved performance.

In the journal paper [START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF], we extended the CIG method, from the estimation of ECD models to the estimation of MECD models. We only considered the online version, because MECD are high-dimensional, and usually used to model more complex, large-scale datasets. However, we developed two sub-versions of the online version, one with decreasing step-sizes, and the other one with adaptive step-sizes, which we called the CIG-DS an CIG-AS methods, respectively.

The CIG-AS method implements an adaptive, fully automatic selection of step-sizes, and is theoretically shown to achieve a so-called linear rate of convergence (the term "linear convergence" means the optimisation error decreases exponentially fast [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]). This makes the CIG-AS method quite attractive, both from a practical and a theoretical perspective.

In computer experiments, with simulated or real data, we implemented a "hybrid method", combining CIG-DS and CIG-AS. The CIG-DS method is preferable in the early stages of optimisation, and the CIG-AS method in the later stages. Impressively, this hybrid method had better performance than state-of-the-art methods, such as Euclidean stochastic gradient, expectation-maximization with fixed-point method (EMFP) [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF], and proximal alternating linearized minimization (inertial stochastic PALM) [START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF]. At the same time, the computation time it requires in order to converge is two to three times less than any other existing method.

The reason why CIG-DS and CIG-AS outperform existing methods for MECD estimation is because they were designed to overcome the most problematic issues, observed with these methods. Most important among these issues were the following,

• offline methods (such as EMFP [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Hasannasab | Alternatives to the EM algorithm for ML estimation of location, scatter matrix, and degree of freedom of the Student t distribution[END_REF][START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF]), when applied to large-scale datasets, require very extensive computational resources (time and memory), making them quite impractical.

• existing online methods (described in [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF][START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF]) are difficult to implement when the shape parameters of ECD are unknown. In this case, they are either inefficient, or they require increasingly large mini-batch sizes.

• in addition, online methods involve the particularly delicate task of selecting stepsizes, which critically influence their rate of convergence [START_REF] Vaswani | Painless stochastic gradient: Interpolation, line-search, and convergence rates[END_REF].

The first of these issues is overcome thanks to the fact that CIG-DS and CIG-AS are online methods. The second issue is overcome thanks to the fact that CIG-DS and CIG-AS use the CIM to "pre-calibrate" stochastic gradients, which greatly improves their efficiency. Finally, the last issue is resolved thanks to the adaptive step-size selection, included in CIG-AS.

In addition to the above-mentioned contributions, based on the journal papers [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF][START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF][START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF], the CIG online method for ECD estimation was applied to change detection in multivariate image time series, in our additional journal paper [START_REF] Bouchard | Riemannian geometry for compound Gaussian distributions: Application to recursive change detection[END_REF].

Organization of the thesis

In addition to the current introductory chapter, the present thesis includes Chapters 2 to 7, as well as Appendices A to C. The contents of these chapters and appendices may be summarised as follows.

• Chapter 2 provides a detailed discussion of existing literature on the geometry of ECD, and on the estimation of ECD and MECD. In particular, Sections 2.3 and 2.4 discuss state-of-the-art methods for ECD and MECD estimation, respectively. These methods are finally summarised in Tables 2.1 and 2.2.

• Chapter 3 lays out our contribution to the geometry of ECD and MECD : the introduction of the component-wise information metric (CIM). Sections 3.3 and 3.4 define the CIM for ECD and MECD, respectively. Each one of these sections gives expressions of the CIM, and of the associated information gradient and retraction maps. These expressions are summarised in Table 3.5. For Chapters 5 and 6, it should be kept in mind that the information gradient associated to the CIM is called the component-wise information gradient (CIG).

• Chapter 4 introduces the Riemannian information gradient method, based on our journal paper [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF]. This method can be applied, in the particular case of ECD models with known shape parameter. Its main properties are stated in Propositions 1 to 5. These propositions show that the Riemannian information gradient method can achieve a fast rate of convergence, as well as asymptotic efficiency, with an automatic choice of step-sizes.

• Chapter 5 introduces the CIG method, for the estimation of complete ECD models (with unknown location and shape parameters). Section 5.2 presents the offline version of this method (Algorithm 1), and states its convergence in Proposition 6. Section 5.3 presents the online version of this method (Algorithm 2), and states its convergence, rate of convergence, and asymptotic normality, in Propositions 7 to 9. Section 5.4 studies the global convergence properties of the CIG method, for MGGD and Student t-distribution models, with known shape parameters. These properties are summarised in Tables 5.1 and 5.2.

• Chapter 6 extends the CIG method, from the estimation of ECD models to the estimation of MECD models. Section 6.2 introduces the decreasing step-size version of the CIG method for MECD estimation (Algorithm 3), while Section 6.3 is concerned with the adaptive step-size version (Algorithm 5). The convergence and rate of convergence of these two methods is given by propositions 11 and 12, respectively.

Introduction

This chapter gives a detailed presentation of the existing estimation methods for ECD and MECD models. Some of these methods rely on the Riemannian geometry of ECD models, which will therefore also be presented in this chapter.

The problem of estimating an ECD model is typically formulated as a mathematical optimisation problem, such as maximising the log-likelihood function (MLE estimation), or minimising the Kullback-Leibler divergence. This optimisation problem is then addressed using some kind of iterative method, such as the fixed-point method or gradient descent method.

To formulate these iterative methods (e.g. to define and compute the gradient), or to study their convergence (e.g. to determine whether a fixed-point iteration is contractive or not), it is first of all necessary to equip the parameter space with a geometry. The most classical choice is just that of a Euclidean geometry.

However, the Euclidean geometry has its limitations in applications. First of all, the parameters of the ECD model are subject to non-linear constraints (such as the scatter matrix being positive-definite). Using the Euclidean geometry means these constraints should be checked and then re-enforced at each iteration, leading to additional computations and time consumption, or even to numerical instability. Second, in Euclidean geometry, the convexity of the cost function (the negative log-likelihood or the Kullback-Leibler divergence) cannot be guaranteed, and global convergence of gradient descent is therefore difficult to obtain.

These issues are the main motivation for the introduction of the Riemannian geometry of ECD. Existing work on this geometry is presented in the following Section 2.2. The state of the art on estimation of ECD and MECD is then presented in Sections 2.3 (for ECD) and 2.4 (for MECD).

Geometry of ECD

Estimation methods for ECD and MECD can be classified into classical Euclidean methods and recent, more refined Riemannian methods. Riemannian methods take into account the intrinsic geometry of the parameter space of ECD and MECD models, and provide a better understanding of the structure of the optimisation problems involved in estimating these models.

For example, Riemannian methods automatically preserve fundamental non-linear constraints, such as positive-definiteness of the scatter matrix. Moreover, even when the negative log-likelihood or the Kullback-Leibler divergence is not convex, with respect to the scatter matrix, in the Euclidean sense, it is often convex in the Riemannian sense (geodesically convex) [START_REF] Wiesel | Geodesic convexity and covariance estimation[END_REF].

Thus, the introduction of a Riemannian metric leads to new Riemannian methods which have both practical and theoretical advantages. Mostly, existing works have focused on deriving a Riemannian metric for the scatter matrix, without considering the other parameters (location and shape parameters). These works consider two possibilities for the Riemannian metric of the scatter matrix.

The first, older metric is the affine-invariant metric introduce by [START_REF] Pennec | A Riemannian framework for tensor computing[END_REF]. The second is the Fisher information metric which was derived in [START_REF] Berkane | Geodesic estimation in elliptical distributions[END_REF]. In fact, the Fisher information metric is intimately related with the statistical properties of ECD models.

The Fisher information metric was first introduced by Rao [START_REF] Rao | Information and the accuracy attainable in the estimation of statistical parameters[END_REF]. It was later made popular by Amari [START_REF] Amari | Methods of Information Geometry[END_REF]. In this thesis, we will be interested in the Fisher information metric, because it allows dramatic improvements in the estimation of ECD and MECD. In the context of maximum-likelihood estimation, the gradient of the log-likelihood function, computed with respect to this metric, is the information gradient (sometimes also called the natural gradient). Estimation methods based on the information gradient are easy to implement and have excellent convergence properties, including fast rate of convergence and asymptotic efficiency. These properties would be very hard, or even impossible, to obtain without using the information gradient.

To begin, let us recall the definition and properties of the affine-invariant metric. The Fisher information metric will be discussed in Subsection 2.2.2.

Affine-invariant metric

Recall that an ECD has three parameters, namely the location parameter µ ∈ R m , the scatter matrix Σ ∈ P m , and the shape parameter β ∈ (0, ∞). Out of these three, due to the non-linear geometry of P m , the estimation of Σ is the most difficult. Therefore, existing contributions have mostly focused on the Riemannian geometry of P m .

Geometrically, P m is a cone, sitting inside the vector space of m × m symmetric matrices. An affine-invariant Riemannian metric on this cone was first introduced by [START_REF] Theil | A riemannian geometry of the multivariate normal model[END_REF]. For any point Σ ∈ P m and any tangent vectors U , V in the tangent space T Σ P m (these will always be real m × m symmetric matrices), this affine-invariant metric is given by

U, V ↑ Σ = tr(Σ -1 U Σ -1 V ) (2.1)
Here, instead of the conventional symbol " •, • ", we added a symbol "↑" in order to distinguish it from the Fisher information metric that will appear in the following section.

In the same way, the norm of this metric is denoted by

U ↑ Σ = U, U ↑ Σ .
The geodesics of the affine-invariant metric (2.1) are given by the exponential map [START_REF] Pennec | A Riemannian framework for tensor computing[END_REF]

Exp Σ (U ) = Σ 1 2 exp(Σ -1 2 U Σ -1 2 )Σ 1 2 = Σ exp(Σ -1 U ) (2.
2)

The geodesic is the generalization of the straight line, which is defined by a starting point Σ and a starting velocity U . For the starting velocity U , there exists a unique geodesic curve γ : [0, 1] → P m , t → γ(t) such that γ(0) = Σ, γ (0) = U and Exp Σ (U ) = γ(1). The introduction of the exponential map allows us to compute this geodesic. That is to say, the non linear constraints are respected by using this exponential map. Its inverse is the Riemannian logarithm

log Σ (U ) = Σ 1 2 log(Σ -1 2 U Σ -1 2 )Σ 1 2 = Σ log(Σ -1 U ) (2.3)
Another important concept is the Riemannian distance. For any pair of points Σ and T = Exp Σ (U ) in P m , we define the Riemannian distance d(Σ, T ) to be the infimum of the lengths of all admissible curves from Σ to T . The Riemannian distance corresponding to the metric (2.1) has the following expression,

d 2 ↑ (Σ, T ) = log Σ (T ) ↑ Σ 2 = tr log(Σ -1 T ) 2 Σ, T ∈ P m (2.4)
Here, all matrix functions, such as matrix power, exponetial and logarithm, are symmetric matrix functions [START_REF] Nicholas | Functions of matrices: theory and computation[END_REF]. Now consider the log-likelihood function p (θ; x) = log p(x|θ) where p(x|θ) is given by (1.1). Using the first order of Tyler's development, the Riemannian gradient of (Σ) is defined by the unique vector field that satisfies

∇ ↑ Σ p , U = d p (Σ)[U ],
where d is the differential operator. The Riemannian gradient of p (θ; x) with respect to Σ, computed using the affine-invariant metric (2.1) has the following expression

∇ ↑ Σ p (θ; x) = - 1 2 Σ - ∂h β (δ) ∂δ (x -µ)(x -µ) † (2.5)
where h β (δ) = log g β (δ) for the density generator function g β in Equation (1.1). The above Riemannian metric, exponential map, distance, and Riemannian gradient, can be used to define the geodesic convexity (g-convexity, for short) with respect to Σ. We say that a function f is geodesically convex, if for any geodesic curve γ(t), the composite function f • γ(t) is a convex function in the classical sense.

Definition 1 If Θ is a Riemannian manifold, a function f : Θ → R is said to be geodesically α-strongly convex if for any θ, θ ∈ Θ,

f (θ ) f (θ) + ∇ θ f (θ), Exp -1 θ (θ ) θ + α 2 d 2 (θ, θ )
where •, • θ denotes the Riemannian metric on Θ. The vector ∇ θ f (θ) means the Riemannian gradient of f , computed using this metric, Exp and d(•, •) mean the Riemannian exponential map and distance (again, with respect to the same Riemannian metric). For example, when θ = (Σ), the metric •, • θ can be the affine-invariant metric (2.1), with the gradient, exponential and distance given by (2.2), (2.4) and (2.5).

Specific to the ECD distribution family, the geodesic convexity of the log-likelihood function p (θ; x), with respect to Σ, has been studied by several authors [START_REF] Sra | Geometric optimisation on positive definite matrices for elliptically contoured distributions[END_REF][START_REF] Sra | Conic geometric optimization on the manifold of positive definite matrices[END_REF][START_REF] Zhang | Multivariate generalized Gaussian distribution: Convexity and graphical models[END_REF].

Fisher information metric

In addition to the affine-invariant metric (2.1), it is very interesting to consider the Fisher information metric of the ECD model, with respect to Σ, for fixed µ and β.

The Fisher information metric, as the affine-invariant metric, is defined by a scalar product U, V * = i,j I * (Σ)U i V j . The matrix [I * (Σ)] i,j is the Fisher information matrix of ECD model, which is a symmetric positive definite matrix and given by its definition

[I * (Σ)] i,j = -E ∂ 2 (Σ) ∂Σ i ∂Σ j
The analytical expression is given as [START_REF] Berkane | Geodesic estimation in elliptical distributions[END_REF] U,

V * Σ = I Σ,1 tr(Σ -1 U Σ -1 V ) + I Σ,2 tr(Σ -1 U )tr(Σ -1 V ) (2.6)
where * is the notation for Fisher information metric, not the conjugate transpose. The information coefficients I Σ,1 and I Σ,2 are defined by

I Σ,1 = 2A m(m + 2) I Σ,2 = A m(m + 2) - 1 4 A = E ∂h β (δ) ∂δ δ 2 (2.7)
with h β (δ) = log g β (δ) as in equation (2.5). This Fisher information metric shares the same exponential map and logarithmic map with the affine-invariant metric in (2.1). That is, both of these metrics have their exponential map given by (2.2) and logarithmic map given by (2.3). Similar to the affine-invariant metric, according to d * (Σ, T ) = Log Σ (T ) * , the Fisher information metric can also derive a Riemannien distance, that is analytically expressed as

d 2 * (Σ, T ) = I Σ,1 tr log(Σ -1 T ) 2 + I Σ,2 tr 2 log(Σ -1 T ) Σ, T ∈ P m (2.8)
From definition (2.6) of the Fisher information metric, it is possible to compute the information gradient with respect to Σ. This turns out to be

∇ * Σ p (θ; x) = J Σ,1 [∇ ↑ Σ p (θ; x)] ⊥ + J Σ,2 [∇ ↑ Σ p (θ; x)] (2.9)
Here, the vector ∇ ↑ Σ p (θ; x) denotes the affine-invariant Riemannian gradient (2.5). Moreover, the coefficients J Σ,1 and J Σ,2 are given by

J Σ,1 = 1 I Σ,1 J Σ,2 = 1 I Σ,1 + mI Σ,2
(2.10) in terms of I Σ,1 and I Σ,2 defined in (2.7), and the symbols ⊥ and denote the following orthogonal decomposition of ∇ ↑ Σ p (θ; x),

[∇ ↑ Σ p (θ; x)] = 1 m tr Σ -1 ∇ ↑ Σ p (θ; x) Σ (2.11a) [∇ ↑ Σ p (θ; x)] ⊥ = ∇ ↑ Σ p (θ; x) -[∇ ↑ Σ p (θ; x)] (2.11b)
The gradient (2.9) is our first example of an information gradient. The information gradient is a central theme in this thesis, and we attempt to use it directly, for gradient descent, whenever possible. This is because the information gradient is "pre-calibrated", thanks to the presence of coefficients such as J Σ,1 and J Σ,2 in (2.9), which greatly simplify the task of choosing step-sizes.

Estimation of ECD

Typically, single ECD models (that is, one ECD rather than a mixture) are not used for modeling large-scale or high-dimensional datasets. Therefore, existing approaches to the estimation of ECD have focused on offline methods, where there is a finite number of datapoints. For a given dataset

X = {x 1 , • • • , x T }, the objective function of MLE is L(θ; X ) = T t=1 p (θ; x t ) = T log c(β) - T 2 log det(Σ) + T t=1 h β (δ t ) (2.12)
where

δ t = (x t -µ) † Σ -1 (x t -µ) and h β (δ) = log g β (δ)
, in the notation of (1.1). The methods of maximizing this objective function can be classified into two categories, the first contains Euclidean (or classical) methods, and the second contains the methods based on Riemannian geometry.

Euclidean methods

In the context of Euclidean geometry, if the dataset is supposed to be centred, the MLE of (Σ, β) for MGGD was studied in [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF]. Also for MGGD, under the framework of block majorization-minimization, the MLE of all the three parameters (µ, Σ, β) is considered in [START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF]. For Student t-distributions, with the degree of freedom parameter being fixed, the MLE of (µ, Σ) is given in [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF], according to a new parameterization of the couple (µ, Σ).

Here is a more detailed description.

The estimation of Σ: In all the three works just mentioned ( [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF][START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF]), the fixed-point (FP) equation of Σ is obtained by differentiating the log-likelihood L(θ; X ) with respect to Σ. This follows from

∂L(θ; X ) ∂Σ = - 1 2 Σ -1 - 1 T T t=1 ∂h β (δ t ) ∂δ t Σ -1 (x t -µ)(x t -µ) † Σ -1
Setting this derivative to zero results in the FP (Fixed-Point method) equation

Σ = F (Σ) where F (Σ) = - 2 T T t=1 ∂h β (δ t ) ∂δ t (x t -µ)(x t -µ) † (2.13)
The sequence generated by the fixed-point iteration

Σ (n+1) = F (Σ (n)
) converges to a stationary point of L(Σ). Note that, in [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF], the scatter matrix Σ is normalized according to tr(Σ) = m and Σ = mM . Also in this same work [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF], for MGGD models, the existence and uniqueness of the estimator of Σ is proved, whenever β ∈ (0, 1).

The estimation of (µ, Σ): When both µ and Σ are unknown, the MLE can be achieved in two ways. The first one is implemented according to a new parameterization

S = Σ + µµ † µ µ † 1 therefore S ∈ P m+1 (2.14)
If the new random vector y is given by y † = (x † , 1), then the log-likelihood function can be reformulated as

L(S) = - T 2 log det(S) - T t=1 hβ (δ yt ) δ yt = y † t S -1 y t (2.15)
where the function hβ is proportional to the original h β . Therefore, the FP equation for S can be obtained by replacing Σ with S and then h β with hβ in (2.13). In [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF], for Student t-distributions, the maximization of this new function L(S) was proven to be equivalent to the maximization of L(µ, Σ). In addition, the existence and uniqueness of the MLE were also proved in [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF].

Another idea is to find the MLE of µ and Σ through an alternating iteration. This idea is realized, in the case of MGGD, in [START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF]. The derivative of L with respect to µ is (recall that, here, θ = (µ, Σ))

∂L(θ) ∂µ = -2 T t=1 ∂h β (δ t ) ∂δ t Σ -1 (x t -µ) its FP equation is derived by setting ∂L(θ) ∂µ = 0 µ = T t=1 ∂h β (δt) ∂δt x t T t=1 ∂h β (δt) ∂δt (2.16)
In [START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF], equations (2.13) and (2.16) are used in an alternating manner to iteratively update estimates of µ and Σ, leading to the MLE of (µ, Σ). Global convergence of this method is guaranteed under the block majorization-minimization framework (BMM framework) [START_REF] Meisam Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF].

The estimation of (µ, Σ, β): For the MGGD model, the estimation of the shape parameter β was achieved by solving a nonlinear equation [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF]. Let the derivative of L with respect to β equal to 0

∂L(θ) ∂β = T ∂ log c(β) ∂β + T t=1 ∂h β (δ t ) ∂β = 0 (2.17)
The update formula for β is then

β = solution of equation (2.17) (2.18)
The Newton-Raphson method is employed to solve equation (2.17) in [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF]. In [START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF], the FP equations (2.13), (2.16) and (2.18) are applied in an alternating manner to find MLE estimates of the complete ECD model. Note that the negative log-likelihood function -L(β; µ * , Σ * ) is proved to be strongly convex for any fixed µ = µ * and Σ = Σ * . Therefore, under the framework of BMM, this alternating procedure converges to a stationary point of the log-likelihood L(θ; X ) [START_REF] Wang | Globally convergent algorithms for learning multivariate generalized Ggaussian distributions[END_REF].

Riemannian methods

In the context of Riemannian geometry, all currently existing methods focus only on the MLE of the scatter matrix Σ, with the other parameters considered to be fixed [START_REF] Boukouvalas | A new Riemannian averaged fixed-point algorithm for MGGD parameter estimation[END_REF][START_REF] Sra | Geometric optimisation on positive definite matrices for elliptically contoured distributions[END_REF][START_REF] Sra | Conic geometric optimization on the manifold of positive definite matrices[END_REF][START_REF] Zhang | Multivariate generalized Gaussian distribution: Convexity and graphical models[END_REF]. The gradient of the MLE objective function (2.12), with respect to the affine-invariant metric (2.1), can be found from (2.5)

∇ ↑ Σ L(θ) = - T 2 Σ - T t=1 ∂h β (δ t ) ∂δ t (x t -µ)(x t -µ) † (2.19)
Setting this gradient to zero, the same equation as (2.13) can be obtained. This FP equation is used for MLE of MGGD in [START_REF] Zhang | Multivariate generalized Gaussian distribution: Convexity and graphical models[END_REF], and of general ECD in [START_REF] Sra | Geometric optimisation on positive definite matrices for elliptically contoured distributions[END_REF]. For MGGD, the global g-convexity of the log-likelihood function is also studied in [START_REF] Zhang | Multivariate generalized Gaussian distribution: Convexity and graphical models[END_REF]. When the shape parameter β is fixed and β ∈ [ 1 2 , ∞), the log-likelihood function is globally geodesically convex, in the sense of Definition 1. Note that this is a different range of β than the one found in [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF] using an Euclidean method.

A Riemannian-averaged fixed-point method was introduced in [START_REF] Boukouvalas | A new Riemannian averaged fixed-point algorithm for MGGD parameter estimation[END_REF], which overcame the problem of instability of the FP iteration for larger values of the shape parameter. Precisely, this method implements successive Riemannian averages of fixed-point iterates. Recall the definition of the Riemannian average of Σ, T ∈ P m . For t ∈ (0, 1), the Riemannian average with ratio t of Σ and T is

Σ# t T = Σ 1 2 (Σ -1 2 T Σ -1 2 ) t Σ 1 2
The Riemannian-averaged algorithm is defined as follows. When Σ (n) is given, instead of defining Σ (n+1) by F (Σ (n) ) according to (2.13), let

Σ (n+1) = Σ (n) # tn F (Σ (n) ), t n = 1 n + 1
Thanks to the stability of the Riemannian average, this method can give an efficient estimate of Σ for any value of the shape parameter β.

In addition to the fixed-point method, many Riemannian line-search methods are studied in [START_REF] Sra | Conic geometric optimization on the manifold of positive definite matrices[END_REF]. The line-search optimization methods on Riemannian manifolds are all based on the following general iterative formula

Σ (n+1) = Exp Σ (n) (η (n) U (Σ (n) , X )) (2.20)
where η (n) > 0 is a step-size and U (Σ (n) , X ) denotes a vector in T Σ (n) P m , which refers to the search direction (for example,

U (Σ (n) , X ) = grad Σ L(θ), defined in (2.19)).
With this general iterative formula, the most common line-search optimization methods, such as steepest descent, BFGS, and conjugate gradient methods, are compared with the FP methods, in [START_REF] Sra | Conic geometric optimization on the manifold of positive definite matrices[END_REF]. The conclusion is that the Riemannian line-search methods and the FP method each have their own advantages. Although the Riemannian linesearch methods have excellent performance, the FP algorithm is still competitive in many scenarios.

Estimation of MECD

Mixtures of ECD (MECD) are quite useful in fitting more complicated, large-scale datasets, of possibly higher dimension. For this reason, most existing methods for the estimation of MECD pay closer attention to online algorithms, which are time-saving and use up less memory. As in section 2.3, existing methods are classified into Euclidean and Riemannian.

Euclidean methods

In the Euclidean framework, the estimation of mixture models is classically addressed using Expectation-Maximization methods (EM) [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. In these methods, a latent variable,

(Z k t ; k = 1, . . . , K) is introduced to indicate the membership of each datapoint x t (then, Z k t = 1 if
x t belongs to the kth mixture component, and Z k t = 0 otherwise). The E (expectation) step computes the marginal log-likelihood of observed data. The M (maximization) step maximizes this marginal log-likelihood. These two steps are repeated until convergence.

To state this in a precise way, recall the MECD density (1.2) and its parameters

θ = (w k , µ k , Σ k , β k ; k = 1, • • • , K).
In the Expectation step, using the current estimates of the parameters

θ (n) = (w (n) k , µ (n) k , Σ (n) k , β (n)
k ; k = 1, • • • , K) the conditional expectation of the Z k
t is determined by Bayes theorem, which gives "membership probabilities"

o k (x t |θ (n) ) = w (n) k p(x t |θ (n) k ) K k=1 w (n) k p(x t |θ (n) k ) (2.21)
Then, the function Q(θ|θ (n) ) is defined to be the conditional expectation of the complete log-likelihood function

Q(θ|θ (n) ) = K k=1 T t=1 o k (x t |θ (n) ) log(w k ) + K k=1 T t=1 o k (x t |θ (n) ) log p(x t |θ k )
In the Maximization step, a Lagrange multiplier is added, for the weights w, in order to guarantee K k=1 w k = 1. After differentiating the function Q(θ|θ (n) ) with respect to w, the update formula for w is found to be

w (n+1) k = 1 T T t=1 o k (x t |θ (n) )
In [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF], the update formula for the other parameters (µ

(n) k , Σ (n) k , β (n)
k ), is performed by repeating the same operation as in equations (2.13), (2.16) and (2.18), for each subdistribution (i.e. each mixture component). The resulting method is called the EM with fixed point method (EM-FP). In the same work [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF], it was implemented for mixtures of MGGD, and then compared to the Euclidean stochastic gradient method.

For this Euclidean stochastic gradient method, the well-known logit parameterization is applied to ensure the constraints K k=1 w k = 1. This parameterization is given by

w k = exp(π k ) 1 + K-1 j=1 exp(π j ) , for k ∈ {1, • • • , K -1} and w K = 1 1 + K-1 j=1 exp(π j ) (2.22)
The logit parameterization yields the reparameterized log-likelihood function of the mixture model

f (θ; x) = log K-1 k=1 exp(π k ) 1 + K-1 j=1 exp(π j ) p(x|θ k ) + 1 1 + K-1 j=1 exp(π j ) p(x|θ K )
The Euclidean gradient with respect to the parameters π k is given by

grad π k f (θ; x) = o k - exp(π k ) 1 + K-1 j=1 exp(π j ) with o k = exp(π k ) 1 + K-1 j=1 exp(π j ) p(x|θ k ) f (x|θ) f (x|θ) = K-1 k=1 exp(π k ) 1 + K-1 j=1 exp(π j ) p(x|θ k ) + 1 1 + K-1 j=1 exp(π j ) p(x|θ K )
The stochastic gradient update formula for the π k is then

π (n+1) k = π (n) k + η (n+1) grad π k f (θ (n) ; x n+1 ) (2.23)
Where the coefficient η (n+1) denotes a step-size, and the vector x n+1 is a new data point (it can also be replaced by a mini-batch). From (2.23), the new weights w (n+1) k can be obtained according to (2.22).

The gradients with respect to the other parameters are

grad µ k f (θ; x) = -2o k ∂h β k (δ k ) ∂δ k Σ -1 k (x -µ k ) (2.24a) grad Σ k f (θ; x) = -o k 1 2 Σ k + ∂h β k (δ k ) ∂δ k (x -µ k )(x -µ k ) † (2.24b) grad β k f (θ; x) = o k ∂ log c(β k ) ∂β k + ∂h β k (δ k ) ∂β k (2.24c)
where

h β k (δ k ) = log g β k (δ k ) as in (1.1)
. Then, the Euclidean stochastic gradient updates for

θ k = (µ k , Σ k , β k ) are µ (n+1) k = µ (n) k + η (n+1) grad µ k f (θ (n) ; x n+1 ) (2.25a) Σ (n+1) k = Σ (n) k + η (n+1) grad Σ k f (θ (n) ; x n+1 ) (2.25b) β (n+1) k = β (n) k + η (n+1) grad β k f (θ (n) ; x n+1 ) (2.25c)
In [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF], the gradients (2.24) are used in the standard formulation of a stochastic gradient descent, (2.25). However, in [START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF], they are used in the framework of an online proximal alternating linearized minimisation method (PALM), applied to the estimation of mixtures of Student T-distributions. In this work, the cost function is reshaped by adding the following surjective mappings to keep β k ∈ (0, ∞) and Σ k ∈ P m after each iteration.

ϕ β k (β k ) = β 2 k + ϕ Σ k (Σ k ) = Σ † k Σ k + I m with > 0
In addition, an inertial gradient is applied to this PALM method to accelerate the iterative process, with a hyper parameter

ρ k ∈ [0, 1], for each sub-component θ k of θ θ(n) k = θ (n) k + ρ k (θ (n) k -θ (n) k ) θ (n+1) k = θ(n) k + η (n) grad θ k f (θ; x)

Riemannian methods

In the Riemannian context, the MLE of a mixture of multivariate Gaussian distributions was studied in [START_REF] Hosseini | An alternative to EM for gaussian mixture models: batch and stochastic Riemannian optimization[END_REF]. This work proposed an offline batch gradient method, and also an online mini-batch gradient method. The expectation µ k and covariance matrix Σ k of Gaussian model are integrated into a single parameter according to (2.14). The geometry of the new integrated matrix S k is considered based on the affine-invariant gradients (2.5) and the exponetial map (2.2). The affine-invariant gradient of the reformulated ˜ f (θ; x) with respect to S k is

∇ ↑ S k ˜ f (θ; x) = -o k 1 2 S k + ∂ hβ k (δ k (y)) ∂δ k (y) yy † (2.26)
where y † = (x, 1) † , and δ k (y) = y † S -1 k y. The o k are the membership probabilities, given in (2.21).

Moreover, the weights of the Gaussian mixture model have the same reformulation as in the in Euclidean method (2.22). Then, an alternating iteration scheme is adopted between the weights and the new parameter (S k ) 1 k K . The resulting Riemannian online method (stochastic gradient method) has the following general form.

1. Update π (n) k according to (2.23) 2. Reparametrize S (n) k according to (2.14) with current (µ (n) k , Σ (n) k ) 3. Compute the Riemannian gradient ∇ ↑ S k ˜ f (θ (n) ; x n+1 ) using (2.26) 4. Update S (n) k
according to the exponential map (2.2)

S (n+1) k = Exp S (n) k η (n+1) ∇ ↑ S k ˜ f (θ (n) ; x n+1 )
Possibly, the new datapoint x n+1 can be replaced with a mini-batch.

Decompose S

(n+1) k → (µ (n+1) k , Σ (n+1) k )
In the experiments of [START_REF] Hosseini | An alternative to EM for gaussian mixture models: batch and stochastic Riemannian optimization[END_REF], the Riemannian gradient method, offline or stochastic with mini-batch, is compared with other standard methods, such as EM-FP, conjugate gradient and BFGS. The Riemannian stochastic gradient method shows remarkable convergence behavior, making it a potential candidate for large scale mixture learning.

In addition to mixtures of multivariate Gaussians, an offline estimation method for MECD models with a fixed shape parameter was proposed in [START_REF] Li | A universal framework for learning the elliptical mixture model[END_REF]. This method is very similar to [START_REF] Hosseini | An alternative to EM for gaussian mixture models: batch and stochastic Riemannian optimization[END_REF], also employing the reparametrization (2.14). Then, it implements a Riemannian conjugate gradient rule, into the alternating iterative scheme of the form just described.

Conclusion

We summarize the existing estimation methods of ECD and MECD, reviewed in Sections 2. [START_REF] Amari | Methods of Information Geometry[END_REF] In the following chapter, we will introduce a new geometry for ECD and MECD models. This will be based on the original concept of component-wise information metric (CIM). Later on, in Chapters 5 and 6, the CIM will be used to define a new online estimation method, which we have called the component-wise information gradient method (CIG). This new method is one of the main contributions of the present thesis. It will be shown to have several significant advantages, in comparison with existing estimation methods, summarized here in Tables 2.1 and2

Introduction

Our contribution to the geometry of ECD and MECD lies in the introduction of the component-wise information metric (CIM). Our objective is to propose a new online estimation algorithm for ECD and their mixtures (that is, MECD). For ECD, if the shape parameter is known, efficient online estimation can be carried out using the Fisher information metric. However, if the shape parameter is unknown, the Fisher information metric does not have a closed-form expression. We specifically introduce the CIM as a computationally advantageous alternative to the Fisher information metric (FIM), in this case, and also in the more difficult case of MECD.

The aim of an online estimation algorithm is to find the true parameter θ * . Here, this is formulated as an optimisation problem on a Riemannian manifold, which is the parameter space Θ of an ECD, or Θ of an MECD model. This is the problem of minimising the Kullback-Leibler (or KL) divergence. For ECD models, the KL divergence is given by

D(θ) = p(x|θ * ) log p(x|θ * ) p(x|θ) dx = E θ * [log p(x|θ * )] -E θ * [log p(x|θ)] (3.1)
for θ ∈ Θ, where p(x|θ) is given by (1.1) and θ * is supposed to be the true parameter of the ECD model. Note that θ * is the unique global minimum of D(θ). Another way of formulating (3.1) is to say that any online estimation algorithm searches for the minimum

arg min θ∈Θ -E θ * [ p (θ; x)] with p (θ; x) = log p(x|θ) (3.2)
For MECD, the density function p and parameter θ from (1.1) are replaced by the density f and the parameter θ from (1.2), respectively. The KL divergence is given by

D(θ) = f (x|θ * ) log f (x|θ * ) f (x|θ) dx = E θ * [log f (x|θ * )] -E θ * [log f (x|θ)] (3.3)
for θ ∈ Θ. Then, online algorithms aim to search for the minimum

arg min θ∈Θ -E θ * [ f (θ; x)] with f (θ; x) = log f (x|θ) (3.4)
Here, unlike the ECD case, the true parameter θ * is a global minimum of D(θ), but it is not unique (this is because it is possible to permute mixture components). The target problems (3.2) and (3.4) are optimisation problems on the Riemannian manifolds Θ and Θ, respectively. The general update rule, for offline (batch) gradient optimisation methods on Riemannian manifolds, takes on the following form [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]

θ (n+1) = Ret θ (n) η (n) U (θ (n) ; X ) (3.5)
with U (θ (n) ; X ) the direction of descent. On the other hand, the general form of an update rule for online (stochastic) gradient optimisation methods is the following

θ (n+1) = Ret θ (n) η (n+1) U (θ (n) ; X n+1 ) (3.6)
Here, the smooth mapping Ret θ from the tangent space T θ Θ to Θ is required to be a retraction, in the sense that it verifies

Ret θ (0 θ ) = θ dRet θ (0 θ ) = Id T θ Θ (3.7)
where 0 θ denotes the zero element in T θ Θ, dRet θ (•) means the differential of retraction map Ret θ : T θ Θ → Θ, and Id T θ Θ denotes the identity mapping on T θ Θ. The positive scalars η (n) are step-sizes. For offline methods (3.5), each vector U (θ (n) ; X ) depends on the entire dataset X . For online methods (3.6), each vector U (θ (n) ; X n+1 ) depends only on a new mini-batch X n+1 , which may reduce to a single datapoint x n+1 . Ideally, U (θ (n) ; X ) or U (θ (n) ; X n+1 ) is the Riemannian information gradient, which can be computed, assuming the Fisher information metric is known.

However, in general, ECD and MECD models do not admit any tractable expression of the Fisher information metric. Therefore, there is no tractable means of evaluating the Riemannian information gradient.

In the present chapter, several specific situations are discussed. The Fisher information metric for θ = (Σ) (with known µ and β) has already been given in Subsection 2.2.2. In the following, Section 3.2 gives the Fisher information metric for θ = (µ, Σ) (with known β), and introduces a retraction map Ret for this case θ = (µ, Σ).

Sections 3.3 and 3.4, respectively, consider the most general cases, of complete ECD models and MECD models, where the shape parameter β is unknown. In these cases, the parameter space does not admit any tractable expression of the Fisher information metric [START_REF] Verdoolaege | On the geometry of multivariate generalized Gaussian models[END_REF]. Sections 3.3 and 3.4 introduce the component-wise information metric (CIM), as a computationally advantageous alternative to the Fisher information metric (FIM).

Geometry of (µ,Σ)

Fisher information metric of (µ, Σ): The Fisher information metric for the case θ = (Σ) (with known and fixed µ and β), has already been discussed, in Subsection 2.2.2.

Furthermore, when the location parameter µ is also unknown, but the shape parameter β is still known, the Fisher information metric for θ = (µ, Σ) can be computed explicitly. In fact, this metric takes on a simple form (this can be seen as a "direct product")

U θ , V θ * θ = U µ , V µ * µ + U Σ , V Σ * Σ (3.8)
where U µ , V µ * µ and U Σ , V Σ * Σ are Fisher information metrics in R m and P m respectively. Here, each tangent vector

U θ , V θ ∈ T θ Θ is of the form U θ = (U µ , U Σ ) † and V θ = (V µ , V Σ ) † .
Moreover, the second term U Σ , V Σ Σ is given by the same formula as in (2.6). For the first term, we have computed it to be [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF],

U µ , V µ * µ = I µ U † µ Σ -1 V µ (3.9)
where the information coefficient I µ is given by

I µ = - 4 m E ∂ 2 h(δ, β) ∂ δ 2 δ -2E ∂ h(δ, β) ∂ δ (3.10)
Adding up (2.6) and (3.9), one gets the full Fisher information metic for θ = (µ, Σ),

U θ , V θ θ = I µ U † µ Σ -1 V µ + I Σ,1 tr(Σ -1 U Σ Σ -1 V Σ ) + I Σ,2 tr(Σ -1 U Σ )tr(Σ -1 V Σ ) (3.11)
where the information constants I Σ,1 and I Σ,2 were given in equation (2.7). The Fisher information metric (3.11) generalises the one studied in the multivariate Gaussian case [START_REF] Besson | On the Fisher information matrix for multivariate elliptically contoured distributions[END_REF]. Retraction map of (µ, Σ): The Fisher information metric (3.11) on the product manifold Θ = R m ×P m , does not admit a tractable expression, for its Riemannian exponential map. However, when it is restricted to the subspace P m , its exponential map is well-known, and was given in equation (2.2). In addition, the exponential map on the subspace R m (a Euclidean space) reduces to vector addition. Accordingly, we proposed to use a retraction map which is defined by the product of these two exponentials [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF],

Ret θ : T θ Θ -→ Θ U θ = U µ U Σ -→ Exp µ (U µ ) Exp Σ (U Σ ) = µ + U µ Σ exp(Σ -1 U Σ ) (3.12)
Both of the exponential maps Exp µ and Exp Σ verify the properties (3.7). Therefore, their direct product (3.12) also verifies these properties, and is a well-defined retraction.

Information distance of (µ, Σ): We also introduced a distance function on the product manifold Θ = R m × P m , defined as the sum of squares of the distances in the subspaces R m and P m . This will be called the component-wise information distance (CID) [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF].

For P m , the information distance is defined as in Equation (2.8). For R m , the information distance is proportional to the usual Euclidean distance. Now, the component-wise information distance is the following distance on Θ [90],

d 2 * (θ 1 , θ 2 ) =d 2 * (µ 1 , µ 2 ) + d 2 * (Σ 1 , Σ 2 ) =I µ (µ 1 -µ 2 ) † (µ 1 -µ 2 ) + I Σ,1 tr log(Σ -1 1 Σ 2 ) 2 + I Σ,2 tr 2 log(Σ -1 1 Σ 2 ) (3.13)
where θ 1 = (µ 1 , Σ 1 ) and θ 2 = (µ 2 , Σ 2 ), with µ 1 , µ 2 ∈ R m and Σ 1 , Σ 2 ∈ P m . The constant I µ was given in (3.10), and the constants I Σ,1 and I Σ,2 were given in (2.7).

Information gradient of (µ, Σ): The information gradient associated to the Fisher information metric (3.11) is just the product of the information gradients, with respect to µ and Σ [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF],

∇ * θ p (θ; x) = ∇ * µ p (θ; x) ∇ * Σ p (θ; x) = I -1 µ Σ grad µ p (θ; x) J Σ,1 [∇ ↑ Σ p (θ; x)] ⊥ + J Σ,2 [∇ ↑ Σ p (θ; x)] (3.14)
where the coefficients J Σ,1 and J Σ,2 of ∇ * Σ p were given in equation (2.10), and the information constant I µ in (3.10). Here, the vector grad µ p (θ; x) is the gradient in the classical Euclidean sense

grad µ p (θ; x) = -2 ∂h β (δ) ∂δ Σ -1 (x -µ) (3.15)
and the affine-invariant gradient ∇ ↑ Σ p (θ; x) is given by Equation (2.5).

Geometry of (µ,Σ,β)

When the shape parameter β is also unknown, the Fisher information metric of the complete ECD model, with θ = (µ, Σ, β), does not have a closed form expression, and it is therefore impossible to derive applicable expressions of the Riemannian exponential map, the information gradient, or any other useful geometric objects.

A similar situation arises in the field of artificial neural networks, where the Fisher information metric is too difficult to compute. The Fisher information metric is then replaced with a quasi-diagonal reduction [START_REF] Marceau-Caron | Practical Riemannian neural networks[END_REF]. Motivated by this idea, from the field of artificial intelligence, we introduce the component-wise information metric of ECD models [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF].

Each one of the three parameters µ, Σ, and β has a closed form expression for its Fisher information metric. These closed form expressions give the three main diagonal blocks of the ECD model Fisher information metric. The other (off-diagonal) blocks being unknown, we just set them to zero. In this way, we obtain the component-wise information metric (CIM) of the ECD model, at θ = (µ, Σ, β),

I (θ) = V µ V Σ V β U µ I µ 0 0 U Σ 0 I Σ 0 U β 0 0 I β
The blocks I µ , I Σ and I β are respectively the Fisher information matrics of R m , P m and R + for ECD.

Component-wise information metric of (µ, Σ, β): The Fisher information metric for (µ, Σ) was given in Equation (3.11). The component-wise information metric is obtained by adding to (3.11) a further term for the shape parameter β. The CIM of the ECD model follows [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF] U

θ , V θ θ = U µ , V µ * µ + U Σ , V Σ * Σ + U β , V β * β =I µ U † µ Σ -1 V µ + I Σ,1 tr(Σ -1 U Σ Σ -1 V Σ ) + I Σ,2 tr(Σ -1 U Σ )tr(Σ -1 V Σ ) + I β U β V β (3.16)
where is the symbol for component-wise information metric,

U θ = (U µ , U Σ , U β ) and V θ = (V µ , V Σ , V β )
are tangent vectors at the point θ = (µ, Σ, β). The information constants I Σ,1 , I Σ,2 and I µ are given in (2.7) and (3.10), respectively. The information constant of β is

I β = -E ∂ 2 log c(β) ∂ β 2 + ∂ 2 h β (δ) ∂ β 2 (3.17)
Retraction map of (µ, Σ, β): We introduce a retraction map for (µ, Σ, β), using the same idea as in (3.12). Recall that the exponential map for µ is the vector addition, and the intrinsic exponential map on P m was given in equation (2.2). Since β belongs to R + , a 1-dimensional version of (2.2) is used as the exponential map for β.

Combining these three exponential maps, The retraction for the complete ECD model is given as [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF],

Ret θ : T θ Θ -→ Θ U θ =   U µ U Σ U β   -→   Exp µ (U µ ) Exp Σ (U Σ ) Exp β (U β )   =   µ + U µ Σ exp(Σ -1 U Σ ) βe U β /β   (3.18)
where U θ is the tangent vector (in the following, U θ will represent the direction of descent). All the three exponential maps Exp µ , Exp Σ and Exp β verify the properties (3.7). Therefore, their direct product (3.18) also verifies these properties, and is a well-defined retraction.

Component-wise information distance of (µ, Σ, β): The information distance of (µ, Σ) was given in (3.13). In order to include β, which belongs to (0, ∞), a onedimensional version of Equation (2.4) is added to (3.13). The component-wise information distance for the complete model (with all three parameters µ, Σ and β), then has the expression [90]

d 2 (θ 1 , θ 2 ) =d 2 * (µ 1 , µ 2 ) + d 2 * (Σ 1 , Σ 2 ) + d 2 * (β 1 , β 2 ) =I µ (µ 1 -µ 2 ) † (µ 1 -µ 2 ) + I Σ,1 tr log(Σ -1 1 Σ 2 ) 2 + I Σ,2 tr 2 log(Σ -1 1 Σ 2 ) + I β log 2 β -1 1 β 2 (3.19)
Where the information constants I µ , I Σ,1 , I Σ,2 and I β are given respectively in (3.10), (2.7) and (3.17).

Component-wise information gradient of (µ, Σ, β): Consider now the Componentwise Information Gradient (CIG for short), derived from the metric (3.16). The CIG is the unique vector field ∇ θ p (θ; x) on Θ which satisfies

d p (θ; x) [U θ ] = ∇ θ p (θ; x), U θ θ recall p = log p (3.20)
where the scalar product on the right-hand side is the CIM (3.16), and d p is the differential of p . Accordingly, the component-wise information gradient has the following form [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF]

∇ θ p (θ; x) =   ∇ * µ p (θ; x) ∇ * Σ p (θ; x) ∇ * β p (θ; x)   =   I -1 µ Σ grad µ p (θ; x) J Σ,1 [∇ ↑ Σ p (θ; x)] ⊥ + J Σ,2 [∇ ↑ Σ p (θ; x)] I -1 β grad β p (θ; x)   (3.21)
where the coefficients I µ , J Σ,1 , J Σ,2 and I β are respectively given in (2.7), (3.10), and (3.17). For the third component in (3.21), the derivative of p with respect to β is

grad β p (θ; x) = ∂ log c(β) ∂β + ∂h β (δ) ∂β (3.22)

Mixtures of ECD

A mixture of ECD is given by the probability density f (x|θ) in (1.2). This density is parameterised by θ

= (w k , θ k ; k = 1, • • • , K) where θ k = (µ k , Σ k , β k )
are the parameters of mixture components p(x|θ k ), which are single ECD of the form (1.1). In this thesis, the mixture weights w = (w k ) k are mapped to a point r = (r k ) k on the unit sphere S K-1 ⊂ R K , by setting w k = r 2 k . In this way, the parameter space Θ of the MECD model is a product manifold, with subspaces S K-1 (for the weights), and Θ k = R m × P m × R + (for each single ECD mixture component). In other words,

Θ = S K-1 × Θ 1 × . . . × Θ K Θ k = R m × P m × R +
This parameter space will be given the structure of a d-dimensional Riemannian manifold (the dimension d is given by d = K(m(m + 3)/2 + 1)).

Ideally, one hopes to equip Θ with the Fisher information metric, and derive the corresponding Riemannian information gradient and exponential map. However, as in the single ECD case (of the previous section), this Fisher information metric does not have a closed form expression, and the same is true of the corresponding exponential map and information gradient.

Therefore, as in the previous section, a component-wise information metric, of quasidiagonal form, is introduced on Θ,

I (θ) = dr V µ k V Σ k V β k       U r I r 0 0 0 U µ k 0 I µ k 0 0 U Σ k 0 0 I Σ k 0 U β k 0 0 0 I β k Its explicit expression is [89] U θ , V θ θ = U r , V r * r + K k=1 U µ k , V µ k * µ k + K k=1 U Σ k , V Σ k * Σ k + K k=1 U β k , V β k * β k (3.23)
where U θ , V θ are tangent vectors, in the tangent space T θ Θ, and all sub-metrics are the Fisher information metrics in the corresponding subspaces. Precisely, the Fisher information metric •, • * r on unit sphere coincides with the scalar product in Euclidean sens [START_REF] Amari | Information geometry and its applications[END_REF]. For each mixture component, the Fisher information metric in each subspace has the same form as in Equation (3.16).

U µ k , V µ k * µ k = I µ k U † µ k Σ -1 k V µ k (3.24a) U Σ k , V Σ k * Σ k = I Σ k ,1 tr(Σ -1 k U Σ k Σ -1 k V Σ k ) + I Σ k ,2 tr(Σ -1 k U Σ k )tr(Σ -1 k V Σ k ) (3.24b) U β k , V β k * β k = I β k U β k V β k (3.24c)
The information constants in Equations (3.24) are also equal to the values given above.

I Σ k ,1 = 2A m(m + 2) I Σ k ,2 = A m(m + 2) - 1 4 A = E ∂h β k (δ k ) ∂δ k δ k 2 (3.25a) I µ k = - 4 m E ∂ 2 h β k (δ k ) ∂ δ 2 k δ k -2E ∂ h β k (δ k ) ∂ δ k (3.25b) I β k = -E ∂ 2 log c(β k ) ∂ β 2 k + ∂ 2 h β k (δ k ) ∂ β 2 k (3.25c)
where

δ k = (x -µ k ) † Σ -1 k (x -µ k ).
To define a retraction on Θ, the retraction map in (3.18) is used for each component θ k .

For the weights, as given by the parameters r = (r k ) k ∈ S K-1 , the intrinsic Riemannian exponential on the sphere S K-1 ⊂ R K has an easy expression

Exp r (U r ) = cos ( U r ) r + sin ( U r ) U r U r (3.26)
The direct product of these maps is chosen as the retraction for MECD [START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF],

Ret θ : T θ Θ → Θ U θ =     U r (U µ k ) k (U Σ k ) k (U β k ) k     →     Exp r (U r ) Exp µ k (U µ k ) k Exp Σ k (U Σ k ) k Exp β k (U β k ) k     (3.27)
here, for any component k, following (3.18),

  Exp µ k (U µ k ) Exp Σ k (U Σ k ) Exp β k (U β k )   =   µ k + U µ k Σ k exp(Σ -1 k U Σ k ) β k e U β k /β k  
The component-wise information gradient (CIG), derived from the metric (3.23), is given as [START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF],

∇ θ f (θ; x) =     ∇ * r f (θ; x) ∇ * µ k f (θ; x) k ∇ * Σ k f (θ; x) k ∇ * β k f (θ; x) k     (3.28)
Each element ∇ * i f (θ; x) in this column vector is the information gradient with respect to the corresponding i-th subspace. Concretely, the information gradient with respect to r is

∇ * r f (θ; x) = ∂ ∂r f (θ; x) - ∂ ∂r f (θ; x) , r × r (3.29)
where •, • denotes the scalar product in R K (the dot product). The information gradients for µ k , Σ k and β k are consistent with (3.21)

  ∇ * µ k f (θ; x) ∇ * Σ k f (θ; x) ∇ * β k f (θ; x)   =   I -1 µ k Σ k grad µ k f (θ; x) J Σ k ,1 ∇ ↑ Σ k f (θ; x) ⊥ + J Σ k ,2 ∇ ↑ Σ k f (θ; x) I -1 β k grad β k f (θ; x)  
Here, all necessary information constants are defined in (3.25). The vectors grad µ k f (θ; x), ∇ ↑ Σ k f (θ; x) and grad β k f (θ; x) are defined in (2.24) and (2.5).

Conclusion

Our main contribution to the geometry of ECD and MECD is the introduction of the component-wise information metric (CIM), in Sections 3.3 and 3.4 of the above. The CIM coincides with the Fisher information metric (FIM) whenever the shape parameter is known (when θ = (Σ) or θ = (µ, Σ)). When the shape parameter is unknown (when θ = (µ, Σ, β)), the CIM has a straightforward "direct product" structure, while the FIM does not admit any closed-form expression.

The following Table 3.5 details the geometric objects introduced in the present chapter, for each special case of ECD or MECD models.

Model Summary

ECD with θ = (Σ)

• the CIM coincides with the FIM (2.6)

• the information gradient is given by (2.9)

• Exp is tractable, given by (2.2)

ECD with θ = (µ, Σ)

• the CIM coincides with the FIM (3.11)

• the information gradient is given by (3.14)

• use Ret given by (3.12), instead of intractable Exp ECD with θ = (µ, Σ, β)

• the CIM is different from the FIM, and given by (3.16)

• the component-wise information gradient is given by (3.21)

• use Ret given by (3.18), instead of intractable Exp

MECD with θ = (w, θ 1 , • • • , θ K )
• the CIM is different from the FIM, and given by (3.23)

• the component-wise information gradient is given by (3.28)

• use Ret given by (3.27), instead of intractable Exp Table 3.1: The CIM and related geometric objects for ECD and MECD

The above table shows that there are certain special cases where the CIM coincides with the FIM, which is known in closed form. In such cases, the Riemannian information gradient can be computed explicitly, and replaced into the general update rules (3.5) or (3.6). The following Chapter 4 will be exactly concerned with such cases, and studies Riemannian online estimation of a general statistical model, whose Fisher information metric is known. More difficult cases, where the CIM is used, because the FIM is unknown in closed form, will be the subject of subsequent chapters 5 and 6.

Chapter 4

Riemannian information gradient method

Introduction

The present chapter reflects our work, published in the journal paper [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF]. Its central theme is the application of the Riemannian information gradient method, to online (also called recursive or stochastic) estimation of statistical models parameterised on Riemannian manifolds.

For these statistical models, the information gradient method can be applied, under the assumption that the Fisher information metric is known in closed form. Then, the information gradient can be computed directly from the Fisher information metric. The information gradient method greatly simplifies the task of selecting step-sizes. Moreover, this method guarantees fast convergence, asymptotic normality, and asymptotic efficiency of online estimation.

Unlike other chapters of this thesis, the present chapter does not focus solely on ECD models, but considers general statistical models, whose parameters belong to a Riemannian manifold. After stating general results on the Riemannian information gradient method, it returns to certain special cases of ECD, in order to illustrate these results.

The focus will be on an online update rule which falls under the general form (3.6). The aim of this update rule is to minimise a statistical divergence function D(θ) (such as the Kullback-Leibler divergence), defined on some Riemannian manifold Θ. This minimisation is equivalent to finding the true value θ * of the statistical parameter θ, as already discussed in Section 3.1, before equation (3.1).

Specifically, the general online update rule considered in the following is (4.1a). The Riemannian information gradient method is a special case of this general rule, which will be given by (4.5). The main problem considered in the present chapter is how to choose the step-sizes in (4.1a), in order to obtain fast and asymptotically efficient online estimates. In solving this problem, several original results are introduced.

First, under mild assumptions on the divergence function, it is proved that, with an adequate choice of step-sizes, the update rule computes online estimates which achieve a fast rate of convergence.

Second, the asymptotic normality of these online estimates is proved, by employing a novel linearisation technique (see Appendix A.1.4).

Third, it is proved that the Riemannian information gradient method is asymptotically efficient, in the sense that it achieves an optimal asymptotic rate of convergence.

These results, while relatively familiar in the Euclidean context (for example, see [START_REF] Duflo | Random iterative models[END_REF][START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF]), were formulated and proved for the first time, in a Riemannian context, in our paper [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF]. They will be illustrated with a numerical application to the online estimation of the scatter matrix of an ECD model, with known location and shape parameters.

The mathematical problem, considered in the present chapter, is formulated in Section 4.2. This involves a parameterised statistical model P of probability distributions P θ , where the statistical parameter θ belongs to a Riemannian manifold Θ. Given independent observations, with distribution P θ * for some θ * ∈ Θ, the aim is to estimate the unknown true parameter θ * .

In principle, this is done by minimising a statistical divergence function D(θ), which measures the dissimilarity between P θ and P θ * . Taking advantage of the observations, there are two approaches to minimising D(θ) : stochastic minimisation, which leads to online estimation, and empirical minimisation, which leads to classical techniques, such as maximum-likelihood estimation [START_REF] Broniatowski | Minimum divergence estimators, maximum likelihood and exponential families[END_REF][START_REF] Broniatowski | Parametric estimation and tests through divergences and the duality technique[END_REF].

The original results, obtained in the present chapter are stated in Section 4.3. In particular, these are Propositions 2, 4, and 5. Overall, these propositions show that online estimation, which requires less computational resources than maximum-likelihood estimation, can still achieve the same optimal performance, characterised by asymptotic efficiency. Recall that asymptotic efficiency means the asymptotic Cramér-Rao lower bound is achieved [START_REF] Ibragimov | Statistical estimation-asymptotic theory[END_REF][START_REF] Aad | Asymptotic statistics[END_REF].

To summarise these propositions, consider a sequence of online estimates θ (n) , computed using a Riemannian online update rule of the form (3.6). Informally, under assumptions which guarantee that θ * is an attractive local minimum of D(θ), and that the updates are neither too noisy, nor too unstable, in the neighborhood of θ * ,

• Proposition 2 states that, with an adequate choice of step-sizes, the θ (n) achieve a fast rate of convergence to θ * . Precisely, the expectation of the squared Riemannian distance between θ (n) and θ * is O (n -1 ). This is called a fast rate, because it is the best achievable, for any step-sizes which are proportional to n -q with q ∈ (1/2, 1] [7,[START_REF] Duflo | Random iterative models[END_REF].

• Proposition 4 states that the distribution of the θ (n) becomes asymptotically normal, centred at θ * , when n grows increasingly large, and also characterises the corresponding asymptotic covariance matrix.

• Proposition 5 is concerned with the Riemannian information gradient method. If the Riemannian manifold Θ is equipped with the Fisher information metric of the statistical model P , then Riemannian information gradient descent computes online estimates θ (n) which are asymptotically efficient. This is illustrated, with a numerical application to the online estimation of the scatter matrix of an ECD model, in Section 4.

The end-result of Proposition 5 is that the Riemannian information gradient method, which uses the Fisher information metric, achieves asymptotic efficiency. In [START_REF] Tripuraneni | Averaging stochastic gradient descent on Riemannian manifolds[END_REF], an alternative route to asymptotic efficiency is proposed, using the averaged Riemannian stochastic gradient method. This method does not require any prior knowledge of the Fisher information metric, but has an additional computational cost, which comes from computing on-line Riemannian averages.

The proofs of Propositions 2, 4, and 5 are detailed in Appendices A.1, A.2 and A.3.

Problem statement

Let P = (P, Θ, X) be a statistical model, with parameter space Θ and sample space X.

To each θ ∈ Θ, the model P associates a probability distribution P θ on X. Here, we consider the case where Θ is a Riemannian manifold, and X is any measurable space.

The Riemannian metric of Θ will be denoted •, • , with its Riemannian distance d (•, •). In general, the metric •, • is not the information metric of the model P .

Let (Ω, F, P) be a complete probability space, and (x n ; n = 1, 2, • • • ) be i.i.d. random variables on Ω, with values in X. While the distribution of x n is unknown, it is assumed to belong to the model P . That is, P • x -1 n = P θ * for some θ * ∈ Θ, which is the true parameter.

Consider the following problem : how to obtain fast, asymptotically efficient, online estimates θ (n) of the true parameter θ * , based on observations of the random variables x n ? The present work proposes to solve this problem through a detailed study of the decreasing-step-size online update rule,

θ (n+1) = Exp θ (n) η (n+1) U (θ (n) , X new mb ) n = 0, 1, • • • (4.1a)
starting from an initial guess θ (0) . Here, X new mb could be just one sample or a mini-batch. When it is one sample, it is also denoted as x n . This update rule is similar to the one in [START_REF] Bonnabel | Stochastic gradient descent on Riemannian manifolds[END_REF], and falls under the general form (3.6). Instead of the general retraction Ret in (3.6), it uses the Riemannian exponential Exp of the metric •, • of Θ. In the following, the step-sizes η (n) are strictly positive, decreasing, and verify the usual conditions for stochastic approximation [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF][START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF]

η (n) = ∞ (η (n) ) 2 < ∞ (4.1b)
Moreover, U (θ, x) is a continuous vector field on Θ for each x ∈ X, which generalises the classical concept of score statistic [START_REF] Christopher C Heyde | Quasi-likelihood and its application: a general approach to optimal parameter estimation[END_REF][START_REF] Ibragimov | Statistical estimation-asymptotic theory[END_REF].

A priori knowledge about the model P is injected into the update rule (4.1a) using a divergence function D(θ). As defined in [START_REF] Amari | Methods of Information Geometry[END_REF], this is a positive function, equal to zero if and only if P θ = P θ * , and with positive-definite Hessian at θ = θ * . Since one expects that minimising D(θ) will lead to estimating θ * , it is natural to require that

E θ * U (θ, x) = -∇D(θ) (4.1c)
In other words, that U (θ, x) is an unbiased estimator of minus the Riemannian gradient of D(θ). With U (θ, x) given by (4.1c), the update rule (4.1a) is a Riemannian stochastic gradient method, of the form considered in [START_REF] Bonnabel | Stochastic gradient descent on Riemannian manifolds[END_REF][START_REF] Tripuraneni | Averaging stochastic gradient descent on Riemannian manifolds[END_REF][START_REF] Zhang | First-order methods for geodesically convex optimization[END_REF].

In practice, a suitable choice of D(θ) is often the Kullback-Leibler divergence [START_REF] Shiryaev | Probability[END_REF],

D(θ) = -E θ * [ (θ; x)] (θ; x) = log p θ (x) p θ * (x) (4.2)
where it is assumed that P θ and P θ * have common support, with probability density functions p θ (x) and p θ * (x), respectively. Indeed, if D(θ) is chosen to be the Kullback-Leibler divergence, then (4.1c) is satisfied by

U (θ, x) = ∇ (θ; x) (4.3)
which, in many practical situations, can be evaluated directly.

Riemannian information gradient method:

When the Fisher information metric •, • * of the model P is known, it is natural to equip Θ with this Fisher information metric. Then, the Riemannian gradient on the right-hand side of equation ( 4.3) coincides with the Riemannian information gradient,

U (θ, x) = ∇ * (θ; x) (4.4)
If this expression of U (θ, x) is replaced into equation (4.1a), this equation becomes

θ (n+1) = Exp θ (n) η (n+1) ∇ * θ (θ (n) ; X new mb ) (4.5)
which will be called the Riemannian information gradient method.

Before stating the main results, in the following Section 4.3, it may be helpful to recall some general background, on online estimation [START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF]. For simplicity, let D(θ) be the Kullback-Leibler divergence (4.2). The problem of estimating the true parameter θ * is equivalent to the problem of finding a global minimum of D(θ). Of course, this problem cannot be tackled directly, since D(θ) cannot be computed without knowledge of θ * . There exist two routes, around this difficulty.

The first route is empirical minimisation, which replaces the expectation in (4.2) with an empirical mean over observed data. Given the first N observations, instead of minimising D(θ), one minimises the empirical divergence D(θ),

D(θ) = - 1 N N n=1 (θ; x n ) (4.6)
where (θ; x) is the likelihood function of (4.2). Now, given the minus sign ahead of the sum in (4.6), it is clear that minimising D(θ) amounts to maximising the sum of log-likelihoods. Thus, this lead to the method of maximum-likelihood estimation.

It is well-known that maximum-likelihood estimation, under general regularity conditions, is asymptotically efficient [START_REF] Ibragimov | Statistical estimation-asymptotic theory[END_REF]. Roughly, this means the maximum-likelihood estimator has the least possible asymptotic variance, equal to the inverse of the Fisher information matrix. On the other hand, as the number N of observations grows very large, it can be especially difficult to deal with the empirical divergence D(θ) of (4.6). In the process of searching for the minimum of D(θ), each evaluation of this function, or of its derivatives, will involve a massive number of operations, ultimately becoming unpractical.

Aiming to avoid this difficulty, the second route, that applies to online estimation, is based on observation-driven updates, following the general scheme of (4.1a). In this scheme, each new online estimate θ (n+1) is computed using only one new observation x n+1 , or eventually a mini-batch. Therefore, as the number of observations grows very large, the overall computational effort required by online estimation, remains the same.

The main results in the following section show that online estimation can achieve the same asymptotic performance as maximum-likelihood estimation, as the number N of observations grows large. As a word of caution, it should be said that online estimation does not, in general, fare better than maximum-likelihood estimation for moderate values of the number N of observations, and models with a small number of parameters. However, it is a very desirable substitute to maximum-likelihood estimation for models with a large number of parameters, which typically require a very large number of observations, in order to be estimated correctly.

Main results

The motivation of the following Propositions 1 to 5 is to provide general conditions, which guarantee that the update rule (4.1a) computes fast, asymptotically efficient, online estimates θ (n) of the true parameter θ * . In the statement of these propositions, it is implicitly assumed that conditions (4.1b) and (4.1c) are verified. Moreover, the following assumptions are considered.

(d1) the divergence function D(θ) has an isolated stationary point at θ = θ * , and Lipschitz gradient in a neighborhood of this point.

(d2) this stationary point is moreover attractive : D(θ) is twice differentiable at θ = θ * , with positive-definite Hessian at this point.

(u1) in a neighborhood of θ = θ * , the function

V (θ) = E θ * U (θ, x) 2 θ is uniformly bounded. (u2) in a neighborhood of θ = θ * , the function R(θ) = E θ * U (θ, x) 4
θ is uniformly bounded.

For Assumption (d1), the definition of a Lipschitz vector field on a Riemannian manifold may be found in [START_REF] Munier | Steepest descent method on a Riemannian manifold: the convex case[END_REF]. For Assumptions (u1) and (u2), • denotes the Riemannian norm.

Assumptions (u1) and (u2) are so-called moment control assumptions. They imply that the noisy nature of the observations x n does not cause the iterates θ (n) to diverge, and are also crucial to proving the asymptotic normality of these iterates.

Let Θ * be a neighborhood of θ * which verifies (d1), (u1), and (u2). Without loss of generality, it is assumed that Θ * is compact and convex (see the definition of convexity in [START_REF] Petersen | Riemannian geometry[END_REF][START_REF] Shiryaev | Probability[END_REF]). Then, Θ * admits a system of normal coordinates (θ α ; α = 1, • • • , d) with origin at θ * . With respect to these coordinates, denote the components of U (θ * , x) by u α (θ * ) and let Σ * = (Σ * αβ ),

Σ * αβ = E θ * u α (θ * ) u β (θ * ) (4.7)
When (d2) is verified, denote the components of the Hessian of D(θ) at θ = θ * by H = (H αβ ),

H αβ = ∂ 2 D ∂θ α ∂θ β θ α =0 (4.8)
Then, the matrix H = (H αβ ) is positive-definite [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. Denote by λ > 0 its smallest eigenvalue.

Propositions 1 to 5 require the condition that the estimates θ (n) are stable, which means that all the θ (n) lie in Θ * , almost surely. The need for this condition is discussed in Remark 3. Note that, if θ (n) lies in Θ * , then θ (n) is determined by its normal coordinates (θ (n) ) α .

Proposition 1 (consistency) assume (d1) and (u1) are verified, and the estimates θ (n) are stable. Then, lim θ (n) = θ * almost surely.

Proposition 2 (mean-square rate) assume (d1), (d2) and (u1) are verified, the estimates θ (n) are stable, and η (n) = a n where 2λa > 1. Then

E d 2 (θ (n) , θ * ) = O n -1 (4.9)
Proposition 3 (almost-sure rate) assume the conditions of Proposition 2 are verified. Then, d 2 (θ (n) , θ * ) = o(n -p ) for p ∈ (0, 1) almost surely (4.10)

Proposition 4 (asymptotic normality) assume the conditions of Proposition 2, as well as assumption (u2), are verified. Then, the distribution of the re-scaled coordinates [n 1/2 (θ (n) ) α ] converges to a centred d-variate normal distribution, with covariance matrix Σ given by Lyapunov's equation

A Σ + Σ A = -a 2 Σ * (4.11)
where A = (A αβ ) with A αβ = 1 2 δ αβ -aH αβ (here, δ denotes Kronecker's delta). Proposition 5 (asymptotic efficiency) assume the Riemannian metric •, • of Θ coincides with the information metric of the model P , and let D(θ) be the Kullback-Leibler divergence (4.2). Further, assume (d1), (d2), (u1) and (u2) are verified, the online estimates θ (n) are stable, and η (n) = a n where 2a > 1. Then, (i) the rates of convergence (4.9) and (4.10) hold true.

(ii) if a = 1, the distribution of the re-scaled coordinates [n 1/2 (θ (n) ) α ] converges to a centred d-variate normal distribution, with covariance matrix Σ * .

(iii) if a = 1, and U (θ, x) is given by (4.4), then Σ * is the identity matrix, and the online estimates θ (n) are asymptotically efficient.

(iv) the following rates of convergence also hold

E D(θ (n) ) = O(n -1 ) (4.12a) D(θ (n) ) = o(n -p ) for p ∈ (0, 1)
almost surely (4.12b)

Remark 1 assumptions (d2), (u1) and (u2) do not depend on the Riemannian metric •, • of Θ. Precisely, if they are verified for one Riemannian metric on Θ, then they are verified for any Riemannian metric on Θ. Moreover, if the function D(θ) is C 2 , then the same is true for assumption (d1). In this case, Propositions 1 to 5 apply for any Riemannian metric on Θ, so that the choice of the metric •, • is a purely practical matter, to be decided according to applications.

Remark 2 the conclusion of Proposition 1 continues to hold, if (4.1c) is replaced by

E θ * U (θ, x),∇D(θ) < 0 for θ = θ * (4.13)
Then, it is even possible to preserve Propositions 2, 3, and 4, provided assumption (d2) is replaced by the assumption that the mean vector field, X(θ) = E θ * U (θ, x), has an attractive stationary point at θ = θ * . This generalisation of Propositions 1 to 4 can be achieved following essentially the same approach as laid out in Section 4.3. However, in the present work, it will not be carried out in detail.

Remark 3 the condition that the online estimates θ (n) are stable is standard in all prior work on stochastic optimisation in manifolds [START_REF] Bonnabel | Stochastic gradient descent on Riemannian manifolds[END_REF][START_REF] Tripuraneni | Averaging stochastic gradient descent on Riemannian manifolds[END_REF][START_REF] Zhang | First-order methods for geodesically convex optimization[END_REF]. In practice, this condition can be enforced through replacing (4.1a) by a so-called projected or truncated update rule. This is identical to (4.1a), except that θ (n) is projected back onto the neighborhood Θ * of θ * , whenever it falls outside of this neighborhood [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF][START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF]. On the other hand, if the θ (n) are not required to be stable, but (d1) and (u1) are replaced by global assumptions, (d1') D(θ) has compact level sets and globally Lipschitz gradient.

(u1') V (θ) C (1 + D(θ)) for some constant C and for all θ ∈ Θ.

then, applying the same arguments as in the proof of Proposition 1, it follows that the θ (n) converge to the set of stationary points of D(θ), almost surely.

Remark 4 from (ii) and (iii) of Proposition 5, it follows that the distribution of n d 2 (θ (n) , θ * ) converges to a χ 2 -distribution with d degrees of freedom. This provides a practical means of confirming the asymptotic efficiency of the online estimates θ (n) .

Application to scatter-matrix estimation

Here, the conclusion of Proposition 5 is illustrated, by applying the update rule (4.1a) to the estimation of the scatter matrix of an ECD, with fixed location and shape parameters. Precisely, in the notation of Section 4.2, let Θ = P m the space of m × m positivedefinite matrices, and X = R m . For each θ ∈ Θ, let P θ be the ECD distribution with density p(x|θ) given by (1.1), with location parameter µ = 0, scatter matrix Σ = θ, and known shape parameter

β = β * , p(x|θ) = c(β * ) 1 det(θ) g β * (δ) (4.14)
where δ = x † θ -1 x. In the following, the familiar notation h β (δ) = log g β (δ) is also used. In addition, let (x (n) ; n = 1, 2, • • • ) be i.i.d. random vectors in R m , with distribution P θ * , where θ * ∈ Θ is the true scatter matrix. The standard approach to estimating the true scatter matrix θ * is based on maximum-likelihood estimation [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Sra | Conic geometric optimization on the manifold of positive definite matrices[END_REF]. An original approach, based on online estimation, is now introduced using the update rule (4.1a).

As in Proposition 5, the parameter space Θ = P m is equiped with the Fisher information metric (2.6). Here, this can be written as follows

U, V * θ = I θ,1 tr θ -1 U θ -1 V + I θ,2 tr θ -1 U tr θ -1 V (4.15)
for U, V ∈ T θ Θ, where the information coefficients I θ,1 and I θ,2 are given by (2.7),

I θ,1 = 2A m(m + 2) I θ,2 = A m(m + 2) - 1 4 A = E θ ∂h β * (δ) ∂δ δ 2 (4.16)
Note that these coefficients I θ,1 and I θ,2 do not depend on θ, but only on β * . The Fisher information metric with respect to θ being given by (4.15), it becomes possible to specify the update rule (4.1a). Indeed, recall from Section 2.2 that the exponential map of the metric (4.15) is given by (2.2),

Exp θ (u) = θ exp θ -1 u (4.17a)
Then, according to (ii) of Proposition 5, consider the step-sizes

η (n) = 1 n (4.17b)
Finally, according to (iii) of Proposition 5, let U (θ, x) be equal to the Riemannian information gradient,

U (θ, x) = ∇ * (θ; x) (4.17c)
This is the gradient of (θ; x) = log p(x|θ) with respect to the Fisher information metric (4.15), and can be computed from (2.9). Now, replacing (4.17) into (4.1a) defines an original update-rule (i.e. an original algorithm) for online estimation of the true scatter matrix θ * . (4.15). This distance is given by (2.8)

d 2 * (θ, θ * ) = I θ,1 tr log θ -1 θ * 2 + I θ,2 tr 2 log θ -1 θ * θ , θ * ∈ Θ (4.18)

Conclusion

This chapter introduced the Riemannian information gradient method, in the form of the update rule (4.5), applying it to online estimation of statistical models parameterised on Riemannian manifolds. The Riemannian information gradient method leads to a fast rate of convergence, and even to asymptotic efficiency. However, it can only be applied if the Fisher information metric is known.

The main results obtained in this chapter describe the convergence behavior of online estimates θ (n) , computed using the the general decreasing-step-size online update rule (4.1a), under Assumptions (d1), (d2), (u1), and (u2), formulated at the beginning of Section 4.3.

• Proposition 1 states that the θ (n) converge almost surely to the true parameter θ * , under Assumptions (d1) and (u1). The rate of almost sure convergence is stated in Proposition 3, with the additional Assumption (d2).

• Proposition 2 states the mean-square rate of convergence of the θ (n) . If the step-sizes η (n) are chosen correctly, this becomes a fast rate O(n -1 ). This proposition requires Assumptions (d1), (d2), and (u1).

• Propositions 4 states that the asymptotic distribution of the θ (n) is a normal distribution, centred at θ * , and characterises its covariance matrix. This proposition requires Assumptions (d1), (d2), (u1), and (u2).

• Proposition 5 states that, when the Fisher information metric is used to implement the update rule (4.1a) (this corresponds to the Riemannian information gradient method), there is an automatic choice of step-sizes (η (n) = 1/n) which guarantees the θ (n) are asymptotically efficient.

Section 4.4 applies Proposition 5 to the estimation of the scatter matrix of an ECD, with known location and shape parameters. It is impossible to extend this application to the case of a general ECD model (with unknown shape parameter), since the Fisher information metric is unknown, in this case. In order to overcome this difficulty, the following Chapter 5 will consider the component-wise information gradient method (CIG method).

Chapter 5

Estimation of ECD with CIG

Introduction

The main body of this chapter is based on our article [START_REF] Zhou | Riemannian information gradient methods for the parameter estimation of ECD: Some applications in image processing[END_REF], currently under review.

As shown in the previous chapter, the Fisher information metric and the related information gradient have significant advantages. The Riemannian information gradient method has a fast rate of convergence, and achieves asymptotic efficiency, with a simple (universal) choice of step-sizes (η (n) = 1/n).

For a complete ECD model, with unknown parameter θ = (µ, Σ, β), including location and shape parameters, the Fisher information metric does not have a closed form expression, so the Riemannian information gradient method cannot be applied.

Here, the Component-wise Information Metric (CIM), introduced in Chapter 3, will be used instead of the Fisher information metric (FIM). The gradient with respect to the CIM is the component-wise information gradient (CIG), which will now be used to define the CIG method for estimation of the complete ECD model.

The CIM is a quasi-diagonal reduction of the FIM. Thus, the CIG method retains some of the good properties of the Riemannian information gradient method (fast convergence and asymptotic efficiency, as seen in Proposition 5 of the previous Chapter 4), but not all of them. For example, the CIG method can still achieve a fast rate of convergence when estimating the true parameter θ * = (µ * , Σ * , β * ), but will fail to achieve asymptotic efficiency (theoretical results about the CIG method are given in Section 5.3).

The CIG method has an offline version, which can be applied to moderate-sized datasets, and an online version, more suited for recursive computations, which may be needed in applications such as change detection [START_REF] Bouchard | Riemannian geometry for compound Gaussian distributions: Application to recursive change detection[END_REF]. Both of these methods compute updated estimates (µ (n+1) , Σ (n+1) , β (n+1) ) from current estimates (µ (n) , Σ (n) , β (n) ), based on the alternating scheme step 1 :

µ (n+1) ← (µ (n) , Σ (n) , β (n) ) step 2 : Σ (n+1) ← (µ (n+1) , Σ (n) , β (n) ) step 3 : β (n+1) ← (µ (n+1) , Σ (n+1) , β (n) ) (5.1)
where each sub-parameter µ, Σ, or β is updated separately, in its own turn.

In the following, Section 5.2 starts with offline version of the CIG method. The online version and its theoretical convergence properties are introduced in Section 5.3. The question of the geodesic convexity of the Kullback-Leibler divergence D(θ), in the case of ECD estimation (recall Equation (3.1)), is studied in Section 5.4.

CIG offline

The CIG offline method is a second-order deterministic gradient method, somewhat similar to a Newton method. In the Newton method, the direction of descent is found by solving the Newton equation [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. In the CIG offline method, the Hessian in the Newton equation is approximated by the component-wise information metric (or matrix) I (θ), from Section 3.3.

For the CIG offline method, the choice of update direction depends on the complete dataset. The cost function (3.1) is reformulated, by replacing the KL divergence, with the empirical average ofp (θ; x), as in equation (4.6). This empirical average will be denoted by D(θ).

If the current estimate is θ (n) = (µ (n) , Σ (n) , β (n) ), the direction of update is given by one of the three components (depending on which sub-parameter is being updated)

∇ * µ D(θ (n) ) = - 1 T T t=1 ∇ * µ p (θ (n) ; x t ) (5.2a) ∇ * Σ D(θ (n) ) = - 1 T T t=1 ∇ * Σ p (θ (n) ; x t ) (5.2b) ∇ * β D(θ (n) ) = - 1 T T t=1 ∇ * β p (θ (n) ; x t ) (5.2c)
where the vectors ∇ * µ p , ∇ * Σ p and ∇ * β p under the empirical averages are given in (3.21). Using the expressions in (5.2), the CIG offline algorithm can now be stated as follows.

Algorithm 1 CIG offline algorithm for ECD estimation

Input: A finite dataset X , an initialization θ (0) = (µ (0) , Σ (0) , β (0) ); Output: The estimate θ;

1: for n = 0, 1, 2 • • • , N do 2:
Update µ:

3: θ current ← (µ (n) , Σ (n) , β (n) ); 4: Compute ∇ * µ D(θ current ); 5:
η µ ← Armijo-Goldstein(X , θ current );

6:

µ (n+1) ← Exp µ (n) -η µ ∇ * µ D(θ current ) ; 7:
Update Σ:

8: θ current ← (µ (n+1) , Σ (n) , β (n) ); 9:
Compute ∇ * Σ D(θ current ); 10:

η Σ ← Armijo-Goldstein(X , θ current );

11:

Σ (n+1) ← Exp Σ (n) -η Σ ∇ * Σ D(θ current ) ; 12:
Update β:

13:

θ current ← (µ (n+1) , Σ (n+1) , β (n) ); 14:
Compute ∇ * β D(θ current );

15:

η β ← Armijo-Goldstein(X , θ current );

16:

β (n+1) ← Exp β (n) -η β ∇ * β D(θ current ) ; 17: end for 18: θ ← (µ (N +1)) , Σ (N +1) , β (N +1) );
In this algorithm, the three exponential maps Exp µ , Exp Σ and Exp β are defined in (3.18), and the gradients are computed based on (5.2). The constant η denotes the step-size, which is selected according to the Armijo-Goldstein rule, thanks to a backtracking procedure (Definition 4.2.2. in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]). The following Proposition 6 states the convergence of Algorithm 1.

In this proposition, θ * denotes a stationary point of the cost function D(θ) (defined as in (4.6)), and Θ * a compact neighborhood of θ * . Proposition 6 Assume that θ * is the unique stationary point of D(θ) in Θ * and let (θ (n) ) n 0 be a sequence generated by Algorithm 1. If the (θ (n) ) n 0 remain within Θ * , then

lim n→∞ θ (n) = θ * .
The proof is given in Appendix B.1.

For the cases θ = (Σ) or θ = (µ, Σ), the CIM coincides with the FIM (recall Table 3.5). However, it is well know that, near the true value θ * , the Hessian of the function D(θ) is approximated by the FIM [START_REF] Amari | Methods of Information Geometry[END_REF]. Therefore, in these two cases, one should expect the θ (n) to converge to θ * with a superlinear rate of convergence, just like the Newton method does (Theorem 6.3.2 in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]).

Precisely, if θ = (Σ) or θ = (µ, Σ), with a fixed shape parameter β * , then, under the assumptions of Proposition 6, one should expect Algorithm 1 to generate a sequence (θ (n) ) n 0 converging superlinearly to θ * . Here, this will not be proved mathematically, but will be observed experimentally in Subsection 7.2.1 (see Figure 7.1).

CIG Online

The CIG online method is a stochastic gradient method, entirely based on the geometry of the component-wise information metric (CIM), described in Section 3.3. This method aims to minimise the original cost function in Equation (3.1).

In the CIG online method, for the current estimate θ

(n) = (µ (n) , Σ (n) , β (n)
), its correspond stochastic information gradients are given in (3.21). Accordingly, the expected direction of descent is is equal to 0 at any stationary point θ * of D(θ) in (3.1).

As in the classic stochastic gradient descent method, the step-size η (n) = a n is strictly positive, decreasing, and verifies the usual conditions (4.1b).

The choice of the coefficient a requires some attention. Whenever the true shape parameter β * is known, optimal performance is obtained by taking a = 1 (see also Corollary 1, below). However, if β * is unknown, then a has to be chosen manually, by trial and error. This will be discussed in connection with Proposition 8 and with Figure 7.2 in Section 7.2.1.

The CIG online algorithm can now be stated as follows.

Algorithm 2 CIG online algorithm for ECD estimation

Input: A dataset X , an initialization θ (0) = (µ (0) , Σ (0) , β (0) ), the coefficient a > 0; Output: The estimate θ;

1: for n = 0, 1, 2, • • • , N do 2:
η (n+1) ← a n+1 ;

3:

Update µ:

4:

θ current ← (µ (n) , Σ (n) , β (n) ); 5:
Pick a new mini-batch X new mb ;

6:

Compute ∇ * µ p (θ current ; X new mb );

7:

µ (n+1) ← Exp µ (n) η (n+1) ∇ * µ p (θ current ; X new mb ) ; 8:
Update Σ:

9: θ current ← (µ (n+1) , Σ (n) , β (n) ); 10:
Pick a new mini-batch X new mb ;

11:

Compute ∇ * Σ p (θ current ; X new mb );

12:

Σ (n+1) ← Exp Σ (n) η (n+1) ∇ * Σ p (θ current ; X new mb ) ; 13:
Update β:

14:

θ current ← (µ (n+1) , Σ (n+1) , β (n) ); 15:
Pick a new mini-batch X new mb ;

16:

Compute ∇ * β p (θ current ; X new mb );

17:

β (n+1) ← Exp β (n) η (n+1) ∇ * β p (θ current ; X new mb ) ; 18: end for 19: θ ← (µ N +1 , Σ N +1 , β N +1 );
Here, the three exponential maps Exp µ , Exp Σ and Exp β are defined in Equation (3.18), and the stochastic gradients are computed based on (3.21). In each update, X new mb could be just one sample, or a mini-batch. If X new mb is a mini-batch (with more than one sample), then an empirical average over this mini-batch should be introduced, for the stochastic gradients.

The following proposition 7 states the convergence of Algorithm 2. In this proposition, θ * is a stationary point of D(θ) in (3.1), and Θ * is a compact neighborhood of θ * , which verify the following assumption.

A1. θ * is the unique stationary point of D(θ) in Θ * . Moreover, the second derivatives

∇ * 2 µ D(θ * ), ∇ * 2 Σ D(θ * ), ∇ * 2 β D(θ *
) are all positive-definite.

Proposition 7 Assume that A1 is verified and that the estimates (θ (n) ) n 0 generated by Algorithm 2 remain within Θ * . Then, lim θ (n) = θ * almost surely.

The proof of this convergence is given in Appendix B.2.

In addition to almost-sure convergence, stated in Proposition 7, we have also studied the convergence rate and asymptotic normality of Algorithm 2. Recall Θ * is a neighborhood of θ * which satisfies the conditions in Proposition 7. This neighborhood Θ * admits a system of normal coordinates (θ i ; i = 1, • • • , d) with origin at θ * , where d is the dimension of the parameter space Θ, d = m(m+1) 2 + m + 1. Since D(θ) has an isolated stationary point at θ = θ * , its Hessian at this point θ = θ * can be expressed in normal coordinates

H ij = ∂ 2 D ∂θ i ∂θ j θ i =0 (5.3)
When θ * is a local minimum of D(θ), the matrix H = (H ij ) is positive-definite [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. With these notations, the rate of convergence is given by the following Proposition 8.

Recall that the step-sizes in Algorithm 2 are of the form η (n) = a/n, where the constant a is chosen by the user.

Proposition 8 Under the assumptions of Proposition 7, if a > 1 2λ , where λ > 0 is the smallest eigenvalue of H,

E d 2 (θ * , θ (n) ) = O(n -1 ) (5.4)
Here, d (•, •) stands for the CID in (3.19), and the "big O" notation means that there exist K > 0 and n 0 > 0 such that

∀n n 0 E[d 2 (θ * , θ (n) )] K n
In terms of the normal coordinates (θ i ), let the component-wise information gradient ∇ θ p (θ * ; x) at the point θ = θ * , given by (3.21), have components

(u i (θ * )). Let G * = (G * ij ), be the matrix G * ij = E θ * u i (θ * )u j (θ * ) (5.5)
Then, the following proposition gives the asymptotic normality of the CIG online algorithm.

Proposition 9 (asymptotic normality) Under the assumptions of Propositions 7 and 8, the distribution of the re-scaled coordinates (n

1 2 θ i ) i∈{1,••• ,d} converges to a centred d-variate normal distribution
, where d is the dimension of Θ, with covariance matrix G given by the following Lyabunov equation

AG + GA = -a 2 G * (5.6)
Here, A = (A ij ) with A ij = 1 2 δ ij -aH ij (δ denotes Kronecker's delta). The proofs of Propositions 8 and 9 and are discussed in Appendix B.3.

For the cases θ = (Σ) or θ = (µ, Σ), the CIM coincides with the FIM (recall Table 3.5). Therefore the component-wise information distance (the CID (3.19)) also coincides with the "true" information distance, which is the distance associated with the FIM. Thus, in these cases, Assumptions (d1) to (u1) in Proposition 5 of the previous Chapter 4 are satisfied, and the following corollary may be obtained [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF].

Corollary 1 For the ECD model, parameterised by θ = (Σ) or θ = (µ, Σ), with a fixed β * , the CIM (3.16) coincides with the Fisher information metric.

1. the rate in equation (5.4) holds, whenever a > 1/2.

2. if a = 1 the distribution of the re-scaled coordinates (n 1/2 θ i ) converges to a centred d-variate normal distribution, with covariance matrix equal to the identity G * = I d , and the recursive estimates θ (n) are asymptotically efficient.

Note that, Item 2) of Corollary 1 implies that the distribution of nd 2 * (θ * , θ (n) ) converges to a χ 2 -distribution with d degrees of freedom.

nd 2 * (θ * , θ (n) ) dist -→ χ 2 m(m + 1) 2 for θ = (Σ) nd 2 * (θ * , θ (n) ) dist -→ χ 2 m(m + 1) 2 + m for θ = (µ, Σ)
This provides a practical means of confirming the asymptotic normality of the estimators θ (n) .

Global convergence analysis

This section studies the global convergence of the CIG algorithm, in both its offline and online versions. The term "global convergence" means that the iterates (θ (n) ) n≥0 , generated by Algorithm 1 or Algorithm 2, always converge to the true parameter value θ * , whatever the initial guess θ (0) . We will mainly consider the two most well-known sub-families of ECD, the Multivariate Generalized Gaussian Distribution (MGGD) and Multivariate Student t-Distribution, and only in the cases θ = (Σ) and θ = (µ, Σ). The main results are stated in the following two tables. In these tables, the shape parameter β is considered fixed and known.

MGGD θ = (Σ)

global convergence holds for β > 0 θ = (µ, Σ) global convergence holds for β > 0 for the online version) has a unique stationary point at the true parameter θ * , which is the global minimizer. This follows from the strict g-convexity (i.e. geodesic convexity of this cost function), which will now be considered.

First, for the case of θ = (Σ) with known µ * and β * , let

f(δ, β) = 1 g β (δ) (5.7)
then, for the MGGD model

f(δ, β) = exp 1 2 δ β β > 0 (5.8a)
and for the Student t model,

f(δ, β) = 1 + δ β β+m 2 β > 0 (5.8b)
The following proposition introduces a sufficient condition for the Kullback-Leibler divergence D(Σ) and its empirical approximation D(Σ) to be geodesically strictly convex.

Proposition 10 Assume that the function f : R + → R + in (5.7) verifies the following condition : for any function ϕ : R → R + ϕ is strictly log-convex ⇒ f • ϕ is strictly log-convex (5.9)

Then, the Kullback-Leibler divergence D(Σ) and the empirical divergence D(Σ) are geodesically strictly convex.

In particular, when the conclusion of Proposition 10 holds, the unique global minimum of D(Σ) is at the true Σ * . Moreover, Σ * is also the unique stationary point of D(Σ). Proposition 10 directly yields the following corollary, for the MGGD and Student T models. By plugging (5.8a) and (5.8b) into (5.9), it is easy to obtain the following. Compared with [START_REF] Zhang | Multivariate generalized Gaussian distribution: Convexity and graphical models[END_REF], this proposition can give a wider range of shape (or degree of freedom) parameters, for geodesic convexity.

Corollary 2 the Kullback-Leibler divergence D(Σ) and its empirical approximation D(θ) are geodesically strictly convex, with unique global minimum (and unique stationary point), in both of the following cases.

1. the dataset X is distributed according to an MGGD model, with scatter matrix Σ * and with shape parameter β > 0.

2. the dataset X is distributed according to a Student T model, with scatter matrix Σ * and degree of freedom β > 0.

Thus, when Σ is unknown and β satisfies the conditions of Corollary 2, this corollary implies the global convergence of Algorithms 1 and 2. Precisely, these algorithms will always converge to the true value θ * = (Σ * ) of the parameter θ = (Σ), whatever the initial guess θ (0) . For the more complicated situation θ = (µ, Σ), global convergence does not always hold. The cost function D(θ) is not geodesically convex, but may be reformulated, using a new matrix argument [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF].

S = Σ + µµ † µ µ † 1
Precisely, if the new random vector y is given by

y † = x † , 1 †
then the cost function can be reformulated as

D(θ) = - 1 2 log det(S) -log f(δ y )
where

δ y = y † S -1 y = (x -µ) † Σ -1 (x -µ) + 1
Here, the function f for MGGD is

f(δ y ) = exp 1 2 (δ y -1) β
and for Student T-distributions

f(δ y ) = 1 - 1 β + δ y β β+m 2
(5.10)

In [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF], the minimization of D(θ) was proven to be equivalent to the minimization of D(θ).

Replacing the new function f into (5.9), the following corollary is obtained.

Corollary 3

The KL divergence D(µ, Σ), and the empirical divergence D(µ, Σ), both have a unique global minimum (and unique stationary point) at (µ * , Σ * ), in the following cases 1. the dataset X is distributed according to an MGGD model, with expectation and scatter matrix (µ * , Σ * ) and with fixed shape parameter β > 0.

2. the dataset X is distributed according to a Student t model, with expectation and scatter matrix (µ * , Σ * ) and with the fixed degree of freedom β > 0.

For these two cases, global convergence is then guaranteed. Finally, for the most complicated case θ = (µ, Σ, β), the cost function is always nonconvex. Moreover, we have verified experimentally that it has multiple stationary points in the parameter space Θ = R m × P m × R + . Therefore, the iterates (θ (n) ) n≥0 generated by Algorithms 1 and 2 cannot be guaranteed to converge to the true parameter θ * , unless the initial guess θ (0) is close enough to θ * .

Conclusion

This chapter introduced the CIG methods, which has two versions, an offline (batch) version and an online (stochastic) version. Both versions aim to estimate the complete ECD model, where all the parameters are unknown, including the shape parameter β.

Propositions 6 and 7 state the convergence of the offline and online CIG methods, respectively. For Proposition 7, the assumption A1 is required.

The online CIG method retains some of the good properties of the Riemannian information gradient method, seen in the previous Chapter 4. For example, Propositions 8 and 9 state the fast convergence and asymptotic normality of this method, provided the step-size η (n+1) = a n+1 is chosen correctly, with a > 1 2λ as explain in Proposition 8. In general, when the shape parameter β is unknown, this choice of step-size is carriedout manually, and may be very time consuming. On the other hand, when β is known, Corollary 1, after Proposition 9, shows that it is possible to recover all of the good properties of the Riemannian information gradient method, by choosing η (n+1) = 1 n+1 . In particular, the online CIG method even turns out to be asymptotically efficient.

In fact, when β is known, there is also a notable improvement in the performance of the offline CIG method, which turns out to have a super-linear rate of convergence.

For the MGGD and Student T models, the global convergence of the CIG methods was considered in Section 5.4, and it summarized in Tables 5.1 and 5.2. For the cases indicated in these tables, the CIG method converges to the true parameter θ * = (µ * , Σ * ) independently of the initial guess θ (0) .

All of the results obtained in this chapter will be illustrated by computor experiments in Chapter 7. The following Chapter 6 will extend the online CIG method to the estimation of mixtures of ECD (MECD).

Chapter 6

Online estimation of MECD

Introduction

The main body of this chapter is based on our article [START_REF] Zhou | Online estimation of mixtures of ECD: Component-wise information gradient method[END_REF], which is currently under review.

An MECD model is a mixture (i.e. a weighted sum) of K ECDs, where K is called the number of components. The dimension of its parameter space can be quite large, as it is more than K times the dimension of the parameter space of a single ECD model. In addition, MECD are typically used to model more complex, large-scale datasets. Accordingly, traditional offline methods may become impractical when applied to the estimation of MECD, due to lack of computational resources (time and memory).

Therefore, it seems more suitable to use online (stochastic) methods. However, Euclidean stochastic gradient methods may need increasingly large mini-batch sizes, to overcome instability [START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF]. Besides, choosing optimal step-sizes, for stochastic gradient methods, is a difficult task, especially for the estimation of MECD.

The present chapter proposes a different approach, by applying the CIG online method to the estimation of MECD. This extends the previous Chapter 5, which focused on estimating a single ECD, rather than a full mixture.

The geometric background, needed to introduce the CIG online method for MECD estimation, was given in Section 3.4 of Chapter 3. This includes the component-wise information metric (CIM) and component-wise information gradient (CIG) for MECD, respectively given by (3.23) and (3.28).

The CIG online method for the estimation of MECD has two versions, a decreasing step-size version (CIG-DS) and an adaptive step-size version (CIG-AS). The decreasing step-size version is a direct extension of the CIG online method of the previous chapter, Algorithm 2. It can be obtained by introducing the retraction map (3.27) and CIG (3.28) of MECD models, instead of the retraction map (3.18) and CIG (3.21) of ECD models.

This decreasing step-size version requires manual selection of the step-sizes, according to a rule similar to the one in Proposition 8 of the previous Chapter 5, for single ECD. This is not always practical for MECD.

The adaptive step-size version involves an automatic selection of the step-sizes, which guarantees a fast rate of convergence. Automatic selection of the step-sizes is very helpful in saving time, since a manual selection involves running the algorithm several times, selecting a suitable step-size by trial and error.

The CIG method turns out to have several significant advantages. First, being an online (i.e. stochastic) method, each iteration requires access to only one mini-batch, of a constant size, instead of the whole dataset (or mini-batch with increasing size). This considerably reduces the required time and memory usage. Second, when shape parameters are unknown, it still converges to accurate estimates, without increasing the size of mini-batches (compare to [START_REF] Hertrich | Inertial stochastic PALM and its application for learning student-t mixture models[END_REF]). Third, the CIG-AS variant (with an adaptive step-size) avoids manual step-size selection.

In the following, Section 6.2 presents the CIG-DS method, as well as a short description of another method, which is an online backtracking method, called CIG-OB. In Section 6.3, the CIG-AS method is stated in detail. Here, the theoretical convergence properties, of the CIG-DS and CIG-AS methods, will be stated in Propositions 11 and 12, respectively. In the following Chapter 7, these properties will be evaluated by numerical experiments, in Section 7.2.2.

CIG-DS

The CIG-DS method uses the most common decreasing step-size, η (n) = a n with a > 0, in order to verify the usual stochastic approximation Conditions (4.1b). The constant a should be selected in advance, before running the algorithm, (just as in Algorithm 2). The choice of this constant has a crucial influence on the rate of convergence.

Recall here the MECD density (1.2) is parameterised by θ = (r, (µ k ) k , (Σ k ) k , (β k ) k ) where k = 1, . . . , K, and r = (r k ) k is related to the mixture weights w = (w k ) k by w k = r 2 k .

Algorithm 3 CIG-DS algorithm for MECD

Input: A dataset X , an initialization θ (0) = r (0) , (µ

(0) k ) k , (Σ (0) k ) k , (β (0)
k ) k , a constant a > 0; Output: The estimate θ;

1: for n = 0, 1, 2, • • • , N do 2: η (n+1) ← a n+1 ; 3:
Update r:

4: θ current ← r (n) , (µ (n) k ) k , (Σ (n) k ) k , (β (n) k ) k ; 5:
Pick a new mini-batch X new mb ;

6:

Compute ∇ * r f (θ current ; X new mb ); 7: r (n+1) ← Exp r (n) η (n+1) ∇ * r f (θ current ; X new mb ) ; 8:
Update µ:

9: θ current ← r (n+1) , (µ (n) k ) k , (Σ (n) k ) k , (β (n) k ) k ; 10: for k = 1, 2, • • • , K do 11:
Pick a new mini-batch X new mb ;

12:

Compute ∇ * µ k f (θ current ; X new mb ); 13: µ (n+1) k ← Exp µ (n) k η (n+1) ∇ * µ k f (θ current ; X new mb ) ; 14:
end for 15:

Update Σ:

16:

θ current ← r (n+1) , (µ (n+1) k ) k , (Σ (n) k ) k , (β (n) k ) k ; 17: for k = 1, 2, • • • , K do 18:
Pick a new mini-batch X new mb ;

19:

Compute ∇ * Σ k f (θ current ; X new mb ); 20: Σ (n+1) k ← Exp Σ (n) k η (n+1) ∇ * Σ k f (θ current ; X new mb ) ; 21:
end for

22:

Update β:

23:

θ current ← r (n+1) , (µ (n+1) k ) k , (Σ (n+1) k ) k , (β (n)
k ) k ;

24:

for k = 1, 2, • • • , K do 25:
Pick a new mini-batch X new mb ;

26:

Compute ∇ * β k f (θ current ; X new mb ); 27: β (n+1) k ← Exp β (n) k η (n+1) ∇ * β k f (θ current ; X new mb ) ; 28:
end for 29: end for 30: θ ← r (N +1) , (µ

(N +1) k ) k , (Σ (N +1) k ) k , (β (N +1) k ) k
In the above algorithm, all the exponential maps are given in (3.27), and all the gradients are computed using (3.28) and the formulas just after.

The following proposition gives the convergence and rate of convergence for Algorithm 3. In this proposition, D(θ) is the Kullback-Leibler divergence (3.3).

Proposition 11 (Convergence rate of CIG-DS) Assume the Kullback-Leibler divergence function D(θ) has a stationary point at θ = θ * , Θ * ⊂ Θ is a compact and convex neighborhood of θ * , such that θ * is the unique stationary point of D(θ) in Θ * , and the estimates (θ (n) ) n 0 , generated by Algorithm 3, remain within Θ * . Then, for the sequence

(θ (n) ) n 0 , lim n→∞ θ (n) = θ * If a > 1 2λ
, where λ > 0 is the smallest eigenvalue of the Hessian matrix of D(θ) at the point θ = θ * , then

D(θ (n) ) -D(θ * ) = O(n -1 ) (6.1)
Here, the "big O" notation means that there exist C > 0 and n 0 > 0 such that

∀n n 0 D(θ (n) ) -D(θ * ) C n
The proof of this proposition is given in Appendix C.1. In Algorithm 3, X new mb denotes a mini-batch of samples drawn at random from the complete dataset X . Theoretically, the optimal value for a can be selected based on the condition a > 1 2λ . In practice, it is very difficult to do so, and a blindly selected a may not guarantee the convergence rate in equation (6.1). One way of overcoming this issue is to introduce an online backtracking line-search technique [START_REF] Vaswani | Painless stochastic gradient: Interpolation, line-search, and convergence rates[END_REF]. When applied to the CIG method, this leads to the following CIG-OB (online backtracking) version.

Algorithm 4 CIG-OB algorithm for MECD

Input: A dataset X , an initialization θ (0) = r (0) , (µ (0) k ) k , (Σ (0) k ) k , (β (0) k ) k ; Output: The estimate θ; 1: for n = 0, 1, 2 • • • , N do 2:
Update r:

3: θ current ← r (n) , (µ (n) k ) k , (Σ (n) k ) k , (β (n) k ) k ; 4:
Pick a new mini-batch X new mb ;

5:

Compute ∇ * r f (θ current ; X new mb );

6: η r ← Armijo-Goldstein(X new mb , θ current ); 7: r (n+1) ← Exp r (n) (η r ∇ * r f (θ current ; X new mb )); 8:
Update µ:

9: θ current ← r (n+1) , (µ (n) k ) k , (Σ (n) k ) k , (β (n) k ) k ; 10: for k = 1, 2, • • • , K do 11:
Pick a new mini-batch X new mb ;

12:

Compute ∇ * µ k f (θ current ; X new mb );

13:

η µ k ← Armijo-Goldstein(X new mb , θ current ); 14: µ (n+1) k ← Exp µ (n) k (η µ k ∇ µ k f (θ current ; X new mb )); 15:
end for 16:

Update Σ:

17:

θ current ← r (n+1) , (µ (n+1) k ) k , (Σ (n) k ) k , (β (n) k ) k ; 18:
for k = 1, 2, • • • , K do 19:
Pick a new mini-batch X new mb ;

20:

Compute ∇ * Σ k f (θ current ; X new mb ); 21: η Σ k ← Armijo-Goldstein(X new mb , θ current ); 22: Σ (n+1) k ← Exp Σ (n) k (η Σ k ∇ Σ k f (θ current ; X new mb )); 23:
end for

24:

Update β:

25:

θ current ← r (n+1) , (µ (n+1) k ) k , (Σ (n+1) k ) k , (β (n) k ) k ; 26:
for k = 1, 2, • • • , K do 27:
Pick a new mini-batch X new mb ;

28:

Compute ∇ * β k f (θ current ; X new mb ); 29: η β k ← Armijo-Goldstein(X new mb , θ current ); 30: β (n+1) k ← Exp β (n) k (η β k ∇ β k f (θ current ; X new mb )); 31:
end for 32: end for 33: θ ← r (N +1) , (µ

(N +1) k ) k , (Σ (N +1) k ) k , (β (N +1) k ) k
In the algorithm above, Armijo-Goldstein(•, •) means the Armijo-Goldstein backtracking criterion [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF][START_REF] Vaswani | Painless stochastic gradient: Interpolation, line-search, and convergence rates[END_REF], for which one needs to provide the current estimator and a mini-batch. The CIG-OB method inherits the properties of classic deterministic line-search methods. With respect to number of iterations, it has a fast rate of convergence, but the actual time consumption is more expensive than CIG-DS. Moreover, the size of the mini-batch significantly affects its accuracy, this will be discussed in detail, in Section 7.2.2.

CIG-AS

The CIG-AS method is similar in structure to the CIG-DS method. The main difference between the two methods is the following. In CIG-DS, the step-size η (n) is chosen in advance (for example, η (n) = a n ). On the other hand, in CIG-AS, η (n) is computed based on the current estimate θ (n) , in order to ensure a faster rate of convergence. In other words, CIG-AS involves an adaptive choice of the step-size η (n) .

The idea of introducing an adaptive step-size is loosely based on [START_REF] Yi Meng | Fast and furious convergence: Stochastic second order methods under interpolation[END_REF], which deals with classical (Euclidean) gradient descent.

All the steps of CIG-AS are similar to algorithm 3, except the step-size is adaptive. This is reflected in the following algorithm (in this algorithm, the adaptive step-size η (n) is set in line 2).

Algorithm 5 CIG-AS algorithm for MECD

Input: A dataset X , an initialization θ (0) = r (0) , (µ (0) k ) k , (Σ (0) k ) k , (β (0) k ) k ; Output: The estimate θ; 1: for n = 0, 1, 2 • • • do 2: η (n+1) ← τ (n) min ρ L τ (n) max ; 3:
Update r:

4: θ current ← r (n) , (µ (n) k ) k , (Σ (n) k ) k , (β (n) k ) k ; 5:
Pick a new mini-batch X new mb ;

6:

Compute ∇ * r f (θ current ; X new mb ); 7: r (n+1) ← Exp r (n) η (n+1) ∇ * r f (θ current ; X new mb ) ; 8:
Update µ:

9:
θ current ← r (n+1) , (µ

(n) k ) k , (Σ (n) k ) k , (β (n) k ) k ; 10: for k = 1, 2, • • • , K do 11:
Pick a new mini-batch X new mb ;

12:

Compute ∇ * µ k f (θ current ; X new mb ); 13: µ (n+1) k ← Exp µ (n) k η (n+1) ∇ * µ k f (θ current ; X new mb ) ; 14:
end for 15:

Update Σ:

16:

θ current ← r (n+1) , (µ (n+1) k ) k , (Σ (n) k ) k , (β (n) k ) k ; 17: for k = 1, 2, • • • , K do 18:
Pick a new mini-batch X new mb ;

19:

Compute ∇ * Σ k f (θ current ; X new mb ); 20: Σ (n+1) k ← Exp Σ (n) k η (n+1) ∇ * Σ k f (θ current ; X new mb ) ; 21:
end for

22:

Update β:

23:

θ current ← r (n+1) , (µ (n+1) k ) k , (Σ (n+1) k ) k , (β (n)
k ) k ;

24:

for k = 1, 2, • • • , K do 25:
Pick a new mini-batch X new mb ;

26:

Compute ∇ * β k f (θ current ; X new mb ); 27: β (n+1) k ← Exp β (n) k η (n+1) ∇ * β k f (θ current ; X new mb ) ; 28:
end for 29: end for 30: θ ← r (N +1) , (µ

(N +1) k ) k , (Σ (N +1) k ) k , (β (N +1) k ) k
In the algorithm above,

τ (n) min = min 1, λ (n) min,k I (n) µ k k , J (n) Σ k ,2 k , 1
I (n) β k k τ (n) max = max    1, λ (n) max,k I (n) µ k 2 k , J (n) Σ k ,1 2 k , 1
I (n) β k 2 k  
 where λ min,k is the smallest eigenvalue of Σ k , and λ max,k is the biggest eigenvalue of Σ k . Moreover, the coefficients I µ k , J Σ k ,1 , J Σ k ,2 , and I β k are computed from (2.10) and (3.25).

The following proposition gives the convergence of this algorithm. The meaning of strong geodesic convexity, for Condition (i), was explained in Definition 1, of Chapter 2. The definitions of Conditions (ii) and (iii) can be found in Equation (C.7) and Equation (C.20), of Appendix C.2.

The L-geodesically smooth assumption gives a sequence that upper bounds the error E Xnew [D(θ (n+1) )]-D(θ * ). Then, this sequence is proved to be convergent to 0 by employing the α-geodesically strongly convexity and ρ-strong growth condition. The step-size which convient to this convergence condition is under the following form

η (n+1) = τ (n) min ρ L τ (n) max
Proposition 12 (Linear convergence of CIG-AS) Assume that θ * is the unique stationary point of the Kullback-Leibler divergence D(θ) in a compact convex neighborhood Θ * . Assume that D(θ) satisfies the following properties on Θ * , (i) α-geodesically strongly convex [START_REF] Zhang | First-order methods for geodesically convex optimization[END_REF],

(ii) L-geodesically smooth [START_REF] Zhang | First-order methods for geodesically convex optimization[END_REF],

(iii) satisfies the ρ-strong growth condition [START_REF] Vaswani | Fast and faster convergence of SGD for over-parameterized models and an accelerated perceptron[END_REF].

If θ (0) belongs to Θ * , then the sequence generated by Algorithm 5 with the step-size

η (n) = τ (n) min ρ L τ (n) max
and a constant batch size b converges to θ * , with the following rate

E Xmb D θ (N ) -D (θ * ) c(N ) D θ (0) -D (θ * ) (6.2)
where

c(N ) = N n=1      1 - α τ (n) min 2 ρ L τ (n) max      ρ = (N -b)(ρ -1) (N -1)b + 1
Remark that, c is a real positive scalar in (0, 1), not the normalizing constant. The constant b is the size of mini-batch, X mb denotes the law for selection of mini-batch, and the constant N denotes the number of iterations.

The proposition above shows that the convergence rate is linear. Linear convergence means that the following definition holds [START_REF] Senning | Computing and estimating the rate of convergence[END_REF].

Definition 2 Suppose that the sequence D(θ (N ) ) converges to the number D(θ * ). The sequence is said to converge linearly to D(θ * ) if there exist constants c, A > 0 such that

log D(θ (N ) ) -D(θ * ) -cN + A (6.3)
Indeed, it can be shown that this definition holds, as a consequence of Proposition 12, simply by taking logarithms in Equation (6.2), and noting that log c(N ) -cN , where

c = min n 1        log     1 1 - α τ (n) min 2 ρ L τ (n) max            > 0
The fact that c > 0 can be proved using the fact that θ (n) → θ * .

Note that, in practice, it may be difficult to select L and ρ. Too big a value of L and ρ will lead to a very slow convergence at early stages of the estimation. Therefore, in numerical experiments and real applications, we used the following step-size, a combination of the CIG-DS and CIG-AS step-sizes,

η (n) = max τ (n) min ρ L τ (n) max , a n + 1 (6.4)
The use of this "hybrid" step-size selection rule will be illustrated in Subsection 7.2.2 of the following Chapter 7.

Conclusion

The present chapter introduced the CIG method for the estimation of MECD. Since MECD models are high-dimensional and typically used to model large datasets, only the online form of the CIG method was developed. In particular, two new algorithms were proposed, CIG-DS (CIG with decreasing stepsize) and CIG-AS (CIG with adaptive step-size), in Sections 6.2 and 6.3 respectively. The rate of convergence of the CIG-DS method was stated in Proposition 11, and the rate of convergence of the CIG-AS method was stated in Proposition 12.

The CIG-AS method carries out an adaptive step-size selection, avoiding the difficult and time-consuming manual step-size selection. It also achieves a faster, linear rate of convergence (in the sense of Definition 2).

The present chapter completes the theoretical part of this thesis. We have now extended CIG estimation methods to the estimation of complete MECD models. In the following chapter, we will present numerical experiments, with simulated and real data, in order to illustrate the performance of CIG methods, theoretically discussed in Chapters 5 and 6.

Introduction

The present chapter is the final chapter, in the main body of this thesis. It aims to illustrate, through applications to simulated and real data, the various theoretical properties of the CIG estimation methods, for ECD and MECD, as introduced in Chapters 5 and 6.

Section 7.2 is devoted to experiments carried out with simulated data. In this section, Subsection 7.2.1 corresponds to the CIG methods of Chapter 5, applied to the estimation of ECD. Subsection 7.2.2 is devoted to the CIG methods of Chapter 6, applied to the estimation of MECD.

Section 7.3 is concerned with applications to real data. Subsection 7.3.1 applies the CIG offline and online methods, from Chapter 5, to colour transformation for image editing. Subsection 7.3.2 applies the CIG methods of Chapter 6 to texture segmentation.

The results of Section 7.2 confirm the theoretical results in Chapters 5 and 6. In particular, they verify the conclusions of Proposition 8, Corollary 1, and Propositions 11 and 12.

In addition, Section 7.2 demonstrates some properties which were not obtained theoretically in Chapters 5 and 6, such as the super-linear convergence rate of the CIG offline method. It also compares CIG methods to other, state-of-the-art methods, such as method of moments, fixed-point, expectation-maximization, and classical stochastic gradient [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF][START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF]. The comparison involves both estimation accuracy and computation time.

The applications of Section 7.3 showcase all the advantages of CIG online methods, when used on real data. For example (see Subsection 7.3.1), the CIG online method is able to perform the crucial estimation step, for a colour transformation application to a pair of images, with over 10 6 pixels, in only 21 seconds. For this this application, other state-of-the-art methods require two hours to do the same job, with a comparable performance. The first experiment confirms the super-linear convergence rate of the CIG Offline method, for a dataset, distributed according to the MGGD model, which contains N = 10 4 samples. The initial value θ (0) is defined as the MM estimate, using 10% of the entire dataset. Figure 7.1(a) presents the case of θ = (Σ) with known (µ * , β * ). The CIG Offline method converges after only two iterations, and if the same accuracy needs to be achieved, the classic deterministic gradient method (not using the information gradient) requires at least 88 iterations. For the case of θ = (µ, Σ) with known (β * ), things are similar. Figure 7.1(b) shows that CIG Offline method, after two iterations, achieves the same accuracy as the traditional deterministic gradient method, after 200 iterations. Here, the traditional deterministic gradient means the classical (Euclidean) deterministic gradient method, with step-sizes also selected according to the Armijo-Goldstein line search criteria [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF].

Computer experiments

In Figure 7.1, d 2 * (θ * , θ (n) ) denotes the information distance, given by Equation (3.13). In Figure 7.2, d 2 (θ * , θ (n) ) denotes the component-wise information distance, given by Equation (3.19).

The second experiment demonstrates the convergence rate of the CIG online method. In this experiment, both MGGD and Student T datasets are used. The initialization θ (0) is randomly chosen. ) where N = 10 5 . These "estimated p.d.f." coincide very closely with the theoretical chi-squared probability density, confirming the fact that θ (N ) is asymptotically normally distributed about θ * .

In the third experiment, we compare the efficiency of the CIG offline and online methods with other common estimation methods, MM (method of moments) and FP (fixedpoint method). In each trial, the dataset is generated from an MGGD model. For MGGD, the MM was given in [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF], and the FP method in [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF]. In The second additional experiment shows that, for the CIG online method, even tough, the gradient of the Kullback-Leibler divergence converge to 0 (as seen in the first additional experiment), the iterates θ (n) do not necessarily converge to the global minimum θ * (the true parameter value). In Figure 7.5(b), two different initial values were used for this method. The blue curve has θ (0) farther away from θ * . This curve shows that the component-wise information distance d (θ * , θ (n) ) converges to a non-zero constant, meaning that θ (n) do not converge to θ * . The red curve has θ (0) closer to θ * . In this case, the component-wise information distance d (θ * , θ (n) ) converges to 0. Of course, this means θ (n) converge to θ * .

In conclusion, for θ = (µ, Σ, β), the convergence to the global minimum θ * can only be guaranteed locally. For the CIG offline methods, if the initial guess θ (0) is chosen in a neighbourhood Θ * of θ * , the θ (n) always converge to θ * , as soon as Θ * satisfies the assumptions of Proposition 6.

The difference between the CIG online and offline methods is that, due to its stochastic nature, CIG online may jump out of the neighbourhood Θ * during the first few iterations. This leads to convergence to a local minimum, different from θ * . This explains why the final averaged accuracy of CIG online is not as good as that of the CIG offline and FP methods, making the variance of the CIG online estimator notably larger.

As a possible remedy to this problem, the mini-batch CIG online was also tested, and compared with other methods, in the Figure 7.4(c). Two sizes of the mini-batch, 10 and 100 samples, were considered. However, the results show that increasing the mini-batch size has no significant effect on the accuracy of CIG online. As for computation time, CIG online method has a significant advantage. The computation time of the CIG online algorithm is similar to that of MM, and is significantly less than that of FP. Meanwhile, its accuracy is significantly better than that of MM. In most experiments, the accuracy of CIG online method is similar to, or even better than, FP. Although the computation time of CIG offline is greater than that of CIG online, it is comparable to that of FP, while, in most cases, CIG offline can achieve the best accuracy, among the four estimation methods considered.

Mixture of ECD

In this subsection, computer simulations are presented, to evaluate the performance of the CIG online method for MECD estimation, as considered in Chapter 6. These simulations will focus on two models, mixtures of MGGD and mixtures of Student T-distriutions.

A comparison of the CIG online method with the already existing EMFP (Expectation-Maximisation with Fixed-Point [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Hasannasab | Alternatives to the EM algorithm for ML estimation of location, scatter matrix, and degree of freedom of the Student t distribution[END_REF]) and SG (Euclidean Stochastic Gradient [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF]) methods is provided.

The dimension of the observed samples (datapoints) is supposed to be m = 5 and the number of mixture components K = 5. The mixture weights w = r 2 are drawn from a uniform distribution on the unit sphere (recall this representation of mixture weights, in terms of the unit sphere, from Subsection 3.4). The true location parameters (µ * k ) k are randomly chosen from a multivariate normal distribution. The true scatter matrices (Σ * k ) k are defined as Σ(i, j) = ρ |i-j| for i, j ∈ {1, • • • , m}, and ρ uniformly distributed in the interval (0.2, 0.8). Then, the shape parameters (degree of freedom for Student T) (β * k ) k are uniformly chosen in the interval [0.2, 10].

For each scenario, we have done 100 Monte Carlo trials. The dataset for each trial has 3 × 10 5 observations. For all scenarios, algorithms are initialized by a k-mean with 0.1% of the complete dataset.

The first point we would like to illustrate is the rate of success. Under the same initialization conditions, we compared the rate of success of five methods (CIG-DS, CIG-OB, CIG-AS, EMFP, and Euclidean SG). Here, 'success' means that the log-likelihood increases with each update and finally converges.

In [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF], expectation-maximisation using the fixed-point algorithm in [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF] (EMFP) is proved to be the most robust method for mixtures of MGGD. Therefore, this method is included in our comparison. Also in [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF], an online classic (that is, Euclidean) SG method was introduced, which we also included for comparison. All gradient methods (CIG and SG) used the same mini-batch size b = 1000.

In the following Table 7.1, the rates of success are presented. We can see that the rates of success of CIG methods are significantly better. This shows that the intrinsic geometric structures on the parameter space provide a more stable estimation process. Second, we confirm the convergence rates of the CIG-DS and CIG-AS methods, stated in Propositions 11 and 12. In Figure 7.7, two log-log plots are presented. In these plots, the x-axis denotes the number of updates and the y-axis denotes the empirically estimated Kullback-Leibler divergence (4.6). In both these plots, the slope of the curve approaches -1, which means the empirical divergence D(θ) is O(n -1 ), as in (6.1).

As mentioned in Proposition 11, the coefficient a, used in computing step-sizes (η (n) = a/n), needs to verify the condition a > 1 2λ . Here, since it is quite difficult to compute λ, the coefficient a was selected manually, through trial and error. This procedure lead to a = 1000 for mixtures of MGGD, and a = 10 for mixtures of Student T.

Then, we compared the performance of the five methods (CIG-DS, CIG-OB, CIG-AS, EMFP, and Euclidean SG). In Figure 7.8, the x-axis represents the number of epochs (one epoch means all samples in the dataset have been traversed exactly once), and the y-axis represents the log-likelihood. Finally, we observed the variation of the log-likelihood versus time consumption. In Figure 7.9, the x-axis represents the time consumption (in seconds), and the y-axis represents the variation of log-likelihood.

For mixtures of MGGD, EMFP converged in 1000 seconds. The classic (Euclidean) SG is similar to EMFP for time consumption. Then, CIG-DS and CIG-AS are the fastest methods, as they converged in 600 seconds, but CIG-OB did not show any advantages in terms of time consumption. We think that the reason is the sensitivity for β in the MGGD case, which requires more time for the best step-size to be selected.

For mixtures of Student T, EMFP took 1750 seconds to converge to the maximum likelihood value, while the SG did not converge in this time interval. The CIG-DS and CIG-AS methods showed the fastest convergence, as they converged in 750 seconds. CIG-OB converged in less than 250 seconds, but it is still not rather inefficient. Based on the Here, we consider a first application of the CIG method for ECD estimation, to colour transformation for image editing [START_REF] Hristova | Transformation of the multivariate generalized Gaussian distribution for image editing[END_REF]. Precisely, this application will focus on the estimation of MGGD. Its goal is to replace the colour distribution of an input image by the colour distribution of a given target image. The main idea is to fit the input and the target distributions, with two different MGGD models. Then, the transformation between these two MGGD is implemented by a linear Monge-Kantorovich transformation for Σ, and a stochastic transformation for β [START_REF] Hristova | Transformation of the multivariate generalized Gaussian distribution for image editing[END_REF][START_REF] Pitié | The linear monge-kantorovitch linear colour mapping for example-based colour transfer[END_REF]. In each image, the pixels correspond to 3-dimensional vectors, in the case of RGB images, and to 5-dimensional vectors, when spatial gradient-field information is added to colour information.

Starting with the 3D RGB case, Figure 7.10 presents the transformed images and some of their details. The detail (a1) clearly shows that the cloud 'drawn' by MM appears too green. Similarly, FP also presents a green appearance, in detail (a2). On the contrary, the CIG offline and CIG online methods show pure white cloud colour in (a3) and (a4). Note also the difference in the amount of blue in the shadows on the grass. Too much blue is mixed with the shadows, in MM's output detail (b1). In details (b2),(b3),(b4), the results of FP and CIG methods lead to a more natural appearance.

From our point of view, the most interesting aspect of this application is in terms of computation time. The CIG online method takes about 10 seconds for two images (input and output). In contrast, FP and CIG offline methods each require more than two hours. In other words, CIG online has a decisive advantage, in terms of time consumption.

Then, gradient-field information was included, so the transformation came to involve 5-dimensional vectors, which consist in three colour components (of CIELAB) and two components of the image spatial gradient field (dx and dy). For this application, the shape parameter of the MGGD model was supposed to fixed. Figure 7.11 presents the four different outputs. It can be observed that the output of the FP and CIG methods is significantly better than that of MM. In the transformed result of MM, the hue is darker and greener. FP and CIG results are better, since the frost on the grass is whiter and appears more natural, and the forest on the mountain in the image also appears darker. The two images in Figure 7.11 have more than 1.2 × 10 6 pixels (i.e. 1.2 × 10 6 samples). The FP and CIG offline methods need more than 4 hours to run, on the these two images. The CIG online method needs only 21 seconds.

We also considered an application to full HD images. In this case, as demonstrated in Figure 7.12, the advantages of the CIG online algorithm were significant. The result of MM failed to achieve the colour of the autumn leads in the target image, showing light green instead of yellow. Since the input image and the target image have more than 4 × 10 6 pixels (that is 4 × 10 6 samples), it was not feasible to run FP and CIG offline, with the entire dataset. Rather, the estimation was done on subsets of the complete dataset. These two subsets have 4×10 5 samples, that are randomly taken from the original images. In the autumn leafs obtained using FP and CIG offline, the yellow colour has obviously been smeared. CIG online is more natural, and the yellow colour is more uniform, and it is closer to the style of the target image. This section applies the CIG online method for MECD estimation, to the problem of texture segmentation. Specifically, we consider mixtures of MGGD, as they have been successfully used to model wavelet statistics of texture images [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF][START_REF] Kwitt | Testing a multivariate model for wavelet coefficients[END_REF].

For this application, we randomly selected five pictures from the VisTex database [START_REF]Vision texture database[END_REF], fabric, white flowers, pink flowers, leafs and water. Each texture image was considered as an RGB 3-dimensional image, and modelled by a mixture of MGGD, whose parameters θ = (w, (µ k ) k , (Σ k ) k , (β k ) k ) k were estimated using five different methods, EMFP, SG, CIG-DS, CIG-OB and CIG-AS, already discussed in Subsection 7.2.2.

In order to evaluate the true performance of these estimation methods, we did not use special filters or feature extraction techniques. Rather, each given picture element was directly classified by its log-likelihood value. The sub-distribution (mixture component) with the largest log-likelihood value is the class to which the picture element is affected.

Each of the five methods was run for the same number of epochs (the number of epochs means the number of times that the complete dataset is traversed), and all were initialized using a k-means method.

For the first application, we considered a scenario with K = 2 textures. Figure 7.13 shows the visual segmentation results, and table 7.3.2 gives the accuracy of segmentation. In terms of accuracy, CIG-DS, CIG-AS and EMFP have the best performance (greater than 99%), followed by SG, while CIG-OB is the worst.

In terms of visual segmentation results, CIG-DS and AS are similar to the EMFPmethod. The two parts (fabric and flower) are completely and clearly segmented. However, in Figure 7.13(c), the results of SG show traces in the fabric were misidentified as flowers, and in Figure 7.13(e), the green leafs in the flowers were misclassified as fabric. For this scenario, it is clear that CIG-DS and CIG-AS are the most effective methods. white flowers are slightly misclassified. For SG method, in Figure 7.14(c), the fabric was completely misjudged as water, consequently, the accuracy was as low as 63.328% (see Table 7.3.2). For CIG methods, DS and AS also share the good performance of EMFP. In figures 7.14(d) and 7.14(f), fabric, leafs and pink flowers were clearly segmented. There are only some minor errors in water and white flowers. What is of note is that, with regard to the leafs texture, the performance of CIG-DS and CIG-AS is even better than EMFP.

Their accuracy is also quite good, with CIG-AS producing a similar accuracy to EMFP, i.e. 92%, and CIG-DS achieving an even better accuracy 96%. CIG-DS has almost achieved the same F1 score as EMFP with less computation time.

Conclusion

This final chapter compared the CIG methods, introduced in Chapters 5 and 6, to other, state-of-the-art methods, for ECD and MECD estimation, such as method of moments, fixed-point, expectation-maximization, and classical stochastic gradient [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF][START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF].

For ECD estimation, the CIG offline and CIG online methods were compared to the method of moments (MM) and to the fixed-point method (FP), as given in [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF] and [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF]. Figures 7.4 and 7.6 present the comparison, in terms of efficiency (roughly speaking, final estimation accuracy) and time consumption, respectively.

Figure 7.4 shows that the CIG offline method provides the best efficiency, being slightly better than the FP method. However, from Figure 7.6, it appears the CIG offline method is less attractive than the FP method, in terms of time consumption.

On the other hand, the same two figures (7.4 and 7.6) show the CIG online method is still comparable to the FP method, in terms of efficiency, but is highly attractive, in terms of time consumption, since it on par with the method of moments, and converges "instantly", even as the size of the dataset increases.

Overall, the CIG online method appears as the best choice for ECD estimation, among all methods under consideration.

For MECD estimation, CIG online methods were compared to two already existing methods, EMFP (Expectation-Maximisation with Fixed-Point [START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF][START_REF] Hasannasab | Alternatives to the EM algorithm for ML estimation of location, scatter matrix, and degree of freedom of the Student t distribution[END_REF]) and SG (Euclidean Stochastic Gradient [START_REF] Najar | Online recognition via a finite mixture of multivariate generalized Gaussian distributions[END_REF]).

First of all, Table 7.1 showed CIG online methods (CIG-DS and CIG-AS) had the highest rate of success among all methods considered. Under the same initialisation conditions, CIG online methods had the highest chance of converging to a local maximum of the log-likelihood function.

Figures 7.8 and 7.9 compared CIG online methods to EMFP and SG, in terms of their performance and time consumption. These figures showed the CIG-DS and CIG-AS methods were able to achieve the same performance as EMFP and SG, while being two to three times faster.

Overall, the CIG online methods also seemed the most attractive for MECD estimation, as they achieved comparable, or even better performance, to other methods, within a notably shorter time.

The applications to real data, discussed in Section 7.3, showed that CIG methods are able to tackle very large datasets (in the form of HD images, for example), and produce highly satisfactory results.

In conclusion, the numerical experiments and real-data applications in this chapter validate the methodology of CIG online methods, which may be regarded as very interesting new methods, to be considered by any users wishing to work with ECD or MECD.

Chapter 8 Conclusion and perspectives

In this thesis, an original geometric method was proposed, for online estimation of ECD and MECD models. This new method was named the CIG method (Component-wise Information Gradient method).

To begin, since the parameter spaces of ECD and MECD are Cartesian product spaces, a novel Riemannian metric was introduced on these spaces, and was named the CIM (Component-wise Information Metric). At each point of the parameter space, the CIM is a direct product of the Fisher Information Metrics, with respect to individual parameters (for example, location, scatter, and shape parameters).

This leads to the component-wise information gradient, which combines the information gradients with respect to these parameters. Based on the component-wise information gradient, the CIG method implements an alternating update scheme for the parameters, where each parameter is updated in its own turn.

The CIM is an easy-to-compute substitute, for the FIM (Fisher Information Metric). Similarly, the CIG method is a computationally-advantageous substitute for the Riemannian information gradient method, which would be very hard to implement, for complete ECD or MECD models.

For some particular cases of ECD models, the CIM and FIM coincide, so the CIG method shares many of the excellent properties of the Riemannian information gradient method. For instance, it achieves a mean-square convergence rate of order O(n -1 ), and it is asymptotically efficient, with an automatic choice of step-size.

For complete ECD or MECD models, the FIM is not tractable, and it is anyway different from the CIM. The CIG method shares some, but not all, of the properties of the Riemannian information gradient method. For example, it can still achieve a meansquare convergence rate of order O(n -1 ), but this requires a manual selection of step-size.

At the cost of some additional computations, this situation can be improved dramatically. The adaptive step-size version of the CIG method (proposed in Proposition 12), implements adaptive selection of step-size which provides a much faster convergence rate O(c -n ) (where c > 1).

Numerical simulations and real applications have proved that the CIG method has significant advantages over existing methods, especially in terms of time consumption. In future work, the problem of order selection will be included into our investigation of MECD estimation, using information criteria (AIC or BIC), or a non-parametric Bayesian approach. In addition, new classes of algorithms (such as stochastic EM and its many variants) will be studied and adapted to the estimation of ECD and MECD. Finally, the idea of CIG may be generalized to online estimation of any distributions of position-scale type, such as Riemannien Gaussian distribution and von Mise-Fisher distributions. Thus, the matrix H is equal to the identity matrix, and its smallest eigenvalue is λ = 1.

Proof of (i) : this follows directly from Propositions 2 and 3. Indeed, since λ = 1, the conditions of these propositions are verified, as soon as 2a > 1. Therefore, (4.9) and (4.10) hold true.

Proof of (ii) : this follows from Proposition 4. The conditions of this proposition are verified, as soon as 2a > 1. Therefore, the distribution of the re-scaled coordinates n 1/2 (θ (n) ) α converges to a centred d-variate normal distribution, with covariance matrix Σ given by Lyapunov's equation (4.11). If a = 1, then (A.19b) implies A α,β = -1 2 δ α,β , so that Lyapunov's equation (4.11) reads Σ = Σ * , as required.

For the following proof of (iii), the reader may wish to recall that summation convention is used throughout the present work. That is [START_REF] Petersen | Riemannian geometry[END_REF], summation is implicitly understood over any repeated subscript or superscript from the Greek alphabet, taking the values 1 , . . . , d .

Proof of (iii) : let (θ) = log L(θ) and assume U (θ, x) is given by (4.3). Then, by the definition of normal coordinates [START_REF] Petersen | Riemannian geometry[END_REF], the following expression holds

U α (θ * ) = ∂ ∂θ α θ * = 0 (A.20a)
Replacing this into (4.7) gives

Σ * α,β = E θ * ∂ ∂θ α ∂ ∂θ β θ * = 0 = -E θ * ∂ 2 ∂θ α ∂θ β θ α = 0 = ∂ 2 D ∂θ α ∂θ β θ α =0 (A.20b)
where the second equality is the so-called Fisher's identity (see [START_REF] Amari | Methods of Information Geometry[END_REF], Page 28), and the third equality follows from (4.2) by differentiating under the expectation. Now, by (4.8) and (A.19b), Σ * is the identity matrix.

To show that the online estimates θ (n) are asymptotically efficient, let (τ α ; α = 1, . . . , d) be any local coordinates with origin at θ * and let τ α n = τ α (θ (n) ) . From the second-order Taylor expansion of each coordinate function τ α , it is straightforward to show that

n 1/2 τ α n = ∂τ α ∂θ γ θ * n 1/2 (θ (n) ) γ + σ α (θ (n) ) n 1/2 d 2 (θ (n) , θ * ) (A.21a)
where the subscript θ * indicates the derivative is evaluated at θ * , and where σ α is a continuous function in the neighborhood of θ * . By (4.10), the second term in (A.21a) converges to zero almost surely. Therefore, the limit distribution of the re-scaled coordinates (n 1/2 τ α n ) is the same as that of the first term in (A.21a). By (ii), this is a centred d-variate normal distribution with covariance matrix Σ τ given by

Σ τ α,β = ∂τ α ∂θ γ θ * Σ * γ,κ ∂τ β ∂θ κ θ * = ∂τ α ∂θ γ θ * ∂τ β ∂θ γ θ * (A.21b)
where the second equality follows because Σ * γ,κ = δ γ,κ since Σ * is the identity matrix. It remains to show that Σ τ is the inverse of the Fisher information matrix I τ as in (A.34). According to (A.32), this is given by

I τ α,β = ∂ 2 D ∂τ α ∂τ β τ α =0 = -E θ * ∂ 2 ∂τ α ∂τ β τ α =0 = E θ * ∂ ∂τ α ∂ ∂τ β τ α =0 (A.21c)
where the second equality follows from (4.2), and the third equality from Fisher's identity (see [START_REF] Amari | Methods of Information Geometry[END_REF], Page 28). Now, a direct application of the chain rule yields the following

I τ α,β = E θ * ∂ ∂τ α ∂ ∂τ β τ α = 0 = ∂θ γ ∂τ α θ * E θ * ∂ ∂θ γ ∂ ∂θ κ θ γ = 0 ∂θ κ ∂τ β θ *
By the first equality in (A.20b), this is equal to

I τ α,β = ∂θ γ ∂τ α θ * Σ * γ,κ ∂θ κ ∂τ β θ * = ∂θ γ ∂τ α θ * ∂θ γ ∂τ β θ * (A.21d) because Σ * γ,κ = δ γ,
• c)(1) = (θ (n+1) ) α and (θ α • c)(0) = (θ (n) ) α . Moreover, dθ α (U n+1 ) = U α n+1 are the components of U n+1 . Replacing into (A.24a), this yields (θ (n+1) ) α = (θ (n) ) α + η (n+1) U α n+1 + [η (n+1) ] 2 π α n+1 (A.24c)
where π α n+1 is given by

π α n+1 = [η (n+1) ] 2 T α n+1 - 1 2 Γ α β,γ (θ (n) ) U β n+1 U γ n+1 (A.24d)
Comparing (A.24c) to (A.13a), it is clear the proof will be complete upon showing

E π α n+1 = O(n -1/2
E π α n+1 η (n+1) K 1 E U n+1 3 + d 2 × K 2 E d(θ (n) , θ * ) U n+1 2 (A.25a)
where d is the dimension of the parameter space Θ. However, using the fact that the x n are i.i.d. with distribution P θ * ,

E U n+1 3 X n = E θ * U (θ (n) , x) 3 R 3/4 (θ (n) ) (A.25b)
by (u2) and Jensen's inequality [START_REF] Shiryaev | Probability[END_REF]. On the other hand, by the Cauchy-Schwarz inequality,

E d(θ (n) , θ * ) U n+1 2 E d 2 (θ (n) , θ * ) 1/2 E U n+1 4 1/2 b n -1/2 E U n+1 4 1/2
for some b > 0 as follows from (4.9). Then, by the same reasoning that lead to (A.25b),

E d(θ (n) , θ * ) U n+1 2 ≤ b n -1/2 E R(θ (n) ) 1/2 (A.25c)
By (u2), there exists a uniform upper bound M on R(θ) for θ ∈ Θ * . Since θ (n) lies in Θ * for all n, it follows by replacing the inequalities (A.25b) and (A.25c) into (A.25a) that

E π α n+1 η (n+1) K 1 M 3/4 + d 2 × K 2 b n -1/2 M 1/2 (A.25d)
Finally, by recalling that η (n) = a n , it is clear that the right-hand side of (A.25d) is O(n -1/2 ), so the proof is complete. This follows immediately from the second-order Taylor expansion of D(θ), since θ * is a minimum of D(θ), by using (A.32). Formula (A.33) shows that the divergence D(θ) is equivalent to half the squared Riemannian distance d 2 (θ, θ * ), at θ = θ * .

The scalar products appearing in (A.32) form the components of the Fisher information matrix I τ of the coordinates τ α ,

I τ αβ = ∂ 2 D ∂τ α ∂τ β τ α =0
In any change of coordinates, these transform like the components of a (0, 2) covariant tensor [START_REF] Petersen | Riemannian geometry[END_REF]. That is, if (θ α ; α = 1, . . . , d ) are any local coordinates defined at θ * ,

I τ αβ = ∂θ γ ∂τ α θ * I θ γκ ∂θ κ ∂τ β θ *
where the subscript θ * indicates the derivative is evaluated at θ * , and where I θ γκ are the components of the Fisher information matrix I θ of the coordinates θ α .

The recursive estimates θ (n) are said to be asymptotically efficient, if they are asymptotically efficient in any local coordinates τ α , with origin at θ * . That is, according to the classical definition of asymptotic efficiency [START_REF] Ibragimov | Statistical estimation-asymptotic theory[END_REF][START_REF] Aad | Asymptotic statistics[END_REF], if the following weak limit of probability distributions is verified [START_REF] Shiryaev | Probability[END_REF],

L (n 1/2 τ α n ) =⇒ N d (0, Σ τ ) Σ τ = (I τ ) -1 (A.34)
where L{. . .} denotes the probability distribution of the quantity in braces, τ α n = τ α (θ (n) ) are the coordinates of the recursive estimates θ (n) , and N d (0, Σ τ ) denotes a centred dvariate normal distribution with covariance matrix Σ τ .

It is important to note that asymptotic efficiency of the recursive estimates θ (n) is an intrinsic geometric property, which does not depend on the particular choice of local coordinates τ α , with origin at θ * . This can be seen from the transformation rule of the components of the Fisher information matrix, described above. In fact, since these transform like the components of a (0, 2) covariant tensor, the components of Σ τ transform like those of a (2, 0) contravariant tensor, which is the correct transformation rule for the components of a covariance matrix. related, note that U (θ (0) ), ∇ θ D(θ (0) ) = -∇ * µ D(µ (0) , Σ (0) , β (0) ) 2 * U (θ (1) ), ∇ θ D(θ (1) ) = -∇ * Σ D(µ (1) , Σ (0) , β (0) )

2 * U (θ (2) ), ∇ θ D(θ (2) ) = -∇ * β D(µ (1) , Σ (1) , β (0) )

2 * and so on, for n 3. Therefore, U (θ (n) ), ∇ θ D(θ (n) ) < 0 whenever U (θ (n) ) = 0. This means that U (θ (i) ), ∇ θ D(θ (i) ) < 0 for the whole subsequence U (θ (i) ) i 0 , and this subsequence is indeed gradient related.

B.2 Proof of Proposition 7

The proof is a direct application of Remark 2, concerning Proposition 1, in [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF]. According to this remark, if U θ (n) , X (n+1) mb denotes the direction of descent, and if

E n U (θ (n) , X (n+1) mb
), ∇ θ D(θ (n) ) < 0, almost surely, for n > 0 (B.1)

where E i denotes conditional expectation with respect to (X

mb , X

mb , • • • , X

mb , • • • ), then lim θ (n) = θ * almost surely. Here (compare to the proof of Proposition 6), the direction of descent is given by U (θ (0) , X

mb ) = ∇ * µ p (µ (0) , Σ (0) , β (0) ; X

mb ) U (θ (1) , X

(2) mb ) = ∇ * Σ p (µ (1) , Σ (0) , β (0) ; X

mb ) U (θ (2) , X

(3) mb ) = ∇ * β p (µ (1) , Σ (1) , β (0) ; X

mb ) and so on. Therefore, the expectations in (B.1) can be found from

E 0 U (θ (0) , X (1)
mb ), ∇ θ D(θ (0) ) = -∇ µ D(θ (0) ) 2 * (B.2a) E 1 U (θ (1) , X A1. θ * is the unique stationary point of D(θ) in Θ * . Moreover, the second derivatives ∇ * 2 µ D(θ * ), ∇ * 2 Σ D(θ * ), ∇ * 2 β D(θ * ) are all positive-definite. We now show that Assumption A1 guarantees it is strictly negative, as required.

By assumption A1, and a direct application of the implicit function theorem (after taking Θ * sufficiently small) [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF], there exist three submanifolds of Θ, which pass through θ * ,

H µ = {θ ∈ Θ * : ∇ * µ D(θ) = 0} H Σ = {θ ∈ Θ * : ∇ * Σ D(θ) = 0} H β = {θ ∈ Θ * : ∇ * β D(θ) = 0}
Moreover, H µ is uniquely parameterised by (Σ, β), H Σ is uniquely parameterised by (µ, β) and H β is uniquely parameterised by (µ, Σ). In particular, this implies that H µ , H Σ and H β are submanifolds of lower dimension, and therefore have zero Riemannian volume.

The expectations in (B.2) are therefore almost surely all strictly negative. Indeed, each θ n with n > 0 has a probability density function with respect to Riemannian volume, and therefore almost surely does not belong to H µ , H Σ or H β .

B.3 Proof of Propositions 8 and 9

As for Proposition 7, this is an application of Remark 2 in [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF]. According to this remark, in order to obtain the mean-square rate and the asymptotic normality, it is enough to show the mean vector field X(θ) = E θ * [U (θ, x)] has an attractive stationary point at θ = θ * . Since U (θ, x) = ∇ θ (θ; x)

E θ * [U (θ, x)] =   -∇ * µ D(θ) -∇ * Σ D(θ) -∇ * β D(θ)   (B.3)
The covariant derivative of this vector field at the point θ = θ * is equal to the Hessian H(θ * ), which is positive-definite. Therefore, the results of Propositions 8 and 9 follow by Remark 2 in [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF].

B.4 Proof of Proposition 10

For the case of θ = (Σ), the geodesic convexity of the cost function D(θ) (or of D(θ)) follows by provingp (θ; x) is geodesically strictly convex in θ = (Σ) for any x.

Recall that, geodesic curves on P m are of the form [START_REF] Pennec | A Riemannian framework for tensor computing[END_REF] γ :R → P m t →A exp(tr)A † (B.4)

where exp denotes the matrix exponential map, A is an invertible matrix, and r is a diagonal matrix, both of same size as Σ. Then,p (θ) is geodesically convex if and only if the composition (p • γ)(t) is always a convex function with respect to t. Moreover, geodesic strict convexity is defined in exactly the same way. The composition (p • γ)(t) can be expressed where µ = m i=1 u 2 i δ r i , and δ r i is the Dirac measure concentrated at r i . Assume that the function f verifies Condition (5.9). Then, since ϕ is strictly logconvex, f • ϕ is strictly log-convex. Thus, the term log [(f • ϕ)(t)] of (B.5) is a strictly convex function of the real variable t. Since the term tr(r) t 2 of (B.5) amounts to an affine function of t, it is now clear that (p • γ)(t) is a strictly convex function of the real variable t, for any geodesic curve γ : R → P m . Finally, since x was chosen arbitrarily, p (θ; x) is geodesically strictly convex in θ = (Σ) for each x. Therefore, D(θ) and D(θ) are both geodesically strictly convex.

(-p • γ)(t) =

B.5 Proof of Corollary 2 and 3

For the case of θ = (Σ), note that ϕ : R → R + is strictly log-convex if and only if ϕ(t) = exp(ψ(t)) where ψ : R → R is strictly convex. 1) plugging equation (5.8a) into equation (5.9),

log(f • ϕ)(t) = 1 2 exp (β (ψ(t))) (B.8)
Therefore, condition (5.9) is verified since β > 0.

2) plugging equation (5.8b) into equation (5.9),

log (f • ϕ) (t) = β + m 2 log 1 + exp(ψ(t)) β (B.9)
Therefore, condition (5.9) is verified since β + m > 0.

For the case of θ = (µ, Σ), as mentioned above, the function f is reformulated. Then, the same strategy is applied for this reformulated f. 1) For MGGD, recall the geodesic curve for reformulated matrix S(t), S(t) = B exp(st)B † where exp denotes the matrix exponential map, B is an invertible matrix, and s is a diagonal matrix, both of same size as S. This function is proved to be log-convex in equation (B.6). Therefore, condition (5.9) is verified since β > 0 for MGGD model.

2) For Student T, plugging (5.10) into (5.9),

log f • ϕ (t) = β + m 2 1 - 1 β + 1 β exp(ψ(t)) (B.13)
Therefore, condition (5.9) is verified since β > 0.

where

θ (0) = r (0) , µ (0) k k , Σ (0) k k , β (0) k k
θ (1) = r (1) , µ

(0) k k , Σ (0) k k , β (0) k k
θ (2) = r (1) , µ

, Σ

, β (0) k k θ (3) = r (1) , µ

, Σ

, β

(0) k k (C.3)
Recall that D(θ) = -E [ f (θ; x)] and the form of ∇ θ D(θ) is defined as

∇ θ D(θ) =     ∇ * r D(θ) ∇ * µ k D(θ) k ∇ * Σ k D(θ) k ∇ * β k D(θ) k     (C.4)
Therefore, for step 0 to 3, the expectation in (C.1) are

E U θ (0) , x , ∇ θ D θ (0) θ = -∇ r D θ (0) * r 2 E U θ (1) , x , ∇ θ D θ (1) θ = - K k=1 ∇ µ k D θ (0) * µ k 2 E U θ (2) , x , ∇ θ D θ (2) θ = - K k=1 ∇ Σ k D θ (0) * Σ k 2 E U θ (3) , x , ∇ θ D θ (3) θ = - K k=1 ∇ r D θ (0) * β k 2 (C.5)
And so on, for t > 3, all these verify condition (C.1). According to remark 2 in [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF], in order to obtain the asymptotically linear convergence rate, it is enough to show the mean vector field X(θ) = E θ * [U (θ, x)] has an attractive stationary point at θ = θ * . Since U (θ, x) = -∇ θ f (θ; x)

E θ * [U (θ, x)] =     ∇ * r D(θ) ∇ * µ k D(θ) k ∇ * Σ k D(θ) k ∇ * β k D(θ) k     (C.6)
Under the same assumptions of B.2, the covariant derivative of this vector field at the point θ = θ * is equal to the Hessian H(θ * ), which is positive-definite. Therefore the linear convergence rate holds.

C.2 Proof of proposition 12

Note Θ * is the neighborhood of θ * which satisfies all assumptions in proposition 12. The objective function D(θ) is geodesically strongly convex and L-lipschitz smooth in Θ * . Recall the geodesically L-lipschitz

D(θ (n+1) ) D(θ (n) ) -η (n) ∇ θ D(θ (n) ), U (θ (n) , X (n) mb ) θ (n) + L 2 d 2 θ (n) , θ (n+1) (C.7)
Where •, • denotes the Riemannian metric in the general sense, therefore ∇ θ D(θ (n) ) denotes the general Riemannian gradient, and d(•, •) denotes the general Riemannian distance. The vector U (θ (n) , X

(n) mb ) denotes the direction of descent, and η (n) is the stepsize.

The set X

(n)

mb is supposed to be a random mini-batch that is uniformly selected from the complete dataset. Therefore

E X mb U (θ (n) , X (n) mb ) = U θ (n) , X (C.8)
here X denotes the complete dataset. Take expectation with respect to X mb

E X mb D(θ (n+1) ) D(θ (n) ) -η (n) ∇ θ (θ (n) ), U (θ (n) , X ) θ (n) + L 2 E X mb d 2 θ (n) , θ (n+1)
(C.9)

Expand the last term

d 2 (θ (n) , θ (n+1) ) =d 2 (r (n) , r (n+1) ) + K k=1 d 2 (µ (n) k , µ (n+1) k ) + d 2 (Σ (n) k , Σ (n+1) k ) +d 2 (β (n) k , β (n+1) k
d 2 (Σ (n) k , Σ (n+1) k ) = Exp -1 Σ (n) k (Σ (n+1) k ) 2 = [η (n) ] 2 U Σ k θ (n) , X (n) mb 2 (C.11c) d 2 (β (n) k , β (n+1) k ) = Exp -1 β (n) k β (n+1) k 2 = [η (n) ] 2 U β k (θ (n) , X (n) mb ) 2 (C.11d) therefore E X mb D θ (n+1) D(θ (n) ) -η (n) ∇ θ D(θ (n) , X ), U (θ (n) , X ) + L[η (n) ] 2 2 E X mb U (θ (n) , X (n) mb ) 2 (C.12)
For r ∈ S K-1 , the classic Euclidean gradient coincides with the Riemannian information gradient (3.29). k . Then for Σ k , its information gradient is

U r (θ (n) , X (n) mb ) 2 = ∇ * r D(θ (n) , X (n) mb ) 2 = ∇ r D(θ (n) , X (n) mb ) 2 (C.13) For µ k in R p , U µ k (θ (n) , X (n) mb ) 2 = I -1 µ k Σ (n) k ∇ µ k D(θ (n) , X (n) mb ) 2 [λ (n) max,k ] 2 I 2 µ k ∇ µ k D θ (n) , X (
U Σ k (θ (n) , X (n) mb ) =J Σ k ,1 ∇ ↑ Σ k D(θ (n) , X (n) mb ) ⊥ + J Σ k ,2 ∇ ↑ Σ k D(θ (n) , X (n) mb ) =J Σ k ,1 ∇ ↑ Σ k D(θ (n) , X (n) mb ) + J Σ k ,2 -J Σ k ,1 m tr [Σ (n) k ] -1 ∇ ↑ Σ k D(θ (n) , X (n) mb ) Σ (n) k (C.15)
Using the affine invariant metric, the square norm of

U Σ k (θ (n) , X (n) mb ) is U Σ k (θ (n) , X (n) mb ) 2 ↑ = tr (Σ (n) k ) -1 U Σ k (θ (n) , X (n) mb ) 2 (C.16)
After some necessary simplifications, the following relation ship could be summarized

U Σ k θ (n) , X (n) mb 2 ↑ J 2 Σ k ,1 ∇ ↑ Σ k D(θ (n) , X (n) mb ) 2 ↑
(C.17)

Finally, for β k , it is easy to obtain

U β k θ (n) , X (n) mb 2 = I -1 β k ∇ β k D(θ (n) , X (n) mb ) 2 = I -2 β k ∇ β k D(θ (n) , X (n) mb ) 2 (C.18)
We have then Then, using Lemma 1 of [START_REF] Yi Meng | Fast and furious convergence: Stochastic second order methods under interpolation[END_REF], we could obtain

E X mb D θ (n+1) D θ (n) -η (n) ∇ θ D(θ (n) , X ), U (θ (n) , X ) + L [η (n) ] 2 τ (n) max 2 E X mb ∇ θ D(θ (n) , X (n
E X mb ∇ θ D(θ (n) , X (n) mb ) 2 ρ ∇ θ D(θ (n) , X ) 2 (C.21)
where ρ = (n-b)(ρ-1)

(n-1)b + 1, and b is the size of mini-batch, n is the size of dataset. The expected decrease becomes

E X mb D(θ (n+1) ) D(θ (n) ) -η (n+1) ∇ θ D(θ (n) , X ), U (θ (n) , X ) + L [η (n+1) ] 2 τ (n) max ρ 2 ∇ θ D(θ (n) , X ) 2 (C.22)
Then, for item ∇ θ D(θ (n) , X ), U (θ (n) , X ) , starting with r ∇ r D(θ (n) , X ), U r (θ (n) , X ) = ∇ r D(θ (n) , X )

2 (C.23) For µ k ∇ µ k D(θ (n) , X ), U µ k (θ (n) , X ) λ min,k I µ k ∇ µ k D(θ (n) , X ) 2 (C.24)
where λ min,k is the smallest eigenvalue of Σ (n)

k . For scatter matrix Σ k ∇ Σ k D(θ (n) , X ), U Σ k (θ (n) , X ) J Σ k ,2 ∇ Σ k D(θ (n) , X ) 2 ↑ (C.25)
Finally, the shape parameter β k Take the negative on both sides of this inequality to get

∇ β k D(θ (n) , X ), U β k (θ (n) , X ) = I -1 β k ∇ β k D(θ (n) , X ) 2 (C.26) Therefore E X mb D(θ (n+1) ) D(θ (n) ) -η (n+1) τ (n) min ∇ θ D(θ (n) , X ) 2 + L [η (n+1) ] 2 τ (n) max ρ 2 ∇ θ D(θ (n) , X )
f (x) -f (x * ) -Exp -1 x (x * ) , ∇ x f (x) x - α 2 Exp -1 x (x * ) 2 x = - α 2 Exp -1 x (x * ) + 1 α ∇ x f (x) 2 x + 1 2α ∇ x f (x) 2 x 1 2α ∇ x f (x) 2 x (C.33)
where there is the result.

2. 1

 1 State of the art: existing methods for ECD estimation . . . . . . . . . . . . 2.2 State of the art: existing methods for MECD estimation . . . . . . . . . . 3.1 The CIM and related geometric objects for ECD and MECD . . . . . . . . 5.1 Global convergence analysis: MGGD . . . . . . . . . . . . . . . . . . . . . 5.2 Global convergence analysis: Student T . . . . . . . . . . . . . . . . . . . . 7.1 MECD: percentage of successful runs . . . . . . . . . . . . . . . . . . . . . 7.2 MECD: Accuracy of segmentation (Texture 1) . . . . . . . . . . . . . . . . 7.3 MECD: Accuracy of segmentation (Texture 2) . . . . . . . . . . . . . . . .

  mb mini-batch p(x|θ) or p θ density of an ECD with parameter θ p log-likelihood function of an ECD θ parameter of an ECD Θ parameter space of an ECD µ location parameter Σ scatter matrix β shape parameter δ (x -µ) † Σ -1 (x -µ) f (x|θ) or f θ density of a MECD with parameter θ f log-likelihood function of an MECD θ parameter of a MECD Θ parameter space of a MECD w the mixture weights θ k the parameter of the k-th component v, U , or U (θ, x) tangent vector, or descent direction P = (P, Θ, X) statistical (parametric) model with parameter θ D(θ) KL divergence D(θ||θ * ) θ *the true parameter•, • or •, • θ Riemannian metric in general •, • ↑ or •, • ↑ θ affine-invariant metric (at point θ) •, • * or •, • * θ Fisher information metric (at point θ) •, • or •, • θ component-wise information metric (at point θ) d(•, •) Riemannian distance in general d ↑ (•, •) affine-invariant distance d * (•, •) information distance d (•, •) component-wise Information distance grad θ classic Euclidean gradient of function with respect to θ ∇ θ Riemannian gradient of with respect to θ in general ∇ ↑ θ affine-invariant gradient of function with respect to θ ∇ * θ information gradient of function with respect to θ ∇ θ

4. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.3 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.4 Application to scatter-matrix estimation . . . . . . . . . . . . . . . . . . . 53 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

Figures 4 . 2 m 2 +

 422 1 and 4.2 below display numerical results from an application to MGGD, which correspond to h(δ) = -δ β 2 in (4.14) and A = m β in (4.16)[START_REF] Berkane | Geodesic estimation in elliptical distributions[END_REF][START_REF] Wang | Symmetric multivariate and related distributions[END_REF]. These figures were generated from 10 3 Monte Carlo runs of the algorithm defined by (4.1a) and (4.17), with random initialisations, for the specific values β = 4 and m = 7. Essentially the same numerical results could be observed for any β 9 and m 50.

Figures 4 .

 4 1 and 4.2 are concerned with the Riemannian information distance d * (θ (n) , θ * ) induced by the Fisher information metric

Figure 4 .

 4 Figure 4.1 confirms the fast rate of convergence (4.9), stated in (i) of Proposition 5. On a log-log scale, it shows the empirical mean E M C d 2 * (θ (n) , θ * ) over Monte Carlo runs, as a function of n. This decreases with a constant negative slope equal to -1, starting roughly at log n = 4.Figure 4.2 confirms the asymptotic efficiency of the online estimates θ (n) , stated in (iii) of Proposition 5, using Remark 4.Figure 4.2 shows a kernel density estimate of n d 2 * (θ (n) , θ * ) where n = 10 5 (solid blue curve). This agrees with a X 2 -distribution with 28 degrees of freedom (dotted red curve), where d = 28 is indeed the dimension of the parameter space Θ = P m for m = 7.

Figure 4 .

 4 Figure 4.1 confirms the fast rate of convergence (4.9), stated in (i) of Proposition 5. On a log-log scale, it shows the empirical mean E M C d 2 * (θ (n) , θ * ) over Monte Carlo runs, as a function of n. This decreases with a constant negative slope equal to -1, starting roughly at log n = 4.Figure 4.2 confirms the asymptotic efficiency of the online estimates θ (n) , stated in (iii) of Proposition 5, using Remark 4.Figure 4.2 shows a kernel density estimate of n d 2 * (θ (n) , θ * ) where n = 10 5 (solid blue curve). This agrees with a X 2 -distribution with 28 degrees of freedom (dotted red curve), where d = 28 is indeed the dimension of the parameter space Θ = P m for m = 7.

Figure 4 .

 4 Figure 4.1 confirms the fast rate of convergence (4.9), stated in (i) of Proposition 5. On a log-log scale, it shows the empirical mean E M C d 2 * (θ (n) , θ * ) over Monte Carlo runs, as a function of n. This decreases with a constant negative slope equal to -1, starting roughly at log n = 4.Figure 4.2 confirms the asymptotic efficiency of the online estimates θ (n) , stated in (iii) of Proposition 5, using Remark 4.Figure 4.2 shows a kernel density estimate of n d 2 * (θ (n) , θ * ) where n = 10 5 (solid blue curve). This agrees with a X 2 -distribution with 28 degrees of freedom (dotted red curve), where d = 28 is indeed the dimension of the parameter space Θ = P m for m = 7.

Figure 4 . 1 :

 41 Figure 4.1: Riemannian information gradient: Fast rate of convergence

Figure 4 . 2 :

 42 Figure 4.2: Riemannian information gradient: Asymptotic efficiency

Table 5 . 1 :Table 5 . 2 :

 5152 Global convergence analysis: MGGD Student θ = (Σ) global convergence holds for β > 0 θ = (µ, Σ) global convergence holds β > 0 Global convergence analysis: Student T For the cases indicated in Tables 5.1 and 5.2, the cost function ( D(θ) for the offline version, and D(θ)

7. 2 . 1 ECD

 21 This subsection presents a set of computer experiments, which confirm the theoretical results of Sections 5.2 and 5.3, and provide a detailed comparison of the CIG estimation methods, with the already existing MM and FP[START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF]. For every experiment, 1000 Monte Carlo trials were carried out. For each trial, the dataset X = {x 1 , • • • , x T } is independent and identically distributed, according to true parameters (µ * , Σ * , β * ). The dimension m of x t is taken equal to 10. The true µ * is randomly chosen from a multivariate normal distribution. The scatter Σ * is defined as Σ(i, j) = ρ |i-j| for i, j ∈ {1, m}, and ρ uniformly distributed in [0.2, 0.8]. The shape parameter β * is uniformly selected from the interval [0.2, 5], for MGGD and for Student t.

Figure 7 . 1 :

 71 Figure 7.1: ECD: The superlinear convergence rate for CIG offline

  Figures 7.2(a), 7.2(b), and 7.2(c) display the mean-square rate of convergence, stated theoretically in Proposition 8. In these log-log plots, the x-axis

Figure 7 . 2 :

 72 Figure 7.2: ECD: O(n -1 ) convergence rate for CIG Online method

( a )Figure 7 . 3 :

 a73 Figure 7.3: ECD: Asymptotic normality of CIG online method For the case of θ = (Σ) and θ = (µ, Σ), Figures 7.3(a) and 7.3(b) demonstrate the asymptotic normality of CIG online method, obtained theoretically in Corollary 1. The samples being vectors of dimension m = 10, the dashed blue curve is the probability density of a chi-squared distribution with 55 and 65 degrees of freedom, for Figures 7.3(a) and 7.3(b), respectively. The solid lines are the smoothed histograms of N d 2 * (θ * , θ (N )) where N = 10 5 . These "estimated p.d.f." coincide very closely with the theoretical chi-squared probability density, confirming the fact that θ (N ) is asymptotically normally distributed about θ * .In the third experiment, we compare the efficiency of the CIG offline and online methods with other common estimation methods, MM (method of moments) and FP (fixedpoint method). In each trial, the dataset is generated from an MGGD model. For MGGD, the MM was given in[START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF], and the FP method in[START_REF] Pascal | Parameter estimation for multivariate generalized Gaussian distributions[END_REF]. In Figures 7.4(a), 7.4(b) and 7.4(c),

3

 3

  (a) and 7.3(b), respectively. The solid lines are the smoothed histograms of N d 2 * (θ * , θ (N )

Figure 7 . 4 :

 74 Figure 7.4: ECD: Efficiency comparison between MM, FP and CIG methods

  Figure 7.5(a) shows that the norm of the component-wise information gradient of the Kullback-Leibler divergence (3.1) always converge to 0, independently of the initial value θ (0) . Here, this component-wise information gradient is estimated empirically based on the Equation (5.2).

Figure 7 . 5 :

 75 Figure 7.5: Complete ECD: existence of local minima

Figure 7 .

 7 6 refers to the same experiment as in Figure 7.4. In Figures 7.6(a), 7.6(b) and 7.6(c), the x-axis denotes the size of the dataset, and the y-axis denotes the computation time necessary to achieve convergence.

Figure 7 . 6 :

 76 Figure 7.6: ECD: Time consumption until convergence

Figure 7 . 7 :

 77 Figure 7.7: MECD: The convergence rate of CIG-DS and CIG-AS

Figure 7 . 8 :

 78 Figure 7.8: MECD: Number of epochs versus log-likelihood 1

Figure 7 . 9 :

 79 Figure 7.9: MECD: Time-consumption versus log-likelihood

7. 3

 3 Applications to real data 7.3.1 Colour transformation with MGGD

Figure 7 . 12 :

 712 Figure 7.10: 3D colour transformation

Figure 7 . 13 :

 713 Figure 7.13: MECD: Segmentation with K = 2 textures For the second application, the number of components is K = 5. The textures of leafs and water were added. The visual results are presented in Figure 7.14(b). For this more complicated cas, the five segmentations are evaluated by the confusion matrices, and the confusion matrices are summarized in Table 7.3.2 with accuracy and F1 score (weighted average).According to Table7.3.2 and Figure7.14(b), EMFP is always rather satisfactory. The fabric and the pink flowers are totally identified, while some traces in leafs, water and

Figure 7 . 14 :

 714 Figure 7.14: MECD: Segmentation with K = 5 textures

( 2 )

 2 mb ), ∇ θ D(θ(1) ) = -∇ Σ D(θ (1) ) 2 * (B.2b) E 2 U (θ (2) , X (3) mb ), ∇ θ D(θ (2) ) = -∇ β D(θ (2) ) 2 * (B.2c)and so on, for n 3. Thus, the expectation in (B.1) is always negative. Recall the assumption A1 in Proposition 7.

  log det(A) + 1 2 tr(r)t + log [(f • ϕ)(t)] (B.5)recall that, the funtion f is defined in equation (5.7) and (5.8), and hereϕ(t) = m i=1 u 2 i exp(-r i t) (B.6)with u = A -1 x, its components u i , and r i are the diagonal elements of r. The function ϕ : R → R + is strictly log-convex, because it is the Laplace transform of a positive measure[START_REF] Shiryaev | Probability[END_REF] ϕ(t) = ∞ 0 exp(-tx)µ(dx) (B.7)

1 ) ) such that m+1 i=1 v 2 i e -s i t = m+1 i=1 w 2 i e -q i t + 1 (

 11 δ y (t) = y † S -1 y = m+1 i=1 v 2 i e -s i t with v = B -1 y (B.10)According to equation (B.10), we have δ y > 1. Therefore, ∃w ∈ R m+1 and ∃q ∈ (0, +∞) m+1 (e.g. w = (u, 0) and q = (r,

)(C. 10 ) 2 =

 102 Recall the square distance in each sub-spaced 2 (r (n) , r (n+1) ) = Exp -1 r (n) r (n+1) 2 = [η (n) ] 2 U r (θ (n) , X [η (n) ] 2 U µ k (θ (n) , X

n) mb 2 (

 2

) mb ) 2 ( 2 Σ

 22 k ,1 k , I -2 β k k . Recall the definition of ρ-strong growth condition E X mb ∇ θ D(θ, X mb ) 2 ρ ∇ θ D(θ, X ) 2 (C.20)

2 (,E 2 (C. 28 ) 2 d 2

 222822 (J Σ k ,2 ) k , I -1 β k k . Since η (n+1) = τ X mb D(θ (n+1) ) D(θ (n) )max ∇ θ D(θ (n) , X ) Using the Polyak-Lojasiewicz inequality in Riemannian context (which is proved in section C.3) ∇ θ D(θ (n) , X ) 2α D(θ (n) ) -D(θ * ) , we have E X mb D(θ (n+1) ) D(θ (n) ) -(n) ) -D(θ * ) (C.29)substracting D(θ * ) from both sidesE X mb D(θ (n+1) ) -D(θ * ) (n) ) -D(θ * ) (C.30)Range these N times of iterationE X mb D(θ (N ) ) -D(θ * ) (0) ) -D(θ * ) (C.31)C.3 Proof of Polyak-Lojasiewicz inequality in Riemannian contextConsider a continuous and differentiable function defined on a Riemannian manifold f : M → R. Note x * = arg inf x∈M f (x) its minimum, there exists a neighborhood M * ⊂ M of x * such that the function f is strongly geodesically convex in M * . For any x ∈ M * , ∃α > 0 such thatf (x * ) -f (x) Exp -1 x (x * ) , ∇ x f (x) x + α (x, x * ) (C.32)

  Proof of Proposition 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Proof of Proposition 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Proof of Propositions 8 and 9 . . . . . . . . . . . . . . . . . . . . . . . . . B.4 Proof of Proposition 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.5 Proof of Corollary 2 and 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . Proof of proposition 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2 Proof of proposition 12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.3 Proof of Polyak-Lojasiewicz inequality in Riemannian context . . . . . . .

	B Proof for chapter 5
	Table of acronyms B.1 C Proofs for Chapter 6
	C.1
	Acronym Meaning
	ECD	Elliptically Contoured Distribution
	MGGD	Multivariate Generalized Gaussian Distribution
	MECD	Mixture of Elliptically Contoured Distributions
	FP	Fixed Point method
	EMFP	Expectation-Maximisation with Fixed Point method
	SG	Stochastic Gradient methode
	FIM	Fisher Information Metric
	CIM	Component-wise Information Metric
	CIG	Component-wise Information Gradient
	CIG-DS	Component-wise Information Gradient with Decreasing Step-size
	CIG-OB	Component-wise Information Gradient with Online Backtracking
	CIG-AS	Component-wise Information Gradient with Adaptive Step-size

  . Les sections 2.3 et 2.4 présentent les méthodes existantes pour l'estimation des lois à contours elliptiques et de leurs mélanges, respectivement. A la fin du chapitre, ces méthodes sont regroupées dans les tableaux 2.1 et 2.2, sous une forme synthétique.

• le chapitre 3 contient la contribution de la thèse à la géométrie des lois à contours elliptiques et de leurs mélanges : l'introduction d'une nouvelle métrique Riemannienne, appelée CIM (component-wise information metric), dans les sections 3.3 et 3.4. Ces deux sections donnent l'expression de la CIM, du gradient riemannien qui lui est associé (le CIG), et de la rétraction à utiliser avec ce gradient, pour l'estimation des lois à contours elliptiques et de leurs mélanges, respectivement. A la fin du chapitre, ces expressions sont réunies dans le tableau 3.5. Elles seront fondamentales pour les nouvelles méthodes d'estimation, introduites dans les chapitres 5 et 6.

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2 Geometry of ECD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.1 Affine-invariant metric . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.2 Fisher information metric . . . . . . . . . . . . . . . . . . . . . . . 28 2.3 Estimation of ECD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.1 Euclidean methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.2 Riemannian methods . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4 Estimation of MECD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.4.1 Euclidean methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.4.2 Riemannian methods . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

	Chapter 2
	• Chapter 7 illustrates, through applications to simulated and real data, the various properties of CIG methods, introduced in Chapters 5 and 6. It also compares State of the art
	CIG methods to existing state-of-the-art methods, in terms of performance and
	computation time. In this chapter, Section 7.2 presents computer experiments with
	simulated data, while Section 7.3 deals with applications to real data, such as image
	editing and texture segmentation.
	2.1
	Appendices A to C are devoted to the proofs of various propositions and corollaries, stated
	in Chapters 4 to 6.

Table 2 . 1 :

 21 and 2.4, in the following two tables. State of the art: existing methods for ECD estimation

	ECD

Table 2 . 2 :

 22 State of the art: existing methods for MECD estimation Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2 Geometry of (µ,Σ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3 Geometry of (µ,Σ,β) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.4 Mixtures of ECD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

	Hertrich, Johannes, and Gabriele Steidl. [38]	Li, Shengxi, Zeyang Yu, and Danilo
	"Inertial Stochastic PALM and its	Mandic. [51]
	Application for Learning Student T	"A universal framework for learning the
	Mixture Models." 2020	elliptical mixture model." 2020
	Method: Inertial PALM with SGD	Method: conjugate gradient
	Unknown parameters: (w, µ, Σ, β)	Unknown parameters: (w, µ, Σ)
	Roizman Violeta, Matthieu Jonckheere,
	and Frédéric Pascal [67]. "A flexible EM
	-like clustering algorithm for noisy data."
	2019

.2. MECD Euclidean Riemannian Najar, Fatma, et al. [57] "Online recognition via a finite mixture of MGGD." 2020 Method: EM-FP, SGD Unknown parameters: (w, µ, Σ, β) Hosseini, Reshad, and Suvrit Sra. [42] "An alternative to EM for Gaussian mixture models: batch and stochastic Riemannian optimization." 2020 Method:EM-FP,steepest gradient, BFGS,conjugate gradient Unknown Parameters: (w, µ, Σ)

  6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.2 CIG-DS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 6.3 CIG-AS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Table 7 .

 7

	1: MECD: percentage of successful runs
	method success rate
	EM-FP	83%
	SG	77%
	CIG-DS	90%
	CIG-OB	89%
	CIG-AS	90%

Table 7 . 3

 73

	method correct estimates
	EM-FP	99%
	SG	97%
	CIG-DS	99%
	CIG-OB	73%
	CIG-AS	99%
	Table 7.2: MECD: Accuracy of segmentation (Texture 1)

.2 with accuracy and F1 score (weighted average).

According to Table

7

.3.2 and Figure

7

.14(b), EMFP is always rather satisfactory. The fabric and the pink flowers are totally identified, while some traces in leafs, water and

Table 7 .

 7

			score (weighted average)
	EMFP	93%	0.93
	SG	63%	0.58
	CIG-DS	96%	0.92
	CIG-OB	85%	0.78
	CIG-AS	92%	0.88

3: MECD: Accuracy of segmentation (Texture 2)

  κ is the identity matrix. Comparing (A.21b) to (A.21d), it is clear that Σ τ is the inverse of the Fisher information matrix I τ as in (A.34). K 1 which does not depend on n, as can be shown by direct calculation. Of course, (θ α

	where T α n+1 satisfies	T α n+1	K 1 U n+1	3	(A.24b)
	for a constant
					Proof of (iv) : (4.12a)
	and (4.12b) follow from (4.9) and (4.10), respectively, by using (A.33). Precisely, it is
	possible to write (A.33) in the form

  ). To do so, note that each Christoffel symbol Γ α β,γ is a C 1 function on the compact set Θ * , with Γ α β,γ (θ * ) = 0 by the definition of normal coordinates[START_REF] Petersen | Riemannian geometry[END_REF].

	Therefore,
	Γ α β,γ (θ)	K 2 d(θ, θ * )	(A.24e)
	for a constant K

2 which does not depend on n. Replacing the inequalities (A.24b) and (A.24e) into (A.24d), and taking expectations, it follows that

• le chapitre 7 met en évidence les propriétés des méthodes CIG, à travers des applications à des données simulées et réelles. Il compare aussi les méthodes CIG aux méthodes existantes dans la littérature, en termes de leur performance et de leur temps de calcul. Dans ce chapitre, la section 7.2 présente les applications aux données simulées, alors que la section 7.3 s'adresse aux données réelles.

In later chapters of this thesis, we will consider an alternative formulation, which focuses on minimising the Kullback-Leibler divergence, rather than maximising the log-likelihood.

The number of epochs means the number of times that the complete dataset is traversed. For the EMFP method, which is a batch (offline) method, the number of epochs is equal to the number of complete parameter updates.

Remerciements

Chapter 7

Appendix A

Proofs of chapter 4

A.1 Proofs of the main results

A.1.1 Proof of Proposition 1

The proof is a generalisation of the original proof in [START_REF] Bonnabel | Stochastic gradient descent on Riemannian manifolds[END_REF], itself modeled on the proof for the Euclidean case in [START_REF] Saad | On-line learning in neural networks[END_REF]. Throughout the following, let X n be the σ-field generated by x 1 , • • • , x n [START_REF] Shiryaev | Probability[END_REF]. Recall that (x n ; n = 1, 2, • • • ) are i.i.d. with distribution P θ * . Therefore, by (4.1a), θ (n) is X n -measurable and x n+1 is independent from X n . Thus, using elementary properties of conditional expectation [START_REF] Shiryaev | Probability[END_REF],

where (A.1a) follows from (4.1c), and (A.1b) from (u1). Let L be a Lipschitz constant for ∇D(θ), and C be an upper bound on V (θ), for θ ∈ Θ * . The following inequality is now proved, for any positive integer n,

once this is done, Proposition 1 is obtained by applying the Robbins-Siegmund theorem [START_REF] Duflo | Random iterative models[END_REF].

Proof of (A.2) : let c(t) be the geodesic connecting θ (n) to θ (n+1) with equation

From the fundamental theorem of calculus,

Since the online estimates θ (n) are stable, θ (n) and θ (n+1) both lie in Θ * . Since θ * is convex, the whole geodesic c(t) lies in Θ * . Then, since ∇D(θ) is Lipschitz on Θ * , it follows from (A.4), D(θ (n+1) ) -D(θ (n) ) η (n+1) U (θ (n) , x n+1 ), ∇D(θ (n) )

Taking conditional expectations in this inequality, and using (A.1a) and (A.1b),

so (A.2) follows since (u1) guarantees V (θ (n) ) C. Conclusion: by the Robbins-Siegmund theorem, inequality (A.2) implies that, almost surely,

In particular, from the first condition in (4.1b), convergence of the sum in (A.7) implies

Now, since the sequence of online estimates θ (n) lies in the compact set Θ * , it has at least one point of accumulation in this set, say

where the third equality follows from (A.8). This means that θ * is a stationary point of D(θ) in Θ * . Thus, (d1) implies θ * = θ * is the unique point of accumulation of θ (n) . In other words, lim θ (n) = θ * almost surely.

A.1.2 Proof of Proposition 2

The proof is modeled on the proofs for the Euclidean case, given in [7,[START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF]. It relies on the following geometric Lemmas 1 and 2. Lemma 1 will be proved in Appendix A.2. On the other hand, Lemma 2 is the same as the trigonometric distance bound of [START_REF] Zhang | First-order methods for geodesically convex optimization[END_REF]. For Lemma 1, recall that λ > 0 denotes the smallest eigenvalue of the matrix H defined in (4.8).

Lemma 1 for any µ < λ, there exists a neighborhood Θ * of θ * , contained in Θ * , with

Lemma 2 let -κ 2 be a lower bound on the sectional curvature of Θ in Θ * , and C κ = Rκ coth(Rκ) where R is the diameter of Θ * . For τ, θ ∈ Θ * , where τ = Exp θ (u),

Proof of equation (4.9) : let η (n) = a n with 2λa > 2µa > 1 for some µ < λ, and let Θ * be the neighborhood corresponding to µ in Lemma 1. By Proposition 1, the θ (n) converge to θ * almost surely. Without loss of generality, it can be assumed that all the θ (n) lie in Θ * , almost surely. Then, (4.1a) and Lemma 2 imply, for any positive integer n,

Indeed, this follows by replacing τ = θ (n+1) and θ = θ (n) in (A.9b). Taking conditional expectations in (A.10a), and using (A.1a) and (A.1b),

Then, by (u1) and (A.9a) of Lemma 1,

where C is an upper bound on V (θ), for θ ∈ Θ * . By further taking expectations

Using (A.10c), the proof reduces to an elementary reasoning by recurrence. Indeed, replacing η (n) = a n into (A.10c), it follows that

Let b be sufficiently large, so (A.11b) is verified and

Then, by recurrence, using (A.11a) and (A.11b), one also has that E d 2 (θ (n) , θ * ) b(n) for all n n o . In other words, (4.9) holds true.

A.1.3 Proof of Proposition 3

The proof is modeled on the proof for the Euclidean case in [START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF]. To begin, let W n be the stochastic process given by

where q ∈ (0, 1 -p) (A.12a)

The idea is to show that this process is a positive supermartingale, for sufficiently large n. By the supermartingale convergence theorem [START_REF] Shiryaev | Probability[END_REF], it then follows that W n converges to a finite limit, almost surely. In particular, this implies

Then, L p must be equal to zero, since p is arbitrary in the interval (0, 1). Precisely, for any ε ∈ (0, 1 -p),

It remains to show that W n is a supermartingale, for sufficiently large n. To do so, note that by (A.10b) from the proof of Proposition 2,

Here, the first term on the right-hand side is negative, since 2µa > 1 > p. Moreover, the third term dominates the second one for sufficiently large n, since q < 1 -p. Thus, for sufficiently large n, the right-hand side is negative, and W n is a supermartingale.

A.1.4 Proof of Proposition 4

the proof relies on the following geometric Lemmas 3 and 4, which are used to linearise Algorithm (4.1a), in terms of the normal coordinates (θ (n) ) α . This idea of linearisation in terms of local coordinates also plays a central role in [START_REF] Tripuraneni | Averaging stochastic gradient descent on Riemannian manifolds[END_REF].

Lemma 3 let θ (n) , θ (n+1) be given by (4.1a) with η (n) = a n . Then, in a system of normal coordinates with origin at θ * ,

where U α n+1 are the components of U (θ (n) , x n+1 ).

Then, in a system of normal coordinates with origin at

where v α n are the components of v n and the H α,β were defined in (4.8).

Linearisation of (4.1a)

Denote the re-scaled coordinates n 1/2 (θ (n) ) α by g α n , and recall η (n) = a n . Then, using the estimate (n + 1)

) is a first-order, inhomogeneous, linear difference equation, for the "vector" g n of components g α n . Study of equation (A.14b) : switching to vector-matrix notation, equation (A.14b) is of the general form

where I denotes the identity matrix, A has matrix elements A α,β , and (ξ n ) is a sequence of inputs. The general solution of this equation is [START_REF] Kailath | Linear systems[END_REF][START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF]

where the transition matrix A n,k is given by

Since 2λa > 1, the matrix A is stable. This can be used to show that [START_REF] Kailath | Linear systems[END_REF][START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF] q > 1 2 and E ξ n = O(n -q ) =⇒ lim g n = 0 in probability (A.15d)

where |ξ n | denotes the Euclidean vector norm. Then, it follows from (A.15d) that g n converges to zero in probability, in each one of the three cases

Indeed, in the first two cases, the condition required in (A.15d) can be verified using (4.9), whereas in the third case, it follows immediately from the estimate of E|π α n+1 | in (A.13a).

Conclusion : by linearity of (A.14b), it is enough to consider the case ξ α n+1 = w α n+1 in (A.15a). Then, according to (A.15b), g n has the same limit distribution as the sums

By (A.1), (w k ) is a sequence of square-integrable martingale differences. Therefore, to conclude that the limit distribution of gn is a centred d-variate normal distribution, with covariance matrix Σ given by (4.11), it is enough to verify the conditions of the martingale central limit theorem [START_REF] Hall | Martingale limit theory and its application[END_REF],

where Σ k is the conditional covariance matrix

Conditions (A.17) are verified in Appendix A.3, which completes the proof.

A.1.5 Proof of Proposition 5

Denote ∂ α = ∂ ∂θ α the coordinate vector fields of the normal coordinates (θ) α . Since •, • coincides with the Fisher information metric of the model P , it follows from (4.8) and (A. [START_REF] Frontera-Pons | A class of robust estimates for detection in hyperspectral images using elliptical distributions background[END_REF],

However, by the definition of normal coordinates [START_REF] Petersen | Riemannian geometry[END_REF], the ∂ α are orthonormal at θ * . Therefore,

where ω is a continuous function in the neighborhood of θ * , equal to zero at θ = θ * . To obtain (4.12a), it is enough to take expectations in (A.22a) and note that ω is bounded above in the neighborhood of θ * . Then, (4.12a) follows directly from (4.9).

To obtain (4.12b), it is enough to multiply (A.22a) by n p where p ∈ (0, 1). This gives the following expression

From (4.10), n p d 2 (θ (n) , θ * ) converges to zero almost surely. Moreover, by continuity of ω, it follows that ω θ (n) converges to ω(θ * ) = 0 almost surely. Therefore, by taking limits in (A.22b), it is readily seen that

almost surely. However, this is equivalent to the statement that D(θ (n) ) = o(n -p ) for p ∈ (0, 1), almost surely. Thus, (4.12b) is proved.

A.2 Proofs of geometric lemmas

A.2.1 Lemma 1

Let c(t) be the geodesic connecting θ * to some θ ∈ Θ * , parameterised by arc length. In other words, c(0) = θ * and c(t θ ) = θ where t θ = d(θ, θ * ). Denote Π t the parallel transport along c(t), from T c(0) Θ to T c(t) Θ. Since the velocity ċ(t) is self-parallel [START_REF] Petersen | Riemannian geometry[END_REF],

Multiplying this identity by -t θ , it follows that

Moreover, recall the first-order Taylor expansion of the gradient ∇D(θ) [START_REF] Chavel | Riemannian geometry: a modern introduction[END_REF][START_REF] Petersen | Riemannian geometry[END_REF]

where φ(θ) is continuous and equal to zero at θ = θ * . Here, ∇ 2 D(θ * ) is the Hessian of D(θ) at θ = θ * , considered as a linear mapping of T θ * Θ [19, 62]

where ∇ w denotes the covariant derivative in the direction of w. By (d1), the first term on the right-hand side of (A.23b) is equal to zero, so that

Taking the scalar product of (A.23a) and (A.23c),

since parallel transport preserves scalar products. In terms of the normal coordinates θ α , this reads [62]

where H = (H α,β ) was defined in (4.8), θα denotes the quotient θ α /t θ , and the φ α denote the components of φ(θ). Note that t 2 θ = d 2 (θ, θ * ) = θ α θ α , so (A.23e) can be written

where ψ(θ) is continuous and equal to zero at θ = θ * . To conclude, let µ = λ -ε for some ε > 0, and Θ * a neighborhood of θ * , contained in Θ * , such that ψ(θ) ε for θ ∈ Θ * . Then, since λ is the smallest eigenvalue of H = (H α,β ),

for θ ∈ Θ * . This is exactly (A.9a), so the lemma is proved.

A.2.2 Lemma 3

To simplify notation, let U n+1 = U (θ (n) , x n+1 ). Then, the geodesic c(t), connecting θ (n) to θ (n+1) , has equation

Each one of the normal coordinates θ α is a C 3 function θ α : Θ * → R, with differential dθ α and Hessian [62] ∇ 2 θ α = -Γ α β,γ (θ) dθ β ⊗ dθ γ where Γ α β,γ are the Christoffel symbols of the coordinates θ α , and ⊗ denotes the tensor product. Then, the second-order Taylor expansion of the functions θ α • c reads

A.2.3 Lemma 4

The lemma is an instance of the general statement : let θ ∈ Θ * and v = ∇D(θ). Then, in a system of normal coordinates with origin at θ * ,

where v α are the components of v. Indeed, (A.13b) follows from (A.26a) after replacing θ = θ (n) , so that v = v n , and setting

To prove (A.26a), recall (A.23c) from the proof of Lemma 1, which can be written

Denote ∂ α = ∂ ∂θ α the coordinate vector fields of the normal coordinates θ α . Note that [START_REF] Chavel | Riemannian geometry: a modern introduction[END_REF][START_REF] Petersen | Riemannian geometry[END_REF] Exp

From the first-order Taylor expansion of the vector fields ∂ α [START_REF] Chavel | Riemannian geometry: a modern introduction[END_REF][START_REF] Petersen | Riemannian geometry[END_REF]

where χ α (θ) is continuous and equal to zero at θ = θ * . However, by the definition of normal coordinates [START_REF] Petersen | Riemannian geometry[END_REF], each covariant derivative ∇∂ α (θ * ) is zero. In other words,

Now, to obtain (A.26a), it is enough to note the decomposition v = v α ∂ α (θ) is unique, and φ(θ) -H αβ θ β χ α (θ) converges to zero as θ converges to θ * .

A.3 Conditions of the martingale CLT

For the verification of Conditions (A.17), the following inequality (A.27) will be useful. Let ν = aλ -1 2 , so -ν is the largest eigenvalue of the matrix A in (A.15a). There exists a constant C A such that the transition matrices A n,k in (A.15c) satisfy [START_REF] Kailath | Linear systems[END_REF][START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF]

where |A n,k | Op denotes the Euclidean operator norm, equal to the largest singular value of the matrix A n,k .

Condition (A.17a) : to verify this condition, note that for arbitrary ε > 0,

) where the second inequality follows from (A.27). However, it follows from (u2) that there exists a uniform upper bound M w on the fourth-order moments of |w k | . Therefore, by Chebyshev's inequality [START_REF] Shiryaev | Probability[END_REF]

Since ν > 0, the right-hand side of (A.28b) has limit equal to 0 as n → ∞, by the Euler-Maclaurin formula [START_REF] Courant | Introduction to Calculus and Analysis[END_REF]. Replacing this limit from (A.28b) into (A.28a) immediately yields Condition (A.17a).

Condition (A.17b) : to verify this condition, recall that (w k ) is a sequence of squareintegrable martingale differences. Therefore, from (A.16)

where Σ k is the conditional covariance matrix in (A.18). Applying (A.27) to each term under the sum in (A.29a), it follows that

where d is the dimension of the parameter space Θ, and |Σ k | F denotes the Frobenius matrix norm. However, it follows from (u1) that there exists a uniform upper bound S on |Σ k | F . Therefore, by (A.29b)

Since ν > 0, the right-hand side of (A.29c) remains bounded as n → ∞, by the Euler-Maclaurin formula [START_REF] Courant | Introduction to Calculus and Analysis[END_REF]. This immediately yields Condition (A.17b).

Condition (A.17c) : to verify this condition, it is first admitted that the following limit is known to hold lim E (Σ k ) = Σ * (A.30a)

where Σ * was defined in (4.7). Then, let the sum in (A.17c) be written

Due to the equivalence A n,k ∼ exp(ln(n/k)A) (see [START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF], Page 125), the first term in (A.30b) is a Riemann sum for the integral [START_REF] Kailath | Linear systems[END_REF][START_REF] Nevelson | Stochastic approximation and recursive estimation[END_REF] a 2 1 0 e -ln(s) A Σ * e -ln(s)

which is known to be the solution Σ of Lyapunov's equation (4.11). The second term in (A.30b) can be shown to converge to zero in probability, using inequality (A.27) and the limit (A.30a), by a similar argument to the ones in the verification of Conditions (A.17a) and (A.17b). Then, Condition (A.17c) follows immediately.

Proof of (A.30a) : recall that

). Since (w k ) is a sequence of square-integrable martingale differences, it is possible to write, in the notation of (A.18),

Since U (θ, x) is a continuous vector field on Θ for each x ∈ X, and θ (k-1) converge to θ * almost surely, it follows that U (θ (k-1) , x) converge to U (θ * , x) for each x ∈ X, almost surely. On the other hand, it follows from (u2) that the functions under the expectation E θ * in (A.31b) have bounded second order moments, so they are uniformly integrable [START_REF] Shiryaev | Probability[END_REF]. Therefore, To obtain (A.30a) it is enough to note, as already stated in the verification of Condition (A.17b), that the Σ k are uniformly bounded in the Frobenius matrix norm. Thus, (A.31d) implies (A.30a), by the dominated convergence theorem.

A.4 Background on the Fisher information metric

Let D(θ) be the Kullback-Leibler divergence (4.2), or any other so-called α-divergence [START_REF] Amari | Methods of Information Geometry[END_REF]. Assume the Riemannian metric •, • of Θ coincides with the Fisher information metric of the model P . Then, for any local coordinates (τ α ; α = 1, . . . , d ), with origin at θ * , the following relation holds, by definition of the Fisher information metric (see [START_REF] Amari | Methods of Information Geometry[END_REF], Page 54),

where ∂ ∂τ α denote the coordinate vector fields of the local coordinates τ α . It is also possible to express (A.32) in terms of the Riemannian distance d(•, •), induced by the Fisher information metric •, • . Precisely,

Appendix B

Proof for chapter 5

B.1 Proof of Proposition 6

Let (θ (n) ) n 0 be a sequence generated by Algorithm 1. Recall the retraction Ret θ is defined in (3.18). Consider the sequence of tangent vectors U (θ (n) ) n 0 where U (θ (n) ) belongs to 1) , Σ (1) , β (0) ) and so on for U (θ (n) ) with n 0.

It is easy to see that either there exist infinitely many vectors U (θ (n) ) such that U (θ (n) ) = 0, or θ (n) = θ * for all n greater than some n 1 . Indeed, if there are only finitely many U (θ (n) ) such that U (θ (n) ) = 0, then there exists n 0 such that U (θ (n) ) = 0 for all n greater than n 0 . This implies there exists n 1 n 0 with U (θ (n 1 ) ) = -∇ * µ D(θ (n 1 ) ) = 0 U (θ (n 1 +1) ) = -∇ * Σ D(θ (n 1 +1) ) = 0 U (θ (n 1 +2) ) = -∇ * β D(θ (n 1 +2) ) = 0

From Algorithm 1, this implies that θ (n 1 ) = θ (n 1 +1) = θ (n 1 +2) and ∇ θ D(θ (n 1 ) ) = 0 (since all three components of ∇ θ D(θ (n 1 ) ) are zero). Since θ (n 1 ) ∈ Θ * , this implies θ (n 1 ) = θ * . Then (again by Algorithm 1) θ (n) = θ * for all n n 1 . Let θ (n i ) , U (θ (n i ) ) i 0 be the subsequence of θ (n) , U (θ (n) ) n 0 which consists of all the couples θ (n) , U (θ (n) ) such that U (θ (n) ) = 0 (i → n i is a function which counts the terms of this subsequence). Clearly, since U (θ (n) ) = 0 means that θ (n) will not be updated, it is enough to restrict attention to the subsequence θ (n i ) , U (θ (n i ) ) i 0 . For simplicity, this will be denoted θ (i) , U (θ (i) ) i 0 .

The subsequence θ (i) i 0 is given as in Algorithm 1 of [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], θ (i+1) = Ret θ (i) t i U (θ (i) ) with step-size t i chosen according to Armijo-Goldstein rule (precisely, t i = t n i where t n is given by t 0 = η µ , t 1 = η Σ , t 2 = η β , etc.). Moreover, (θ (i) ) i 0 remains within the compact neighborhood Θ * of θ * . According to Corollary 4.3.2 in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], if (η (i) ) i 0 is gradient-related then lim i→∞ ∇ θ D(θ (i) ) = 0.

Then, since θ * is the only stationary point of the cost function (4.6) in Θ * , it follows that lim i→∞ θ (i) = θ * , as required. To show that the subsequence U (θ (i) ) i 0 is gradient-Appendix C Proofs for Chapter 6

C.1 Proof of proposition 11

This proof is based on Remark 2 of [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a Riemannian manifold[END_REF]. Denote U θ (n) , x the direction of descent. According to this remark, if

almost surely, we have lim n→∞ θ (n) = θ * . In this situation, for each iteration, the direction of descent could be considered as