Member, IEEE Nora Ouzir

Senior member, IEEE Frédéric Pascal

Jean-Christophe Pesquet

Ieee Fellow

Convex Parameter Estimation of Perturbed Multivariate Generalized Gaussian Distributions

Keywords: Multivariate Generalized Gaussian Distribution, Multivariate Exponential Power Distributions, Convex Optimization, Robust Estimation, Precision Matrix Estimation, Elliptically Contoured Distributions

The multivariate generalized Gaussian distribution (MGGD), also known as the multivariate exponential power (MEP) distribution, is widely used in signal and image processing. However, estimating MGGD parameters, which is required in practical applications, still faces specific theoretical challenges. In particular, establishing convergence properties for the standard fixed-point approach when both the distribution mean and the scatter (or the precision) matrix are unknown is still an open problem. In robust estimation, imposing classical constraints on the precision matrix, such as sparsity, has been limited by the non-convexity of the resulting cost function. This paper tackles these issues from an optimization viewpoint by proposing a convex formulation with well-established convergence properties. We embed our analysis in a noisy scenario where robustness is induced by modelling multiplicative perturbations. The resulting framework is flexible as it combines a variety of regularizations for the precision matrix, the mean and model perturbations. This paper presents proof of the desired theoretical properties, specifies the conditions preserving these properties for different regularization choices and designs a general proximal primaldual optimization strategy. The experiments show a more accurate precision and covariance matrix estimation with similar performance for the mean vector parameter compared to Tyler's M -estimator. In a high-dimensional setting, the proposed method outperforms the classical GLASSO, one of its robust extensions, and the regularized Tyler's estimator.

I. INTRODUCTION

P PROBABILISTIC distribution parameters, particularly the mean and covariance -or scatter/precision matrix-, are central in statistical signal processing. In machine learning applications, they are building blocks of clustering, classification, dimension reduction, or detection models [START_REF] Roizman | A flexible EM-like clustering algorithm for noisy data[END_REF], [START_REF] Houdouin | Robust classification with flexible discriminant analysis in heterogeneous data[END_REF], [START_REF] Zhang | A novel m-estimator for robust pca[END_REF]. These parameters are rarely known in practice, and estimating them has been a long-standing statistical problem. A Maximum Likelihood (ML) approach relying on the standard Gaussian assumption has long been the solution of choice, bringing simplicity and tractability to performance analysis. The statistical properties of the Gaussian ML estimators are also well-known: the estimated mean is Gaussian, while the covariance matrix estimator follows a Wishart distribution.

The standard Gaussian assumption is no longer sufficient to handle the complexity of real-world datasets today. With datasets becoming increasingly large, high-dimensional, and heterogeneous, and often containing significant outliers, the limitations of the Gaussian probabilistic model are becoming more apparent. In domains such as radar [START_REF] De Maio | Modern radar detection theory[END_REF], financial signal [START_REF] Feng | A signal processing perspective on financial engineering[END_REF], and image processing [START_REF] Pascal | Parameter estimation for Multivariate Generalized Gaussian Distributions[END_REF], where these limitations are well-known, alternative models are being sought to address these challenges. Robust estimation theory is one of those alternatives that can provably deal with data perturbations [START_REF] Maronna | Robust M-estimators of multivariate location and scatter[END_REF]. It models observations with elliptically-contoured (EC) (complex and real Elliptically Symmetric (ES)) distributions. This broad-ranging model encompasses well-known multivariate distributions such as the Gaussian, Generalized Gaussian, t-, or k-distributions (see [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF] for a review).

Regularizing the parameter estimation problem has been another way of tackling current data challenges, such as high dimensionality. For example, sparse ℓ 1 -regularization of the precision matrix is a popular choice for Gaussian distributions [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. Works in the standard Gaussian framework have also laid the grounds for regularizing robust estimators of EC models. In [START_REF] Ollila | Regularized M-estimators of scatter matrix[END_REF], the scatter matrix has been regularized through the trace of its inverse. At the same time, closely related shrinkage Tyler's estimators have been studied in [START_REF] Pascal | Generalized robust shrinkage estimator and its application to stap detection problem[END_REF], [START_REF] Sun | Regularized Tyler's scatter estimator: Existence, uniqueness, and algorithms[END_REF]. Although these works have addressed the problem in the general statistical context of EC models, they remain limited from the regularization viewpoint. For example, unlike existing works in the Gaussian case, they have not considered ℓ 1 -regularization of the precision matrix due to the nonconvexity of the resulting problem in robust estimation. It is also worth noting that mostly plug-in (as opposed to modelbased) robust estimation strategies have been proposed in the Gaussian case with ℓ 1 -regularization [START_REF] Lafit | Robust and Sparse Estimation of Graphical Models Based on Multivariate Winsorization[END_REF], [START_REF] Tarr | Robust estimation of precision matrices under cellwise contamination[END_REF], [START_REF] Öllerer | Robust High-Dimensional Precision Matrix Estimation[END_REF]. In [START_REF] Wiesel | Geodesic convexity and covariance estimation[END_REF], a framework based on geodesic convexity has been proposed to regularize the perturbed covariance matrix of a zero-mean Gaussian distribution. Furthermore, previous works have not explored regularizing other parameters, such as the mean or data perturbations. As explained in the following, this work proposes an original model-based approach that copes with the non-convexity of the original problem and provides a general framework for exploring a larger variety of regularizations, along with associated theoretical conditions.

This work tackles the above issues for a popular subclass of EC models: the MGGD [START_REF] Pascal | Parameter estimation for Multivariate Generalized Gaussian Distributions[END_REF], also called MEP distribution [START_REF] Gómez | A multivariate generalization of the power exponential family of distributions[END_REF]. This distribution family is instrumental for image processing because of its ability to model a wide range of features (e.g., wavelet coefficients, gradients) [START_REF] Pascal | Parameter estimation for Multivariate Generalized Gaussian Distributions[END_REF], [START_REF] Kumar | Semi-supervised robust mixture models in rkhs for abnormality detection in medical images[END_REF]. It is also increasingly used to model weight distributions in (deep) neural networks [START_REF] Wang | Deep global generalized Gaussian networks[END_REF], [START_REF] Chai | Using generalized gaussian distributions to improve regression error modeling for deep learning-based speech enhancement[END_REF]. From a statistical point of view, its muchneeded flexibility is obtained thanks to an extra shape parameter β modelling lighter (β > 2) and heavier than the Gaussian (β < 2) tails. 1 The multivariate Gaussian distribution is a special case of MGGDs when β = 2. Furthermore, the parameter estimation problem is addressed in this work in the particular context of perturbed MGGDs. We are specifically interested in jointly estimating the mean vector, covariance matrix, and perturbation parameters, a complex problem introduced in [21] 2 , with broad applicability in real-world applications. Our main contribution is a convex reformulation of the resulting statistically robust joint parameter estimation problem. The proposed formulation offers all estimators theoretical guarantees of existence and establishes conditions for their uniqueness. In addition to modelling heterogeneity and outliers, the proposed formulation is flexible and can adapt to various reallife data structures. For example, the sparsity of the precision matrix can be incorporated through ℓ 1 -regularization to deal with high dimensional data with relatively few observations. Unlike previous approaches, the proposed framework enables regularizing all the involved parameters beyond the precision matrix. We also study the impact of these regularizations and specify the theoretical conditions they must satisfy to preserve convexity. Our second contribution is designing a proximal primal-dual algorithm tailored to solve the parameter estimation problem with guaranteed convergence [START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF], [START_REF] Combettes | Fixed point strategies in data science[END_REF]. The proposed algorithm is flexible, and associated proximity operators can easily incorporate various regularization choices.

The paper is organized as follows. Section II introduces the proposed perturbed MGGD model, starting with the underlying statistical framework. Section III presents the proposed convex parameter estimation approach, which is the main contribution of this paper. The proposed proximal primal-dual algorithm is detailed in Section IV. In Section V, we compare the performance of the proposed approach to different state-of-the-art robust and non-robust estimators in various experimental scenarios. Finally, concluding remarks and perspectives are provided in Section VI.

Notations: S K denotes the space of symmetric real matrices of size K × K, S + K is the cone of positive semi-definite matrices, and S ++ K the cone of positive definite matrices. Id denotes the identity matrix (whatever its size), 1 the vector (the dimension of which is understood from the context) whose components are all equal to 1, and tr(.) the trace of a matrix. • denotes the Euclidean norm, . S is the operator norm, and • r with r ∈ [1, +∞[ denotes the element-wise ℓ r norm (the same notation will be used for a matrix or for a vector whatever the dimension). Γ 0 (H) denotes the class of lower-semicontinuous convex functions from some Hilbert space H to ]-∞, +∞] which are proper (i.e., finite at least at one point). The domain of f is dom f = x ∈ H f (x) < +∞ . ι D denotes the indicator function of D ⊂ H, which is equal to 0 on this set and +∞ out of it. Finally, a function f

: H → ]-∞, +∞] is coercive if lim x →+∞ f (x) = +∞.

II. PROBLEM FORMULATION

The MGGD [START_REF] Gómez | A multivariate generalization of the power exponential family of distributions[END_REF] belongs to a broad subclass of EC distributions that can model various uni-modal probability density functions (p.d.f) with heavier or lighter tails than the Gaussian one. A random vector following an EC distribution can be defined using its stochastic representation [START_REF] Yao | A representation theorem and its applications to sphericallyinvariant random processes[END_REF]:

x d = µ + R A u, (1) where d
= stands for "is distributed as". The modular variate R is a positive random variable (with an unknown p.d.f), u is a random vector uniformly distributed on the unithypersphere {u ∈ R K | u = 1}, with R and u (statistically) independent. In (1), the mean vector µ and the scatter matrix factorization A (and the associated matrix C = A A ⊤ ) are the unknown parameters of the EC model. The stochastic representation described above is of practical interest for simulating EC-distributed random vectors. In particular, for the MGGD subclass, R β ∼ Γ(1/2, K/β), where Γ(a, b) is the univariate Gamma distribution with parameters a and b (see [START_REF] Norman L Johnson | Continuous univariate distributions[END_REF] for a definition).

The zero-mean MGGD of dimension K, denoted by MGGD K (β, 0, C), is characterized by its p.d.f

p x (•) = C K,β (det C) -1/2 exp - 1 2 (•) ⊤ C -1 (•) β/2 , (2)
where

C K,β = K Γ(K/2) π K/2 Γ(1+K/β)2 1+K/β is a constant, Γ(•)
denotes the Gamma function, and C ∈ S ++ K is (up to a multiplicative factor) the associated covariance matrix. 3 With this notation, when β = 2, (2) corresponds to the standard multivariate Gaussian distribution. When β < 2, distributions with heavier tails than the Gaussian one are obtained. In the remainder of this paper, we assume that the exponent β > 1 is known (see Section II-C).

A. Observation Model with Multiplicative Perturbations

This subsection introduces the proposed MGGD model under multiplicative perturbations. Let (x n ) 1 n N be N realizations of independent and identically distributed (i.i.d.) random K-dimensional vectors generated according to a zeromean MGGD with shape parameter β. Let us now consider a scenario where noisy observations (y n ) 1 n N result from a multiplicative perturbation τ of the unperturbed (x n ) 1 n N . The corresponding observation model reads

(∀n ∈ {1, . . . , N }) y n = τ n x n + µ, (3)
where 3) is equivalent to the assumption that the y n samples follow a MGGD K (β, µ, τ 2 n C). As (3) accounts for possible outliers in the observations (i.e., τ n values larger than 1), it can also be seen as a general EC distribution where τ n 's are realizations of an unknown positive p.d.f. Interestingly, similar models have been widely studied in the particular context of perturbed Gaussian distributions (see, e.g., [START_REF] Pascal | Covariance structure maximum-likelihood estimates in compound gaussian noise: Existence and algorithm analysis[END_REF]). Starting from (3), we aim to jointly estimate the noise and the unknown distribution parameters µ and C. Thus, N scalar parameters (τ n ) 1 n N need to be estimated in addition to the mean and covariance matrix. This problem has been addressed in [START_REF] Pascal | Covariance structure maximum-likelihood estimates in compound gaussian noise: Existence and algorithm analysis[END_REF] for the particular case of a centred Gaussian distribution (µ = 0 and β = 2.) The following subsection recalls the existing ML-based approaches in the more general case.

τ = (τ n ) 1 n N ∈ ]0, +∞[ N and β ∈]1, +∞[. The τ - perturbed model (

B. Previous Work on Estimating (C, µ, τ )

ML estimation of the unknown parameters (C, µ, τ ) has been broadly studied for EC models, including the MGGD or compound Gaussian models. (For non-perturbed EC models, one can refer to [START_REF] Maronna | Robust M-estimators of multivariate location and scatter[END_REF], [START_REF] Chitour | Exact maximum likelihood estimates for sirv covariance matrix: Existence and algorithm analysis[END_REF].) For the proposed τ -perturbed MGGD model, estimating C, µ, and τ can be achieved by minimizing the negative log-likelihood function arising from (3) (up to the normalizing constant that does not depend on the unknown parameters), i.e.,

L(C, µ, τ ) = 1 2 N n=1 (y n -µ) ⊤ C -1 (y n -µ) β/2 τ β n + N 2 log det C + K N n=1 log τ n . (4)
Function ( 4) is non-convex, and its minimization has been intensively investigated [START_REF] Chen | Robust shrinkage estimation of high-dimensional covariance matrices[END_REF], [START_REF] Pascal | Parameter estimation for Multivariate Generalized Gaussian Distributions[END_REF]: the standard approach being to first minimize with respect to τ , then plug the optimal value τ (C, µ) into (4). Precisely, the first step of the standard approach yields

τ (C, µ) = β[(y n -µ) ⊤ C -1 (y n -µ)] β/2 2K 1 β 1 n N , (5)
and by assuming none of the vectors (y n ) 1 n N is equal to µ (which is true with probability 1 for continuous random vectors for a given µ), plugging (5) into (4) leads to

L(C, µ, τ (C, µ)) = K 2 N n=1 log (y n -µ) ⊤ C -1 (y n -µ) + N 2 log det C, (6)
where the constant term KN β 1+log 2K β has been omitted. Note that ( 6) is (up to constant terms) the log-likelihood function of the (central) angular Gaussian distribution [START_REF] Tyler | Statistical analysis for the angular central Gaussian distribution on the sphere[END_REF]. For a given value of µ, minimizing [START_REF] Pascal | Parameter estimation for Multivariate Generalized Gaussian Distributions[END_REF] with respect to C (or C -1 ) on S ++ K thus leads to Tyler's estimator [START_REF] Tyler | A distribution-free M-estimator of multivariate scatter[END_REF], defined as the unique (up to a scale factor) solution to the fixed-point equation:

C = K N N n=1 (y n -µ)(y n -µ) ⊤ (y n -µ) ⊤ C -1 (y n -µ) . ( 7
)
Considering both C and µ unknown leads to the joint fixedpoint equations [START_REF] Frontera-Pons | Hyperspectral anomaly detectors using robust estimators[END_REF]:

                       µ = N n=1 y n (y n -µ) ⊤ C -1 (y n -µ) N i=1 1 (y n -µ) ⊤ C -1 (y n -µ) , C = K N N n=1 (y n -µ)(y n -µ) ⊤ (y n -µ) ⊤ C -1 (y n -µ) . (8)
Tyler's estimator has been extensively studied from a computational and statistical perspective ( [START_REF] Tyler | A distribution-free M-estimator of multivariate scatter[END_REF], [START_REF] Pascal | Covariance structure maximum-likelihood estimates in compound gaussian noise: Existence and algorithm analysis[END_REF], [START_REF] Pascal | Performance analysis of covariance matrix estimates in impulsive noise[END_REF]). Note that when one redefines Tyler's estimator as a Huber-type estimator [START_REF] Tyler | A distribution-free M-estimator of multivariate scatter[END_REF], an exponent 1/2 is added to the denominator in [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF]. The following remarks are of importance.

Remark II.1 On the joint fixed-point algorithm: (i) First, the existence, uniqueness, and convergence of the recursive algorithm associated with ( 8) has yet to be proved, although these estimators have been introduced in the 80s [START_REF] Tyler | A distribution-free M-estimator of multivariate scatter[END_REF]. Establishing convergence properties when both µ and C are unknown is still an open issue. If existence/uniqueness is assumed, then the asymptotic joint distribution of the corresponding estimators ( μ, Ĉ) is entirely characterized [START_REF] Tyler | A distribution-free M-estimator of multivariate scatter[END_REF]. (ii) The resulting estimators ( μ, Ĉ) do not depend on β. It can be shown (see [START_REF] Roizman | A flexible EM-like clustering algorithm for noisy data[END_REF]) that this is even true for any τperturbed EC distribution expressed by [START_REF] Zhang | A novel m-estimator for robust pca[END_REF], where extraparameters, such as shape, degree of freedom, or scale, only impact the estimation of τ .

In this work, we prove the existence, uniqueness, and convergence of all three µ, C, and τ estimators in a regularized setting.

C. Shape Parameter β

As mentioned previously, the shape parameter β only plays a role in estimating the perturbation τ . One could compensate for an approximate β value and retrieve the original distribution by varying τ , regardless of the true shape of the MGGD. Furthermore, estimating the (extra-)shape parameter of the EC model is a difficult problem that has been studied intensively in the non-perturbed MGGD case [START_REF] Sharifi | Estimation of shape parameter for generalized gaussian distributions in subband decompositions of video[END_REF], [START_REF] Song | A globally convergent and consistent method for estimating the shape parameter of a generalized gaussian distribution[END_REF], [START_REF] Dominguez-Molina | A practical procedure to estimate the shape parameter in the generalized gaussian distribution[END_REF], [START_REF] Krupiński | Approximated fast estimator for the shape parameter of generalized gaussian distribution[END_REF]. Other extra-parameters, such as the degree of freedom of the tdistribution [START_REF] Pascal | Improved estimation of the degree of freedom parameter of multivariate t-distribution[END_REF] or the shape parameter of a k-distribution [START_REF] Douglas | Reliable methods for estimating the k-distribution shape parameter[END_REF] have also been studied.

With this in mind, estimating β for the presented τperturbed MGGD model appears out of the scope of this work. We will assume that the shape parameter is known or previously estimated in a non-perturbed MGGD context for the experiments.

III. PROPOSED CONVEX FORMULATION

This section introduces a convex alternative to minimizing (4). The main idea of the proposed approach is to convexify the cost function through suitable variable changes and regularization. The purpose of regularizing the cost function is two-fold; first, prior information about the sought variables may be used to improve their estimation (e.g., sparsity of the precision matrix C -1 [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]). Secondly, regularization is critical in tackling the non-convexity of the cost function. The following subsection explains the details of these transformations.

A. A Regularized Cost Function

As explained above, convexifying (4) requires a series of transformations. First, let us re-parameterize L by setting

Q = C -1/2 (9) m = Qµ (10) θ = (θ n ) 1 n N = (τ β/(β-1) n ) 1 n N , (11)
where Q ∈ S ++ K . Using the re-parametrization ( 9)-( 11), we can rewrite the original cost function as

L(C, µ, τ ) = L(Q, m, θ), where L(Q, m, θ) = 1 2 N n=1 Qy n -m β θ β-1 n -N log det Q + K(1-1/β) N n=1 log θ n . ( 12
)
If we exclude the last term, the obtained function becomes convex with respect to (Q, m, θ). In Subsections III-B and III-E, we show that carefully choosing a regularization function for θ can counteract the concavity of the last term in [START_REF] Sun | Regularized Tyler's scatter estimator: Existence, uniqueness, and algorithms[END_REF] and lead to a fully convex formulation. In the following, we first introduce the general regularized form of ( 12) that accounts for prior information on the variables:

f (Q, m, θ) =      L(Q, m, θ, d)+g Q (Q)+g m (m)+g ϑ (θ) if θ ∈ ]0, +∞[ N and Q ∈ S ++ K +∞ otherwise, (13)
where g Q : S K → ]-∞, +∞], g m : R K → ]-∞, +∞], and g ϑ : R N → ]-∞, +∞] are regularization functions on Q, m, and θ, respectively. 4 From a Bayesian perspective, minimizing the regularized cost function [START_REF] Lafit | Robust and Sparse Estimation of Graphical Models Based on Multivariate Winsorization[END_REF] amounts to a Maximum A Posteriori (MAP) estimation of Q, m, and θ. The regularization functions introduced in ( 13) can be viewed as the potentials associated with (possibly improper) prior probability density functions proportional to exp(-g Q (•)), exp(-g m (•)), and exp(-g ϑ (•)), respectively. The regularized cost function f is suitable for different robust estimation problems depending on the choice of the regularization functions. For example, in graph processing applications, it is known that the precision matrix C -1 is sparse [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. Prior information on the nature of the perturbation τ may also be accessible (e.g., bounds), or one may seek to restrict the mean to a specific set (e.g., known values for a restricted subset of components). Examples of different practical regularizations will be discussed in more detail in Subsections III-D and III-E.

B. Study of Cost Function f

This subsection discusses key properties of the proposed cost function f . We will rely on these properties to choose the regularization functions and develop a suitable minimization strategy. Namely, we introduce two central propositions describing the assumptions that will ensure convexity and the existence of a minimizer. (Proof outlines are provided in Subsection III-C, and detailed proofs of these propositions can be found in Appendices A-B.) To simplify the analysis, we start by introducing the following functions:

ψ : R → ]-∞, +∞] : ξ → -log ξ if ξ > 0 +∞ otherwise, ( 14
)
g ϑ : R N → ]-∞, +∞] : θ → g ϑ (θ) -K(1 -1/β) N n=1 ψ(θ n ), (15)
Ψ : S K → ]-∞, +∞] : Q → -N log det Q if Q ∈ S ++ K +∞ otherwise. ( 16
)
We then have the following results.

Proposition III.1 Assume that g Q ∈ Γ 0 (S K ) and that there exists an invertible matrix

Q ∈ S ++ K such that g Q (Q) < +∞. Assume that g m ∈ Γ 0 (R K ) and g ϑ ∈ Γ 0 (R N ). Then, f is a proper lower-semicontinuous convex function on S K × R K × R N . Proof. See Appendix A.
Proposition III.2 In addition to the conditions stated in Proposition III.1, let us make the following assumptions:

(i) g m 0;

(ii)

g ϑ = g ϑ,0 + g ϑ,1 where g ϑ,0 ∈ Γ 0 (R N ), g ϑ,1 ∈ Γ 0 (R N ), dom g ϑ,1 = ]0, +∞[ N ,
and

(∀ θ ∈ [0, +∞[ N \ ]0, +∞[ N ) lim θ→ θ g ϑ,1 (θ) = +∞; (17)
(iii) Ψ + g Q and g ϑ are coercive functions. Then f admits a minimizer. Such a minimizer is unique if g Q and g ϑ are strictly convex.

Proof. See Appendix B.

It stems from Propositions III.1 and III.2 that the convexity of f and the existence of a minimizer depend on the choice of the regularizations. First, one can easily notice that the conditions on g m and g Q are satisfied by a wide range of typical regularization functions. Secondly, the assumptions on g ϑ are less usual but critical regarding the convexity of f . More precisely:

• The assumptions on g m are mild. They are satisfied both when no information is available on parameter m (g m = 0) or when this vector is known (g m = ι {m} with m ∈ R K ). In practice, restricting the mean to a specific set (e.g., known values for a restricted subset of components) satisfies these conditions. • The assumptions on g Q are rather classical. In particular, the core assumptions are satisfied by standard convex penalization promoting sparsity or group sparsity [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. The strict convexity assumptions can be satisfied by adding a quadratic term leading to an elastic net-like penalization [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. Subsection III-D will explicit some useful regularization functions for Q. • The key conditions for convexity are related to the choice of g ϑ . Note that this function contains the regularization that acts implicitly on the perturbation parameter τ . In Subsection III-E, we propose a regularization satisfying the required conditions for g ϑ and analyse its impact on estimating MGGD parameters.

C. Proof Outlines 1) Proof of Proposition III.1: The proof first shows that f is a sum of proper lower-semicontinuous convex functions, which makes it a lower-semicontinuous convex function. More precisely, because of the form of the first term in L, f involves the lower-semicontinuous envelope of a perspective function [START_REF] Combettes | Perspective functions: Proximal calculus and applications in high-dimensional statistics[END_REF], which leads to the first part of the function f being lower-semicontinuous and convex according to [START_REF] Bauschke | Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Space[END_REF]Proposition 9.42]. In addition, the assumptions made on the different regularizations in Proposition III.1 allow us to deduce that (the entire) function f is convex lower-semicontinuous. In a second step, based on the assumptions made on g Q and the form of the function g ϑ , we show that f is also proper (i.e., its domain is non-empty), which completes the proof.

2) Proof of Proposition III.2: The proof uses three main steps. First, we show the existence of a unique infimum m based on strict convexity and coercivity of f with respect to m (only). Then, using [START_REF] Bauschke | Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Space[END_REF]Proposition 8.35], we restrict the analysis to a new function composed of the marginal function and remaining parts of f independent of m. Writing the domain of this new function explicitly, we show that the necessary properties for [41, Proposition 9.33] are satisfied; thus, it is proper convex lower-semicontinuous. Based on Assumption (iii), we finally show that this new function is coercive. Thus, a minimizer exists and is unique if strict convexity is achieved for g Q and g ϑ .

D. Regularization on m and Q

As seen in Subsection III-B, the assumptions on the regularizations g m and g Q are relatively straightforward to satisfy. In the following, we let g m = ι {0} and consider a frequently encountered scenario where the precision matrix C -1 is sparse. Note that the standard Graphical LASSO (GLASSO) problem [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] can be obtained in the particular case where β = 2 by setting g m = ι {0} , g ϑ = ι {1} , and g Q : Q → λ Q 2 1 with λ ∈ ]0, +∞[. However, with the proposed re-parametrization, it is more direct to impose sparsity on the (new) variable Q by choosing g Q : Q → λ Q 1 . We will use the latter function for the experiments in Section V. Another typical choice for enforcing sparsity is elastic net regularization where

g Q = λ • 1 + ǫ 2 • 2
F with λ and ǫ ∈ ]0, +∞[. It should be noted that elastic net regularization is strictly convex, and if g ϑ is also strictly convex, f will admit a unique minimizer, as stated in Proposition III.2.

E. Regularization on θ

1) Choice of g ϑ : The choice of the regularization on θ is central to our analysis. The goal is to ensure the convexity of the global cost function f by compensating for the last term in [START_REF] Sun | Regularized Tyler's scatter estimator: Existence, uniqueness, and algorithms[END_REF]. The concavity of this term will be counterbalanced by choosing a regularization leading to a coercive function. A simple choice for the regularization function g ϑ allowing us to achieve this goal (while satisfying the requirements in Propositions III.1 and III.2) is the potential of a generalized Gamma distribution with scale parameter η, shape parameter (κ + 1), and exponent parameter α, i.e.,

(∀θ ∈ R N ) g ϑ (θ) = 1 η α θ α α + κ N n=1 ψ(θ n ) (18)
where ψ has been defined in [START_REF] Tarr | Robust estimation of precision matrices under cellwise contamination[END_REF],

α ∈ [1, +∞[, η ∈ ]0, +∞[, and κ ∈]K(1 -1/β), +∞[. We can then set g ϑ,0 (θ) = 1 η α θ α α g ϑ,1 (θ) = κ -K(1 -1/β) N n=1 ψ(θ n ),
which are functions satisfying Assumptions (ii) and (iii) in Proposition III.2. As emphasized above, coercivity is achieved by adjusting κ such that the weight κ -K(1 -1/β) of the logarithmic term is positive. In addition, since g ϑ,1 is then strictly convex, g ϑ is also strictly convex. It follows that the resulting cost function has a unique minimizer, provided that the chosen g Q is also strictly convex.

2) Impact on Cost Function f : To gain better insight into the impact of regularization [START_REF] Kumar | Semi-supervised robust mixture models in rkhs for abnormality detection in medical images[END_REF], we study the average behaviour of the cost function f with respect to θ and the associated regularization parameters η, κ, and α. 5 Let us first isolate the terms of interest by introducing functions

(f n ) 1 n N such that, for every (Q, m, θ) ∈ S ++ K × R K × ]0, +∞[ N , f (Q, m, θ) = N n=1 f n (θ n , Q, m) + Ψ(Q) + g Q (Q) + g m (m), (19)
and for every n ∈ {1, . . . , N },

f n (θ n , Q, m) = Qy n -m β 2θ β-1 n + θ α n η α -κ-K(1-1/β) log θ n . (20)
To focus our analysis on the variable θ, we make the simplifying assumption that Q and m have been perfectly identified with respect to the original statistical model. We can then restrict the analysis to the above function f n by looking at its average behaviour:

E{f n (θ n , Q, m)} = E{ Qy n -m β } 2θ β-1 n + θ α n η α - κ -K(1 -1/β) log θ n . (21)
In ( 21), (Qy nm) follows an MGGD K (β, 0, τ 2 n Id ), where τ n > 0 is the true value of the multiplicative perturbation factor in Model (3) (with re-parametrizations ( 9)-( 11)). We then have

E{ Qy n -m β } = 2K β θ β-1 n (22)
with

θ n = τ β/(β-1) n
, and

E{f n (θ n , Q, m)} = K(1 -1/β)f (θ n ), (23)
where

f (θ n ) = θ β-1 n (β -1)θ β-1 n + βθ α n K(β -1)η α -κ -1 log θ n ( 24
)
and κ = κβ/(K(β -1)).

With the conditions imposed on κ, one has κ > 1. The derivative of f at θ n is

f ′ (θ n ) = 1 θ n 1 - θ n θ n β-1 + αβθ α n K(β -1)η α -κ . (26)
Moreover, one has θ n = 1 for unperturbed data. In this case, if the scale parameter

η = αβ K(β -1)κ 1/α or equivalently, η = α κ 1/α , (27)
then f is decreasing over ]0, 1] and increasing over [1, +∞[ (see Appendix C). As expected, the minimum of this function is thus attained at θ n = 1. For a given value of κ, the choice for η in (27) makes function f independent of the vector dimension K, and we have, in this case,

f (θ n ) = θ β-1 n (β -1)θ β-1 n + κθ α n α -κ -1 log θ n . (28)
In Appendix C, we show that the above function has a unique minimizer θ n . If θ n > 1, θ n ∈]1, θ n [, which means that a bias is introduced by the regularization. As θ n is an increasing function of θ n , the method remains effective in reducing the influence of outliers. In addition, when θ n > 1, θ n is a decaying function of α (resp. κ). In other words, a stronger regularization comes at the expense of a larger bias in estimating the perturbation parameter τ n . This suggests choosing both α and κ close to 1.

Plots of the graph of the function f for different values of α are shown in Figure 1. As expected, when α increases, stronger perturbations θ n are penalized more heavily. Figure 2 shows the variations of θ n as a function of θ n for the same parameter values. The introduced bias can be seen in the increasing gap between the target θ n and θ n as α increases. We can also notice that the bias increases as the perturbations become larger for a given value of α. Note that although a larger value of α introduces more bias, it can also benefit the convergence speed of optimization algorithms. For example, when α = 2, f is a strongly convex function for which we can typically expect a linear convergence rate of optimization algorithms. This suggests that a trade-off between acceptable bias and reasonable computational speed is to be considered in practice.

IV. PROXIMAL PRIMAL-DUAL ALGORITHM

Minimizing cost function [START_REF] Lafit | Robust and Sparse Estimation of Graphical Models Based on Multivariate Winsorization[END_REF] requires handling various regularizations. Proximal algorithms can efficiently tackle such minimization problems: regularization terms are dealt with using proximity operators, usually in closed form [START_REF] Bauschke | Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Space[END_REF], [START_REF] Combettes | Perspective functions: Proximal calculus and applications in high-dimensional statistics[END_REF], [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF]. In this section, we will show that Problem (13) can be reformulated as the minimization of a sum of compositions of convex functions and linear operators, which can be solved using a proximal primal-dual algorithm. The primal space is here H = S K × R K × R N , endowed with the norm:

(∀p = (Q, m, θ) ∈ H) p H = Q 2 F + m 2 + θ 2 . (29) The dual space is G 1 × G 2 with G 1 = (R K ) N × R N and G 2 = S K × R N .
The first space is endowed with the norm

∀v 1 = (u n ) 1 n N , θ) ∈ G 1 v 1 G1 = N n=1 u n 2 + ω 1 θ 2 , (30)
while the second one is equipped with the norm

∀v 2 = Q, θ) ∈ G 2 v 2 G2 = Q 2 F + ω 2 θ 2 , ( 31
)
where

(ω 1 , ω 2 ) ∈ ]0, +∞[ 2 .
In the proposed algorithm, the latter parameters make the update of the θ estimates faster by introducing more flexibility. Note that all the described product spaces are finite-dimensional Hilbert spaces.

We reformulate the optimization problem as

minimize p∈H F (p) + G 1 (L 1 p) + G 2 (L 2 p), (32)
where L 1 and L 2 are linear operators such that

L 1 : H → G 1 : (Q, m, θ) → (T (Q, m), θ) L 2 : H → G 2 : (Q, m, θ) → (Q, θ), ( 33
)
and T is given by [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. With the choice of the regularisation functions made in Sections III-D and and III-E, F and G 2 are convex functions such that

F (Q, m, θ) = Ψ(Q) + g m (m) + g ϑ,1 (θ) G 2 (Q, θ) = g Q (Q) + g ϑ,0 (θ), (34)
and G 1 = Φ is given by [START_REF] Bauschke | Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Space[END_REF]. In addition, the existence of a solution to the minimization problem follows from Proposition III.2. Since F and G 2 are separable functions, their proximity operators can be calculated componentwise.

To solve [START_REF] Pascal | Performance analysis of covariance matrix estimates in impulsive noise[END_REF], a standard solution is the Chambolle-Pock algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], which reads

For k = 0, 1, . . .         p k+1 = prox γF (p k -γ(ζ 1 L * 1 v 1,k + ζ 2 L * 2 v 2,k )) p k = 2p k+1 -p k v 1,k+1 = Id -prox ζ -1 1 G1 (v 1,k + L 1 p k ) v 2,k+1 = Id -prox ζ -1 2 G2 (v 2,k + L 2 p k ), (35)
starting from initial values p 0 ∈ H, v 1,0 ∈ G 1 , and v 2,0 ∈ G 2 .

The variables (v i,k ) k 0 with i ∈ {1, 2} are dual variables. Hereabove L * 1 (resp. L * 2 ) is the adjoint of L 1 (resp. L 2 ), and (γ, ζ 1 , ζ 2 ) are positive scalar parameters. To ensure the algorithm convergence, these parameters must be chosen such that

γ(ζ 1 L 1 2 S + ζ 2 L 2 2 S ) < 1, ( 36
)
where • S is the operator norm. From the calculations of the norms of L 1 and L 2 in Appendix D, we deduce that a sufficient condition for [START_REF] Krupiński | Approximated fast estimator for the shape parameter of generalized gaussian distribution[END_REF] to be satisfied is

ζ 1 max{ Y 2 S , ω 1 } + ζ 2 max{1, ω 2 } < γ -1 , (37)
where

Y = y 1 . . . y N 1 . . . 1 . ( 38
)
We can then rewrite Algorithm [START_REF] Dominguez-Molina | A practical procedure to estimate the shape parameter in the generalized gaussian distribution[END_REF] to make the role of each involved variable more explicit. Let us set, at each iteration number k,

p k = (Q k , m k , θ k ), p k = ( Q k , m k , θ k ), v 1,k = ((u k,n ) 1 n N , θ # 1,k ), and v 2,k = (Q # k , θ # 2,k
), and let us define

ζ 3 = ω 1 ζ 1 and ζ 4 = ω 2 ζ 2 . θ # 1,k (resp. θ # 2,k ) corresponds to the dual variable of θ k in the space G 1 (resp. G 2 ), while Q # k stands for the dual variable of Q k in the space G 2 .
From the expressions of the adjoint operators of L 1 and L 2 derived in Appendix D, we obtain the following iterative algorithm:

For k = 0, 1, . . .                                  Q k = Q k -γ ζ 1 N n=1 u k,n y ⊤ n + ζ 2 Q # k Q k+1 = prox γΨ ( Q k + Q ⊤ k )/2 m k+1 = prox γgm m k + γζ 1 N n=1 u k,n θ k+1 = prox γ g ϑ,1 θ k -γ(ζ 3 θ # 1,k + ζ 4 θ # 2,k ) Q k = 2Q k+1 -Q k m k = 2m k+1 -m k θ k = 2θ k+1 -θ k For n = 1, . . . , N (u k+1,n , θ # 1,k+1,n ) = Id -prox ζ -1 1 ,ζ -1 3 ϕ (u k,n + Q k y n -m k , θ # 1,k,n + θ k,n ) Q # k+1 = Id -prox ζ -1 2 g Q (Q # k + Q k ) θ # 2,k+1 = Id -prox ζ -1 4 g ϑ,0 (θ # 2,k + θ k ).
The required proximity operators are provided in Appendix E.

V. EXPERIMENTS

This section evaluates the proposed method in comparison with various state-of-the-art approaches. For the proposed method, the experiments are carried out by imposing sparsity on the precision matrix and regularization [START_REF] Kumar | Semi-supervised robust mixture models in rkhs for abnormality detection in medical images[END_REF] on the perturbations (see Subsections III-D and III-E). We present results with varying N (sample size) and K (dimension) values. In particular, we study the standard case N > K and the high-dimensional setting where N ≪ K. The impact of sparse regularization on the precision matrix is investigated by considering dense and sparse matrices. Finally, we study the behaviour of the estimators with varying perturbation levels.

resulting in a tri-diagonal sparse matrix. The entries of the dense matrix C -1 2 are given by

(∀(i, j) ∈ {1, . . . , K} 2 ), C -1 2 (i, j) = ρ |i-j| , (40)
where ρ = 0.5 is the correlation coefficient used to generate both precision matrices. Fig. 3 shows the evolution of the MSE and consistency of the obtained mean, covariance, and precision matrices for different values of N in the case of a sparse precision matrix. As expected, the best mean estimates correspond to the empirical ones in the ideal case. The least accurate/consistent are those obtained with the empirical estimator in the perturbed case. Comparing the proposed method and Tyler's estimates reveals an almost identical performance. This is interesting because, unlike Tyler's method, the proposed method guarantees convergence in the unknown mean case. On the other hand, the results obtained by the proposed method for covariance and precision matrix estimation show a notable improvement over all other approaches. These improvements in MSE and consistency are even more noticeable for small sample sizes, indicating that the proposed method performs relatively well with a limited number of observations. The superior performance of the proposed method can be explained by the direct estimation of the precision matrix, strengthened by an appropriate sparse regularization. The benefits of sparse regularization in this first experimental scenario can be understood by studying the performance of dense precision matrix estimation (Table I). Here, the MSE gap between different methods is less pronounced, probably due to the minor impact of sparse regularization. However, it is worth noting that the performance of the proposed method remains superior to the other approaches for covariance and precision matrix estimation. As in the sparse case, the MSEs of the mean are similar for the proposed and Tyler's methods. We investigated the different estimator behaviours with varying perturbation levels and fixed K = 20 and N = 100. Fig. 4 shows the MSEs of the empirical (perturbed when p τ > 0), Tyler's, and proposed dense covariance matrix estimators 7 with different τ proportions. The experiments are conducted with β = 2 and β = 1.5, corresponding to a Gaussian and heavier-tailed distribution. Note that the empirical statistics are expected to perform ideally in the non-perturbed Gaussian case. As expected, the results show the robustness of the proposed method and Tyler's estimates compared to the empirical statistics when the perturbations increase. We can also see that the proposed method outperforms Tyler's estimates in all cases. Interestingly, the proposed method achieves the ideal performance of the empirical statistics in the non-perturbed Gaussian case (p τ = 0, β = 2). In the non-perturbed heavy-tailed scenario (p τ = 0, β = 1.5), the proposed method outperforms the empirical statistics. In contrast, Tyler's estimates are the least accurate when both distributions are not perturbed. D. High-Dimensional Case: N ≪ K For the higher dimensional experiments, we run N MC = 10 3 Monte Carlo runs for K = 100 and N -values between 10 and 80. We focus on an experimental scenario for which sparse regularization is suitable. A uniformly sparse precision matrix is generated with approximately 90% of the off-diagonal elements equal 0. The sparse regularization parameters providing the best MSEs are used for all methods, and we use the same values as before for all other parameters of the proposed method. Table II shows the obtained MSEs for the precision matrices estimated with the regularized Tyler, GLASSO, GlassoGaussQn, and proposed methods. Recall that GLASSO, its robust version GlassoGaussQn, and the proposed method employ a sparse ℓ 1 -regularization of the precision matrix. This could explain the gap in performance between these methods and the regularized Tyler's method. The performance of GLASSO is similar to its robust version GlassoGaussQn, although slightly improved for higher N -values. We can see that the proposed method results in the smallest errors for all experiments.

VI. CONCLUSION

This paper has introduced a novel robust and convex parameter estimation method for the perturbed multivariate Generalized Gaussian distribution. This optimization-based approach judiciously combines a reparametrization of the original likelihood with regularization. Unlike existing methods, theoretical convergence guarantees are obtained for all parameters, and experiments show improvements in performance in various simulation scenarios, including the high-dimensional case. Sparse precision matrix estimation is notably improved compared to regularized robust and non-robust plugin methods. These improvements can be useful for various (highdimensional) machine learning and graph processing applications. Forthcoming work will investigate such applications in the context of real noisy data.

APPENDICES

A. Proof of Proposition III.1

Let us define the following functions:

Φ : (R K ) N × R N : ((v n ) 1 n N , θ) → N n=1 ϕ(v n , θ n ), (41)
ϕ : R K × R : (u, ξ) →        u β 2ξ β-1 if u = 0 and ξ > 0 0 if u = 0 and ξ = 0 +∞ otherwise, (42)
and the linear operator

T : S K × R K → (R K ) N : (Q, m) → (Qy n -m) 1 n N . ( 43
)
Function ϕ is the lower-semicontinuous envelope of the perspective function of

(R K ) N × R N and the function (Q, m, θ) → Φ T (Q, m), θ ,
which corresponds to the composition of a linear operator with a proper lower-semicontinuous convex function is also proper, lower-semicontinuous, and convex. In addition,

(∀Q ∈ S K ) Ψ(Q) = N K k=1 ψ(σ k ), (44)
where (σ k ) 1 k K are the eigenvalues of Q. Ψ is a spectral function (i.e., a symmetric function of the eigenvalues of its matrix argument in S K ), which is associated with ψ ∈ Γ 0 (R). Therefore Ψ ∈ Γ 0 (S K ).

The considered cost function can be re-expressed as

f (Q, m, θ) = Φ T (Q, m), θ + Ψ(Q) + g Q (Q) + g m (m) + g ϑ (θ) (45)
and it follows from the previous observations on the two first functions Φ T •, • and Ψ, and the assumptions made on the three last ones that f is lower-semicontinuous and convex. Then, it is left to show that the domain of f is non-empty, i.e., f is proper. First, ∅ = dom g ϑ and, according to (15), g ϑ (θ) is not defined if θ ∈ ]0, +∞[ K , which means that dom g ϑ ⊂ ]0, +∞[ N . Since we have assumed that g Q is also finite at least for one symmetric positive definite matrix, f is proper. Thus, we have shown that f is a proper lower-semicontinuous convex function on

S K × R K × R N .
B. Proof of Proposition III.2

We have shown above that f is proper, lowersemicontinuous and convex. Now, we will show that a minimizer exists and provide the conditions for its uniqueness.

Proof. Since • β with β > 1 is a strictly convex function, for given values of (Q, θ) ∈ S ++ K × ]0, +∞[ N , m → Φ T (Q, m), θ + g m (m) (46)
is a strictly convex function. Assumption (i) yields

(∀m ∈ R K ) Φ T (Q, m), θ + g m (m) 1 2 N n=1 Qy n -m β θ β-1 n . (47)
This allows us to deduce from the coercivity of • β that the function in ( 46) is a coercive lower-semicontinuous strictly convex function. Its infimum is thus reached for a unique vector m(Q, θ).

Let function Θ be defined as

(∀(Q, θ) ∈ S K × R N ) Θ(Q, θ) = inf m∈R K Φ T (Q, m), θ + g m (m) + Ψ(Q) + g ϑ,1 (θ). (48)
It follows from [START_REF] Bauschke | Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Space[END_REF]Proposition 8.35] that the marginal function (Q, θ) → inf m∈R K Φ T (Q, m), θ + g m (m) is convex and from the convexity of the last two terms that Θ is a convex function. In addition, by incorporating the infimum m(Q, θ), and by using Definition [START_REF] Wiesel | Geodesic convexity and covariance estimation[END_REF] and Assumption (ii), for every

(Q, θ) ∈ S ++ K × ]0, +∞[ N we have Θ(Q, θ) = Φ T (Q, m(Q, θ)), θ + g m ( m(Q, θ)) + Ψ(Q) + g ϑ,1 (
The function

(Q, θ) → Θ(Q, θ) + g Q (Q) + g ϑ,0 (θ) (52)
is lower-semicontinuous and convex as a finite sum of lowersemicontinuous convex functions. It is also proper since, according to the assumptions in Proposition III.1 and Assumption (ii), there exists at least one point (Q, θ) ∈ S ++ K × ]0, +∞[ N such that g Q (Q), + g ϑ,0 (θ) is finite. In addition, for every (Q, θ) ∈ S K × ]0, +∞[ N ,

Θ(Q, θ) + g Q (Q) + g ϑ,0 (θ) Ψ(Q) + g Q (Q) + g ϑ (θ). ( 53
)
Using now Assumption (iii), we deduce that function (52) is coercive. Since it is a coercive proper lower-semicontinuous convex function, it admits a minimizer ( Q, θ). Hence, a minimizer of f is Q, m( Q, θ), θ). In addition, if g Q and g ϑ are strictly convex, function (52) is strictly convex, which ensures that it has a unique minimizer.

C. Study of Function f

The first and second-order derivatives of f in (28) read

(∀θ n ∈ ]0, +∞[) f ′ (θ n ) = 1 θ n 1 - θ n θ n β-1 + κ(θ α n -1) f ′′ (θ n ) = 1 θ 2 n β θ n θ n β-1 + (α -1)κθ α n + κ -1 . ( 54
)
The second derivative being positive, f is a strictly convex function. Since f (θ n ) → +∞ as θ n → 0 and θ n → +∞, f has a unique minimizer. This minimizer is a function of θ n , κ, α, and β but, for simplicity's sake, we will not make it explicit in our notation. The minimizer θ n satisfies the firstorder optimality condition

1 - θ n θ n β-1 + κ( θ α n -1) = 0, (55)
and we have

f ′ (1) = 1 -θ β-1 n (56) f ′ (θ n ) = κ θ n (θ α n -1). ( 57
)
We deduce that, if θ n > 1, then f ′ (1) < 0 and f ′ (θ n ) > 0, which shows that θ n ∈]1, θ n [. Similarly, if θ n < 1, then θ n ∈ ]θ n , 1[. In addition, by implicit derivation of (55),

q( θ n ) ∂ θ n ∂θ n = (β -1)θ β-2 n θ 1-β n (58) q( θ n ) ∂ θ n ∂α = -κ θ α n log θ n (59) q( θ n ) ∂ θ n ∂κ = 1 -θ α n (60) q( θ n ) ∂ θ n ∂β = ln θ n θ n θ n θ n β-1 (61)
where q( θ n ) = (β -1)θ

β-1 n θ -β n + κα θ α-1 n
. Since ∂ θ n /∂θ n > 0, θ n is an increasing function of θ n . If θ n > 1, then ∂ θ n /∂α < 0, ∂ θ n /∂κ < 0, and ∂ θ n /∂β > 0, which shows that θ n decreases w.r.t. α and κ, and increases w.r.t. β. Conversely, when θ n < 1, θ n is an increasing function w.r.t. α and κ, and a decreasing one with respect to β.

Note that, since θ n increases w.r.t. θ n , θ n → θ ∞ as θ n → +∞, where the limit θ ∞ is either finite or +∞. Eq. ( 55) allows us to discard the former case and to conclude that lim θn→+∞ θ n = +∞.

D. Properties of Operators

Let us now evaluate the norm of T . For every Q ∈ S K and m ∈ R K ,

T (Q, m) 2 = N n=1 Qy n -m 2 = N n=1 [Q -m] y n 1 2 = tr([Q -m]YY ⊤ [Q -m] ⊤ ), ( 64
)
where Y is given by [START_REF] Douglas | Reliable methods for estimating the k-distribution shape parameter[END_REF]. We also have

[Q -m]YY ⊤ [Q -m] ⊤ YY ⊤ S [Q -m][Q -m] ⊤ = Y 2 S [Q -m][Q -m] ⊤ . (65)
Combining (64) and (65) yields

T (Q, m) 2 Y 2 S tr([Q -m][Q -m] ⊤ ) = Y 2 S [Q -m] 2 F = Y 2 S ( Q 2 F + m 2 ). (66)
This shows that T S Y S . Using (62) and (63), we deduce that (37) is a sufficient condition for [START_REF] Krupiński | Approximated fast estimator for the shape parameter of generalized gaussian distribution[END_REF] to be satisfied.

Figure 1 .Figure 2 .

 12 Figure1. Graph of function f in[START_REF] Yao | A representation theorem and its applications to sphericallyinvariant random processes[END_REF] when η is given by (27), β = 1.7, and θn = 1. In blue, the non-regularized case when η → +∞ and κ = 0. Other plots correspond to κ = 1.1 with increasing values of α, i.e., α = 1 in red, α = 1.5 in yellow, and α = 2 in purple. The ordinate axis is in the log scale.

Figure 3 .

 3 Figure 3. MSEs and consistencies of different mean, covariance and precision matrix estimators for the experiments with a sparse AR precision matrix.

Figure 4 .

 4 Figure 4. MSEs of the estimated dense covariance matrices for varying perturbation proportions pτ for (a) an MGGD (β = 1.5) and (b) Gaussian distribution (β = 2).

  θ) < +∞.(49) This shows that the domain of Θ is the open set S ++ K × ]0, +∞[ N . It thus follows from[START_REF] Bauschke | Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Space[END_REF] Corollary 8.39] that Θ is continuous on this domain. In addition let( Q, θ) ∈ S + K × [0, +∞[N be a point on the border of dom Θ. One of the eigenvalues of Q or one of the components of θ is thus equal to zero. As a consequence of the form of Ψ in[START_REF] Wiesel | Geodesic convexity and covariance estimation[END_REF] and Assumption (ii),lim (Q,θ)→( Q, θ) (Q,θ)∈dom Θ Θ(Q, θ) lim (Q,θ)→( Q, θ) Q∈S ++ K ,θ∈]0,+∞[ N Ψ(Q) + g ϑ,1 (θ) = +∞.(50)It thus follows from[START_REF] Bauschke | Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Space[END_REF] Proposition 9.33] that Θ is a proper lower-semicontinuous convex function. Besides,inf Q∈SK ,m∈R K ,θ∈R N f (Q, m, θ) = infQ∈SK ,θ∈R N Θ(Q, θ) + g Q (Q) + g ϑ,0 (θ).

L 1 and L 2 L 1

 21 From the definitions of L 1 , L 2 , and the norms equipping the primal and dual spaces H, G 1 , and G 2 , we deduce that L 1 S = sup p∈H\{0}

Table I MSES

 I OBTAINED FOR THE EXPERIMENTS WITH A DENSE PRECISION MATRIX FOR INCREASING SAMPLE SIZES N .

	N	µ	Empirical C 2 C -1 2	µ	Tyler C 2 C -1 2	µ	Proposed C 2 C -1 2
	21 42.51 1.29e12 3.925 37.94 1.10e12 3.950 40.54 1.90e4 0.629
	30 33.23 7.03e3 0.732 29.37 8.05e3 0.791 29.98 1.65e3 0.330
	50 16.73 766.25 0.262 14.75 651.59 0.265 14.85 419.14 0.191
	80 11.43 235.91 0.127 9.99 177.13 0.118 10.04 147.80 0.100
	100 8.74 145.11 0.100 7.77 119.10 0.092 7.80 99.53 0.079
	200 4.40 47.14 0.045 4.07 38.55 0.040 4.06 35.51 0.037
	300 2.99 27.85 0.029 2.59 22.31 0.025 2.60 21.50 0.024
	400 2.23 20.14 0.022 1.98 15.49 0.019 1.98 15.25 0.018
	500 1.93 15.00 0.018 1.66 12.00 0.015 1.67 11.85 0.015
	1000 0.97 6.89 0.009 0.84 5.43 0.007 0.84 5.42 0.007

  Similar behaviour is observed for the precision matrix estimates, while the mean estimates follow previous findings and are comparable to Tyler's ones for varying perturbation levels.

			Table II
	MSES OBTAINED FOR SPARSE PRECISION MATRIX ESTIMATION WHEN
			K = 100 AND N ≪ K .
	N	Regularized GLASSO GlassoGaussQn Proposed Tyler
	10	0.80	0.77	0.76	0.64
	20	0.88	0.74	0.76	0.61
	40	0.65	0.48	0.48	0.42
	60	0.78	0.51	0.48	0.44
	80	0.89	0.57	0.53	0.44

7

  • β /2 [41, Example 9.43]. Thus, function Φ defined in (41) belongs to Γ 0

Note that β in this work corresponds to

2β in most MGGD-related literature. This choice allows us to simplify notations.[START_REF] Houdouin | Robust classification with flexible discriminant analysis in heterogeneous data[END_REF] A preliminary version of this work was presented in[START_REF] Ouzir | A convex formulation for the robust estimation of multivariate exponential power models[END_REF] with reduced theoretical results and limited experimental validation.

The covariance matrix of the MGGD is equal to 2 2/β Γ((K + 2)/β)(K Γ(K/β)) -1 C. For large values of K, the following approximation can be used (for computational purposes):Γ(a K + b) ∼ √ 2 π exp(-a K)(a K) a K+b-1/2 , leading to 2 2/β Γ((K + 2)/β)(K Γ(K/β)) -1 C ∼ 1/K (2 K/β) 2/β C.

These functions can take infinity values to model potential hard constraints on the variables. For example, if one seeks to restrict the vector m to some set D ⊂ R K (e.g., some hypercube or some ball) a suitable choice for gm is the indicator function ι D of D, equal to 0 on D and +∞ elsewhere.

We will denote the true target value in this subsection by θ (resp. τ ) to distinguish it from the parameter or variable θ (resp. τ ).

The observed samples (yn)1 n N follow an MGGD K (β, µ, τ 2 n C 1 ) (or an MGGD K (β, µ, τ 2 n C 2 )) according to (3).

The work of Jean-Christophe Pesquet was supported by BRIDGE-ABLE ANR Chair in AI. DGA has partially funded the research of Frédéric Pascal under grant ANR-17-ASTR-0015.

Our comparisons aim to evaluate the estimators' consistency and mean square error (MSE) with respect to the true parameters. The consistency of an estimated matrix  is quantified by  -A F , whereas the MSE is E  -A 2 F , where A is the true matrix. In our experiments, the expectation is computed empirically by averaging N MC = 10 4 Monte Carlo runs except for the high-dimensional case where we use N MC = 10 3 .

A. State-Of-The-Art Methods

In the standard case when N > K, the proposed method is compared to the empirical statistics and robust Tyler's estimates in both perturbed and ideal non-perturbed ((∀n) τ n = 1) scenarios. Note that a fixed-point update estimates the mean in Tyler's method. Although this approach has no convergence guarantees, it is frequently used in practice [START_REF] Frontera-Pons | Hyperspectral anomaly detectors using robust estimators[END_REF]. In the high-dimensional setting (N ≪ K), where Tyler's estimates do not exist, we show the results of the regularized Tyler's method [START_REF] Pascal | Generalized robust shrinkage estimator and its application to stap detection problem[END_REF], [START_REF] Ollila | Regularized M-estimators of scatter matrix[END_REF]. In addition, we compare the proposed method to classical GLASSO [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] (see Subsection III-D) and its robust extension proposed in [START_REF] Öllerer | Robust High-Dimensional Precision Matrix Estimation[END_REF]. Note that GLASSO and its robust version also use a sparse regularization on the precision matrix. In [START_REF] Öllerer | Robust High-Dimensional Precision Matrix Estimation[END_REF], several robust GLASSO extensions are proposed by replacing the classical GLASSO input with a robust covariance matrix estimate. Among the estimators presented in [START_REF] Öllerer | Robust High-Dimensional Precision Matrix Estimation[END_REF], we investigate GlassoGaussQn (the bestperforming one), which computes the robust covariance matrix using the Qn-estimator of scale and a robust correlation matrix.

B. Simulation Parameters

Unless indicated otherwise, β = 1.5, the maximum perturbation level is τ max = 5, and the proportion of corrupted data is set to p τ = 0.3 (when invariant). The mean vector is drawn from the standard normal distribution in all experiments. The parameters of the gamma prior are set to κ = 1.1 K(1 -1/β) and α = 1, which are parameters satisfying the conditions in Subsection III-E. The sparse regularization parameter λ is a user-defined parameter; in this work, it is adjusted automatically by providing the desired sparsity level (i.e., the proportion of zero entries in the generated matrices).

C. Experiments with N > K

The first experiments study two different cases of i.i.d. K = 20-dimensional data vectors (x n ) 1 n N distributed first according to an MGGD K (β, 0, C 1 ) where C -1 1 is a sparse precision matrix, and then, an MGGD K (β, 0, C 2 ) with a dense precision matrix C -1 2 . 6 The precision matrix C -1 1 is modelled by an auto-regressive (AR) process of order 3 such that

The adjoint of L 1 is

The adjoint of T can be deduced from the identity

which shows that

(70) In turn, the adjoint of L 2 is simply expressed as

E. Proximity Operators

The expressions of the proximity operators of the functions involved in (45), up to a positive scaling parameter γ, are provided below. 8 • Function ψ: For every ξ ∈ R,

• Function Ψ: For every Q ∈ S K , let us perform the eigenvalue decomposition of Q as U Diag(σ 1 , . . . , σ K )U ⊤ where (σ k ) 1 k K are the eigenvalues of Q and the columns of U form the associated orthonormal basis of eigenvectors. Then,

• Function Φ: As shown by [START_REF] Bauschke | Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Space[END_REF], Φ is a separable function of the components of its arguments. Thus,

The expression of the proximity operator of the perspective function ϕ has been derived in [START_REF] Combettes | Perspective functions: Proximal calculus and applications in high-dimensional statistics[END_REF]Example 3.7]. Let

8 See http://proximity-operator.net.

and, for every u ∈ R K \ {0} and ξ ∈ R such that

where t(u, ξ) is the unique solution in ]0, +∞[ of the equation

Then, for every u ∈ R K and ξ ∈ R,

As shown by the developments in Section IV (see ( 30) -( 31)), one may be interested in deriving the expression of the proximity operator in a metric that weights the two variables in an unbalanced manner. This means that, for every (γ

with γ = γ β 1 /γ β-1 2 .


domain: shs info

text:
(AI) Ethical issues are everywhere now: carbon footprint [START_REF] Strubell | Energy and policy considerations for deep learning in NLP[END_REF] Note: these numbers only concern 1 source of C02 out of 4 [START_REF] Bannour | Evaluating the carbon footprint of NLP methods: a survey and analysis of existing tools[END_REF] ⇒ largely under-estimated

Have you used them?

Guidelines, guidelines everywhere! [START_REF] Hagendorff | The ethics of ai ethics: An evaluation of guidelines[END_REF] Guidelines and checklists are great, but won't fix this "Currently, AI ethics is failing in many cases. Ethics lacks a reinforcement mechanism. Deviations from the various codes of ethics have no consequences. And in cases where ethics is integrated into institutions, it mainly serves as a marketing strategy. Furthermore, empirical experiments show that reading ethics guidelines has no significant influence on the decision-making of software developers." [START_REF] Hagendorff | The ethics of ai ethics: An evaluation of guidelines[END_REF] Beyond Guidelines

Examples of recommendations

x) Recommend that you formally clarify to the operator if the operator's role changes with regard to checking tread alignment, e.g. no longer has to check or checks less frequently. If the operators will no longer check, clarify this to the operators.

x.x-x) (Partner X) Regarding labeling of good and bad images for alarm: Recommend that you clarify who will label the images and estimate how many images need to be labelled and how long it will take.

Not published yet (submitted)

Guidelines, charters and good practices will (not) save us Ethics by design in an industry 4.0 project

Evaluating the Acceptability of Ethical Recommendations

Annotation process for the 120 recommendations 1. Author #1 categorized the recommendations → annotation guidelines 2. Author #2 read the guidelines and annotated 3. Inter-annotator agreement is computed:

→ disagreed on the categories for 60 recommendations (50%) → fully agreed in 38 cases (32%) → hesitated in 22 cases (18%) → strict observed agreement is 31.66%

Improving the guidelines:

  far from enough:" Neither the risk analysis informed by engineering practice, nor the socially informed engineering practice can be replaced by the other."[START_REF] Gurses | Engineering privacy by design[END_REF] AI-Proficient: an Industry 4.0 European project(2020)(2021)(2022)(2023) France, Belgium, and Germany AI researchers, heavy industry partners, technical partnersProject Ethics OfficerEthics by designThe AI-Proficient Approach: Industry is Time and Space uncover the spatial and temporal shop floor context what are the people (operators) doing already? where, when, how long? implicit workplace customs and unspoken responsibilities and hierarchies

  review definitions merge some categories: Human centering+Feedback, Training+Adaptation create new categorie: GDPR Participation from partners: results by category + partner Ethics team Partners "Evaluation" recommendations → the developers were probably worried that if they can't pass the benchmark, it means they can't use it → the "technical solution is always possible" notion is heavily embedded in the tech developer worldview Conclusions the partners had a range of interest in the ethics from 0 to lots there is no clear split between industrial and technical partners some types of recommendations seem easier to follow e.g. Human Centering: to be checked in a different context (robust?) → a large majority of the recommendations were at least partially implemented! Thank you for your attention! Questions ?

Ethical issues are everywhere now: water consumption Ethics by design in an industry 4.0 project Evaluating the Acceptability of Ethical Recommendations About Ethics Frameworks and Guidelines 80% of ethics documents are less than 6 years old [START_REF] Jobin | The global landscape of ai ethics guidelines[END_REF] AI Ethics Guidelines Global Inventory https://inventory.algorithmwatch.org: 167 AI Ethics Guidelines (as of 2020)

Ethics by design in an industry 4.0 project

Evaluating the Acceptability of Ethical Recommendations

The AI-Proficient Approach: Don't just leave it to Partners

The ethical issues appear through questioning Go beyond strictly AI ethical issues → other ethical issues are integrated

Participate in technical meetings

The AI-Proficient Approach: Be practical What have we achieved?

A flexible methodology [Anderson and Fort, 2022a] Quantitative results [this, submitted]

A New Scale Defining Human Involvement in Technology [Anderson and Fort, 2022b]