A Hippert-Ferrer
email: alexandre.hippert-ferrer@centralesupelec.fr

M N El Korso

A Breloy

G Ginolhac

Robust low-rank covariance matrix estimation with a general pattern of missing values

Keywords: Missing data, covariance matrix, mixture of scaled Gaussian, low-rank, expectation-maximization algorithm, classification

This paper tackles the problem of robust covariance matrix estimation when the data is incomplete. Classical statistical estimation methodologies are usually built upon the Gaussian assumption, whereas existing robust estimation ones assume unstructured signal models. The former can be inaccurate in real-world data sets in which heterogeneity causes heavy-tail distributions, while the latter does not profit from the usual low-rank structure of the signal. Taking advantage of both worlds, a covariance matrix estimation procedure is designed on a robust (mixture of scaled Gaussian) low-rank model by leveraging the observed-data likelihood function within an expectation-maximization algorithm. It is also designed to handle general pattern of missing values. The proposed procedure is first validated on simulated data sets. Then, its interest for classification and clustering applications is assessed on two real data sets with missing values, which include multispectral and hyperspectral time series.

Introduction

Missing data appear when no value of the data is available for a given variable and a given observation. This classical problem [START_REF] Anderson | Maximum likelihood estimates for a multivariate normal distribution when some observations are missing[END_REF][START_REF] Afifi | Missing observations in multivariate statistics I. Review of the literature[END_REF] is a pitfall in statistical signal processing and its related fields, as statistical inference [START_REF] Little | Statistical analysis with Missing Data[END_REF][START_REF] Schafer | Analysis of incomplete multivariate data[END_REF] and data analysis [START_REF] Benzécri | L'analyse des données[END_REF][START_REF] Van Buuren | Flexible imputation of missing data[END_REF]. To name a few applications where missing data has drawn significant attention, we can cite biomedical studies, chemometrics [START_REF] Walczak | Dealing with missing data : Part I[END_REF] or remote sensing where missing values created by poor atmospheric conditions or sensor failure can dramatically hamper the understanding of the physical phenomenon under observation [START_REF] Shen | Missing information reconstruction of remote sensing data: A technical review[END_REF]. Covariance matrix (CM) estimation theory, which is a fundamental issue in signal processing and machine learning problems, has witnessed particular efforts focusing in the case of incomplete data.

In this scope, it is known that any efficient estimation algorithm should be able to exploit the source of information offered by missing values, which is formerly called informative missingness [START_REF] Rubin | Inference and missing data[END_REF]. One approach to estimate the CM with missing values is to rely on maximum likehood (ML) estimation with a prior assumption on the probability distribution of the data. Within this framework, the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is a handy iterative procedure to obtain ML estimates as it is based on the expectation of the conditional probability p(z|x, θ) of the latent (missing) variables z given the observed variables x and the parameters θ under estimation. Extensive work has been put into CM estimation with missing values using the EM algorithm by assuming independent and identically distributed (iid) samples drawn from the Gaussian distribution when the sample size n exceeds the dimension p (n > p) [START_REF] Jamshidian | An EM algorithm for ML factor analysis with missing data[END_REF][START_REF] Liu | Efficient ML estimation of the multivariate normal distribution from incomplete data[END_REF] or in highdimensional regime [START_REF] Schneider | Analysis of incomplete climate data: Estimation of mean values and covariance matrices and imputation of missing values[END_REF][START_REF] Lounici | High-dimensional covariance matrix estimation with missing observations[END_REF][START_REF] Städler | Pattern alternating maximization algorithm for missing data in high-dimensional problems[END_REF] (p > n). Other models have been considered in which the covariance is assumed to have a low-rank plus identity structure [START_REF] Lounici | High-dimensional covariance matrix estimation with missing observations[END_REF], which will be referred in the following to as LR structure 1 . This structure is closely related to signal subspace inference or principal component analysis (PCA) with missing values [START_REF] Chen | Robust probabilistic PCA with missing data and contribution analysis for outlier detection[END_REF][START_REF] Josse | Handling missing values in exploratory multivariate data analysis methods[END_REF].

To overcome robutsness issues associated with the classical Gaussian assumption, more general distributions have been considered, such as the multivariate t-distribution [START_REF] Little | Robust estimation of the mean and covariance matrix from data with missing values[END_REF][START_REF] Liu | ML estimation of the t distribution using EM and its extensions[END_REF] and its regularized version for small sample size using an improved EM algorithm [START_REF] Liu | Regularized robust estimation of mean and covariance matrix for incomplete data[END_REF]. These distributions are encompassed by elliptical symmetric (ES) distributions, which are directly linked via ML estimation to the so-called M-estimators in the complete data case [START_REF] Tyler | A distribution-free M-estimator of multivariate scatter[END_REF][START_REF] Conte | Recursive estimation of the covariance matrix of a compound-Gaussian process and its application to adaptive CFAR detection[END_REF][START_REF] Pascal | Covariance structure maximumlikelihood estimates in compound Gaussian noise: Existence and algorithm analysis[END_REF]. The robustness of these tools to non-Gaussianity [START_REF] Zoubir | Robust statistics for signal processing[END_REF] has been illustrated in a wide range of applications including radar processing, hyperspectral imagery and classification [START_REF] Gini | Performance analysis of two adaptive radar detectors against non-gaussian real sea clutter data[END_REF][START_REF] Greco | Statistical analysis of high-resolution SAR ground clutter data[END_REF][START_REF] Theiler | Characterizing non-gaussian clutter and detecting weak gaseous plumes in hyperspectral imagery[END_REF][START_REF] Formont | Statistical classification for heterogeneous polarimetric SAR images[END_REF][START_REF] Ollier | Robust distributed calibration of radio interferometers with direction dependent distortions[END_REF][START_REF] Abdallah | Bayesian signal subspace estimation with compound Gaussian sources[END_REF]. Interestingly, this family can be accounted for in a robust (distribution-free) manner by considering the so-called mixture of scaled Gaussian (MSG) [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF] , which models the data as Gaussian conditionally to an monotone general random white=missing): monotone, general and random. Note that the random pattern is a special case of the general pattern, where missing values appear on individual observations rather than blocks. See [START_REF] Little | Statistical analysis with Missing Data[END_REF] (p.

5) for a precise description of missing data patterns.

unknown deterministic scale for each sample.

This paper proposes to take advantage of both robust estimation and LR structure models by proposing a EM-based procedure to estimate the CM in the case of incomplete data drawn from a MSG distribution. As cited above, existing robust covariance estimation algorithms based on MSG distributions are not designed to deal with missing values, and, to the best of our knowledge, existing works which use the LR structure are only based on Gaussian assumptions. In the full rank case, the work of [START_REF] Frahm | A generalization of Tyler's M-estimators to the case of incomplete data[END_REF] has extended generalized elliptical distributions (GES), from which MSG distributions are a sub-class, to incomplete data by proposing an adapted form of the Tyler's M-estimator.

Furthermore, existing estimation algorithms that handle missing data assume that the missingness pattern, i.e., the pattern describing which values are missing with respect to the observed data, is a monotone pattern [START_REF] Little | Robust estimation of the mean and covariance matrix from data with missing values[END_REF][START_REF] Liu | Efficient ML estimation of the multivariate normal distribution from incomplete data[END_REF][START_REF] Liu | Regularized robust estimation of mean and covariance matrix for incomplete data[END_REF][START_REF] Frahm | A generalization of Tyler's M-estimators to the case of incomplete data[END_REF] (see the illustration of missing data patterns in Fig. 1). This pattern can be of interest, e.g., in longitudinal studies [START_REF] Little | Statistical analysis with missing data[END_REF] or sensor failure [START_REF] Larsson | High-resolution direction finding: the missing data case[END_REF]. However, in other applications such as remote sensing, missing values can take very diverse patterns because of unpredictable events (clouds, snowfall, etc.) [START_REF] Shen | Missing information reconstruction of remote sensing data: A technical review[END_REF], which leads to the so-called general pattern [START_REF] Little | Statistical analysis with Missing Data[END_REF].

With the previous points in mind concerning i) robustness to non-Gaussianity, ii) low-rank structured models and iii) missingness patterns, the contributions of this paper are summarized below:

1) A generic algorithm for CM estimation is developed for structured signals with a non-Gaussian distribution. A procedure for non-structured / Gaussian is also obtained as a special case. In the structured configuration, the covariance matrix is supposed to have a LR structure as in [START_REF] Goodman | On clutter rank observed by arbitrary arrays[END_REF], which will be detailed in the next section;

2) The analysis of missingness patterns is extended to the general pattern, which fills a gap in robust estimation;

3) The proposed estimators are tested on both simulated and real data sets with two applications in machine learning, namely supervised classification and unsupervised clustering with missing values, which has rarely been unfolded (outside of simulated data) in the literature of CM estimation with missing values.

The rest of the paper is organized as follows: Section 2 formulates the problem by framing the MSG distribution, the LR CM structure and the missing data model.

Section 3 describes the proposed robust CM estimation procedures for unstructured and structured models. Section 4 illustrates the performance of the proposed method on simulated data sets in terms of CM estimation and data imputation. Finally, Section 5

shows the interest of the proposed procedure in real data applications using covariancebased machine learning methods.

Notations. a indicates a scalar quantity, a represents a vector quantity and A a matrix.

{a i } n i=1 denotes the set of elements a i with i ∈ [1, n].

Problem formulation

Mixture of scaled Gaussian distributions

Most of covariance matrix estimation procedures with incomplete data use the Gaussian assumption. However, this assumption can be inaccurate in remote sensing applications, where images often include heterogeneous areas. The mixture of scaled Gaussian

2 S p ++ = {Σ ∈ Sp : ∀x ∈ R p \{0}, x Σx > 0}. 3 Stp,r = {U ∈ R p×r : U U = Ir}.
(MSG) distributions can tackle this issue by managing heavier tails, which offers a better fit to empirical data [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]. A real p-vector y ∈ R p follows a zero-mean multivariate MSG distribution 4 if it admits the stochastic representation

y = √ τ n (1)
with n ∼ N (0, Σ) and for some scalar τ ∈ R >0 , called the texture, which is strictly positive, deterministic and unknown for each sample.

Let us now define a rectangular data set Y ∈ R p×n represented by Y = {y i = (y 1,i , y 2,i , . . . , y p,i ) } where {y i } n i=1 are modeled as n iid vectors of dimension p drawn from a MSG distribution. This leads to the following model:

y i |τ i ∼ N (0, τ i Σ), Σ ⊆ S p ++ , τ i > 0 (2)
The loglikelihood function of model ( 2) is given by

log ({y i }|Σ, {τ i }) ∝ -n log |Σ| -p n i=1 log τ i - n i=1 y i (τ i Σ) -1 y i (3)
In this model, the texture can be seen as a scale setting of the Gaussian model [START_REF] Wiesel | Unified framework to regularized covariance estimation in scaled Gaussian models[END_REF].

Such model is also intertwined with the class of compound Gaussian (CG) distributions (which is a subclass of ES distributions), which assumes the texture independent from n in (1) and with a given PDF f τ (.). Note that the type of distribution within the class of CG distributions is mainly guided by the assumption on the PDF function [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF].

Considering deterministic {τ i } instead of assuming a PDF f τ has proven its convenience in terms of robustness 5 [21, 23, 31]. This distribution is also more robust than the purely Gaussian one because the scales allow flexibility in the presence of heterogeneous data, e.g., noisy data, possible outliers or inconsistencies in the data (see, e.g., [START_REF] Mian | New robust statistics for change detection in time series of multivariate SAR images[END_REF] for Synthetic Aperture Radar data).

Remark 1. As it is clear that τ i = 1, ∀i in the deterministic MSG distribution [START_REF] Afifi | Missing observations in multivariate statistics I. Review of the literature[END_REF] gives the Gaussian distribution, results regarding CM estimation will be given by considering the Gaussian distribution as a special case of the MSG distribution.

Here we assume that Σ is characterized by a LR structure which can be modeled by the well-known factor model [START_REF] Ruppert | Statistics and data analysis for financial engineering[END_REF] (also known as spiked model [START_REF] Johnstone | On the distribution of the largest eigenvalue in principal components analysis[END_REF]):

           Σ = σ 2 I p + H H 0 rank(H) = r (4)
where I p denotes the p-dimensional identity matrix and H is a p × p low-rank signal covariance matrix of rank r. As in many works, the latter is considered to be known from prior physical assumptions [START_REF] Goodman | On clutter rank observed by arbitrary arrays[END_REF] or pre-estimated as in model order selection techniques [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF]. Model ( 4) is directly related to principal component analysis (PCA) and subspace recovery [START_REF] Tipping | Probabilistic principal component analysis[END_REF].

Data model

As each vector y i might have missing elements, it is necessary to design a model that takes into account incompleteness. Thereby, each of its observed and missing elements can be grouped in vectors denoted by y o i and y m i respectively, which are stacked in vectors y i such that

y i = P i y i =   y o i y m i   , i = 1, . . . , n (5)
where P i ∈ R p×p denotes a permutation matrix 6 . As illustrated by Fig. 2, this set of operations consists in permuting the missing elements at the bottom of each vector y i .

Then, the covariance matrix of y i becomes

Σ i =   Σ i,oo Σ i,mo Σ i,om Σ i,mm   = P i ΣP i (6)
where In a similar manner, the recent work of [START_REF] Aubry | Structured covariance matrix estimation with missing-(complex) data for radar applications via expectationmaximization[END_REF] models (5) by y i = A i y i , where A i is a n i × p selection matrix constructed from the extraction of the n i ≤ n rows of I p corresponding to the available observations at the i-th snapshot. This is only an alternative view of the problem. However, the derived results of [START_REF] Aubry | Structured covariance matrix estimation with missing-(complex) data for radar applications via expectationmaximization[END_REF] hold for the Gaussian distribution. In our case, it is worth mentioning that since each y i follows a MSG distribution, their permuted version y i also follow a MSG distribution because only elements within each p-vector are permuted.

Σ i,mm , Σ i,mo , Σ i,

The EM algorithm: a brief reminder

The EM algorithm is a widely employed iterative scheme for ML estimation in incomplete data problems [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. This algorithm offers a rigorous and formal approach to the intuitive ad hoc idea of filling in missing values7 : at the E-step, the following conditional expectation of the complete data loglikelihood is found given the observed data and the current deterministic estimated parameters θ:

L c ( Y |θ) = log ({y i }|Σ, {τ i }) ( 7
)
where θ is the vector of unknown parameters (which are τ i and Σ in our case) and

Y = [ y 1 , . . . , y n ]. Note that L c ( Y |θ) is equivalent to L c (Y |θ) since both Y and Y
contain the same missing values, ordered and unordered, respectively. Y (or Y ) are the so-called complete data which are the combination of observed and missing data.

At the M-step, the parameters are updated by maximizing the expected complete data loglikelihood. To summarize, if θ (t) is the current estimate of the parameter θ and f (•)

the density function, the E-step computes

Q(θ|θ (t) ) = L c ( Y |θ)f (y m |y o , θ = θ (t) )dy m (8)
and the M-step find θ (t+1) by maximizing (8):

Q(θ (t+1) |θ (t) ) ≥ Q(θ|θ (t) ) ( 9
)
The goal is then to repeat E and M-steps until a stopping criteria, such as the distance

||θ (m+1) -θ (m)
||, converges to a pre-defined threshold.

Covariance estimation under non-Gaussian environment

In this section, covariance matrix estimation procedures are developed in the presence of missing data under the MSG distribution. The rank constraint given by ( 4) is also resolved. For both configurations, the EM algorithm is adopted. In the following, we consider the determinant-based normalization [START_REF] Paindaveine | A canonical definition of shape[END_REF] which leads to the estimation of the normalized CM Σ/|Σ| 1 p , called the shape matrix [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF].

Robust full-rank estimation

In this case, the unknown parameters to estimate are θ = [ζ , {τ i } n i=1 ] , where ζ contains the elements of the lower triangular matrix of Σ including its diagonal and where {τ i } are the texture parameters. As we shall see hereafter, the differences with the purely Gaussian case (τ i = 1) are that 1) the unknown scales τ i must be taken into account in the formulation and computation of the expectation at the E-step and 2) a closed-form expression τ i must be derived to update Σ at the M-step. These estimators can be found by maximization the loglikelihood of the incomplete data, which is formulated onward.

Replacing y i by its permuted version y i in (3), the complete data loglikelihood of the MSG distribution is alternatively given by:

L c ( Y |θ) ∝ -n log |Σ| -p n i=1 log τ i - n i=1 y i τ i Σ i -1 y i (10)
At the E-step of the algorithm, the expectation of the complete loglikelihood ( 10) is computed by using the so-called Q-function. This function is the expectation of the missing data conditioned by the estimation of the parameters at the t-th iteration of the algorithm, that is:

Q i (θ|θ (t) ) = E y m i |y o i ,θ (t) L c (y o i , y m i |θ) (11)
Due to the iid of the observations, one obtains:

Q(θ|θ (t) ) = n i=1 Q i (θ|θ (t) ) (12)
where

Q i (θ|θ (t) ) = E y m i |y o i ,θ (t) -n log |Σ| -p log τ i -y i τ i Σ i -1 y i = E y m i |y o i ,θ (t) -n log |Σ| -p log τ i -   y o i y m i   τ i Σ i -1   y o i y m i   (13)
The computation of ( 13) is hastened as both first and second terms in the expectation are deterministic, which means that only the expectation of the last term has to be computed. To calculate this expectation, the trace tr{•} is used:

E y m i |y o i ,θ (t)   y o i y m i   τ i Σ i -1   y o i y m i   = E y m i |y o i ,θ (t) tr   y o i y m i   y o i y m i τ i Σ i -1 = τ -1 i tr E y m i |y o i ,θ (t)   y o i y o i y o i y m i y m i y o i y m i y m i   Σ -1 i = τ -1 i tr B (t) i Σ -1 i ( 14
)
where

B (t) i =   D (t) i E (t) i F (t) i G (t) i 
 is a p × p matrix at iteration t of the EM algorithm, with blocks given by

D (t) i = E y m i |y o i ,θ (t) y o i y o i = y o i y o i (15)
E (t) i = E y m i |y o i ,θ (t) y o i y m i = y o i E y m i |y o i ,θ (t) y m i (16)
F (t) i = E (t) i (17)
G (t) i = E y m i |y o i ,θ (t) y m i y m i (18)
Thus, the expectation to compute are

E y m i |y o i ,θ (t) [y m i ] and E y m i |y o i ,θ (t) [y m i y m i ],
which are the expectations of the missing data conditioned by the observed data. Note that y m i and y m i y m i are the sufficient statistics for the complete data Y , which have a distribution from the regular exponential family [START_REF] Little | Statistical analysis with Missing Data[END_REF]. Here, using a classical result on conditional distributions (see Theorem 2.5.1 in [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF], p. 35), the conditional distribution of y m i given

y o i is y m i |y o i ∼ N ( µ i,m|o , Σ i,mm|o ) (19)
where µ i,m|o = Σ i,mo Σ -1 i,oo y o i and Σ i,mm|o = τ i ( Σ i,mm -Σ i,mo Σ -1 i,oo Σ i,om ) are the conditional mean and covariance matrix, respectively. Note that τ i is absent from µ (t) i,m|o as it is annihilated by its inverse. Following [START_REF] Little | Statistical analysis with Missing Data[END_REF], at the t-th iteration, the E-step of the EM algorithm consists in calculating

E y m i |y o i ,θ (t) y m i = µ (t) i,m|o (20)
E y m i |y o i ,θ (t) y m i y m i = Σ (t) i,mm|o + µ (t) i,m|o µ (t) i,m|o (21)
Finally, θ (t+1) is obtained at the M-step of the algorithm as the solution of the following maximization problem:

max θ Q i (θ|θ (t) ) subject to Σ 0 τ i > 0, ∀i (22)
Proposition 3.1. The ML estimates Σ and τ i of problem [START_REF] Conte | Recursive estimation of the covariance matrix of a compound-Gaussian process and its application to adaptive CFAR detection[END_REF] are given by the following closed-form expressions:

τ i = tr B (t) i Σ-1 i p for i ∈ [1, n] (23)
Σ = p n n i=1 C (t) i tr C (t) i Σ -1 ∆ = H( Σ) (24)
with Σi = P i ΣP i , C

i = P i B (t)
i P i and where H(•) is the fixed point equation.

Proof. Let us start by writing the expectation of the loglikelihood L c ( Y |θ) given by ( 10):

E y m i |y o i ,θ (t) L c ( Y |θ) ∝ E y m i |y o i ,θ (t) -n log |Σ| -p n i=1 log τ i - n i=1 y i τ i Σ i -1 y i (25)
Notwithstanding equation ( 25) is similar to [START_REF] Schneider | Analysis of incomplete climate data: Estimation of mean values and covariance matrices and imputation of missing values[END_REF], it differs as it is the expectation conditioned on all observations contained in Y , whereas the latter stands for only one observation. By using the result ( 14) obtained from [START_REF] Schneider | Analysis of incomplete climate data: Estimation of mean values and covariance matrices and imputation of missing values[END_REF], it follows that:

E y m i |y o i ,θ (t) L c ( Y |θ) ∝ -n log |Σ| -p n i=1 log τ i - n i=1 τ -1 i tr B (t) i Σ -1 i (26)
Again, the block matrix

B (t)
i contains the expectation of the sufficient statistics y m i y m i to be computed at the E-step using [START_REF] Städler | Pattern alternating maximization algorithm for missing data in high-dimensional problems[END_REF]. Now, let differentiate [START_REF] Greco | Statistical analysis of high-resolution SAR ground clutter data[END_REF] with respect to τ i and resolve the equality:

∂E y m i |y o i ,θ (t) L c ( Y |θ) ∂τ i = 0
This calculus is trivial and finally gives

τ i = tr B (t) i Σ-1 i p
with Σi = P i ΣP i , which is the desired expression for τ i . To find Σ, we first replace τ i by τ i in [START_REF] Greco | Statistical analysis of high-resolution SAR ground clutter data[END_REF] and obtain the following loglikelihood:

E y m i |y o i ,θ (t) L c ( Y |θ) ∝ -n log |Σ| -p n i=1 log tr B (t) i Σ -1 i p -np (27)
Then, we derive (27) w.r.t. Σ, which leads to the resolution of the following equation:

∂E y m i |y o i ,θ (t) L c ( Y |θ) ∂Σ = 0
By using the definition Σ i = P i ΣP i , one obtains first:

-nΣ -1 -p n i=1 Σ -1 C (t) i Σ -1 tr C (t) i Σ -1 = 0 where C (t) i = P i B (t)
i P i . Arranging terms, then multiplying right and left terms by Σ gives the desired closed-form expression of Σ and concludes the proof.

Remark 2. For n > p, which is our case, this estimator can be computed using the fixed point algorithm Σ m+1 = H(Σ m ) where m refers to the iteration index of the fixed point.

Remark 3. In the Gaussian case, finding Σ that maximizes this expression requires to derive Q with respect to Σ and then to solve

∂Q ∂Σ = ∂E y m i |y o i ,Σ (t) Lc( Y |Σ) ∂Σ = 0 where E y m i |y o i ,Σ (t) L c ( Y |Σ) ∝ -n log |Σ| - n i=1 tr B (t) i Σ -1 i
. Based on the above and after some basic derivation calculus, the solution to be computed at each iteration of the EM algorithm is given by:

Σ = 1 n C (t) (28)
with

C (t) = n i=1 P i B (t)
i P i . One can notice that in the case of no missing values, this result leads to the classical Sample Covariance Matrix (SCM) with P i = I p and

B (t) i = E[y i y i ].
The complete estimation procedure of θ is given in Algorithm 1. The stopping condition of the EM algorithm is ensured by the evaluation of the quantity ||θ (t+1)θ (t) || 2 F at each iteration, while the convergence of the fixed point algorithm relies on

||Σ (t) m+1 - Σ (t) m || 2
F . At the step t = 0, the estimate Σ (0) is initialized with Tyler's estimator from available observations in their full dimension p, denoted Σ Tyl-obs . Unlike Σ, incomplete observations make the direct estimation of τ with the fixed point impracticable: thus, as an initialization, all τ (0) i are set to one.

Algorithm 1 EM-Tyl: Estimation of θ under MSG distribution with missing values.

Require: { y i } n i=1 ∼ N (0, τ i Σ), {P i } n i=1 Ensure: Σ, { τ i } n i=1
1: Initialization:

Σ (0) = Σ Tyl-obs τ (0) = 1 N 2: repeat EM loop, t varies 3: Compute E (t) i = y o i Σ (t) i,mo Σ -1(t) i,oo y o i G (t) i = τ (t) i Σ (t) i,mm -Σ (t) i,mo Σ -1(t) i,oo Σ (t) i,om + Σ (t) i,mo Σ -1(t) i,oo y o i Σ (t) i,mo Σ -1(t) i,oo y o i 4: Compute B (t) i =   y o i y o i E (t) i E (t) i G (t) i   5: Compute C (t) i = P i B (t) i P i 6:
repeat fixed point, m varies (optional loop)

7:

Σ (t) m+1 = H( Σ (t) m ) 8:
until

||Σ (t) m+1 -Σ (t) m || 2 F converges 9:
Compute τ

i , i = 1, . . . , n 10:

t ← t + 1 11: until ||θ (t+1) -θ (t) || 2 F converges

Robust low-rank estimation

Let us now consider the case of data whose covariance matrix lives in a subspace of dimension r < p. As we recall, this configuration is useful in many real signal applications and is closely related to principal component analysis (PCA) and signal subspace estimation [START_REF] Tipping | Probabilistic principal component analysis[END_REF].

Following model ( 4), the parameters to estimate are now given by θ

= [ζ , σ 2 , {τ i } n i=1 ]
, where ζ contains the elements of the lower triangular matrix of H. The maximization problem [START_REF] Conte | Recursive estimation of the covariance matrix of a compound-Gaussian process and its application to adaptive CFAR detection[END_REF] to find θ (t+1) at the M-step of the EM algorithm becomes the following low-rank estimation problem:

max θ Q i (θ|θ (t) ) subject to Σ = σ 2 I p + H rank(Σ) = r σ > 0, τ i > 0, ∀i (29)
A general solution to this problem was found in the seminal work of [START_REF] Tipping | Probabilistic principal component analysis[END_REF], which is recalled subsequently using the formulations of [START_REF] Kang | Rank-constrained maximum likelihood estimation of structured covariance matrices[END_REF] and [START_REF] Sun | Robust estimation of structured covariance matrix for heavytailed elliptical distributions[END_REF]. Firstly, let

Σ (t) EVD = p i=1 λ (t) i u (t) i u (t) i
be the eigenvalue decomposition of Σ (t) at the t-th iteration of the EM algorithm, where λ

(t) 1 < • • • < λ (t)
p are the eigenvalues of Σ (t) and u

(t) 1 , . . . , u (t)
p the corresponding eigenvectors. Then, the solution of ( 29) is given by:

Σ = σ2 I p + r i=1 λi u (t) i u (t) i (30) where σ2 = 1 p -r p i=r+1 λ (t) i (31) λi = λ (t) i -σ2 , i = 1, . . . , r (32)
The procedure to estimate θ under the LR assumption is given in Algorithm 2. It uses the same form than Algorithm 1 where equations ( 30)-( 32) are applied just after the fixed point algorithm to each newly estimated Σ. As stated in [START_REF] Sun | Robust estimation of structured covariance matrix for heavytailed elliptical distributions[END_REF], the proof of monotonicity of the low-rank estimation algorithm is ensured by standards convergence results of the majorization-minimization (MM) algorithm [START_REF] Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF]. Finally, the estimation of the scales τ does not change from the former full-rank algorithm.

Algorithm 2 EM-Tyl-r: low-rank estimation of θ under MSG distribution with missing values.

Require:

{ y i } n i=1 ∼ N (0, τ i Σ), {P i } n i=1 , rank r < p Ensure: Σ, { τ i } n i=1 1: Initialize Σ (0) , τ (0) as in Algorithm 1. 2: repeat EM loop, t varies 3: Compute E (t) i , G (t)
i , B

(t) i and C

(t) i as in Algorithm 1.

4:

repeat fixed point, m varies (optional loop)

5:

Σ (t) m+1 = H( Σ (t) m ) 6: Σ (t) m+1 EVD = p i=1 λ (t) i u (t) i u (t) i 7: Update Σ (t)
m+1 by computing ( 30), ( 31), [START_REF] Frahm | A generalization of Tyler's M-estimators to the case of incomplete data[END_REF] 8:

Σ (t) m+1 = Σ (t) m+1 /tr Σ (t) m+1 9: until ||Σ (t) m+1 -Σ (t) m || 2 F converges 10: Compute τ (t)
i , i = 1, . . . , n 11:

t ← t + 1 12: until ||θ (t+1) -θ (t) || 2 F converges

Implementation details

In both Algorithm 1 and Algorithm 2, the set of permutations matrices {P i } n i=1 is needed as an input to compute the transformed covariance matrices Σ i = P i ΣP i . Indeed, if a vector y i is fully observed, one finds that P i = I p , which makes the computation of Σ pointless. To avoid these extra computations, one can decompose the loglikelihood [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] in the following way:

L c ( Y |θ) ∝ -n log |Σ| -p n i=1 log τ i - i∈No y i τ i Σ i -1 y i - i∈Nm y i τ i Σ i -1 y i ( 33
)
where N o = {i l } no l=1 and N m = {i l } nm l=1 are the sets of indices corresponding to fully observed {y i } and partially observed { y i }, respectively, with n = n o + n m . Thenceforth, only the set of permutations matrices {P i l } nm l=1 need to be computed. An important consequence is that if y i is fully observed,

B (t) i = y o i y o i = y i y i . C (t) i
becomes:

C (t) i =      I p B (t) i I p = y i y i , i ∈ N o P i B (t) i P i , i ∈ N m (34)
Empirical experiments have shown that such care in the computation of C Furthermore, the fixed point loop is optional as it can be seen as an inner EM where {τ i } is the set of latent variables (hence a single update still increases the likelihood).

Our empirical experiments evidenced that both versions of the algorithm achieve similar performance, while performing only a single fixed point iteration tends to achieve a faster convergence.

Numerical simulations

This section illustrates the validation of the proposed algorithms with numerical experiments on simulated data drawn from the multivariate Gaussian and mixture of scaled Gaussian distributions. The performances of the proposed covariance estimation procedure are evaluated in regard to three aspects corresponding to different experiments:

1) The missing data ratio and patterns in subsection 4.1;

2) The quantity of outliers corrupting the data in subsection 4.2;

3) The possibility to perform CM-based data imputation in the presence of outliers and missing data in subsection 4.3.

For the sake of the experiment, the CM and scales parameters are known. The CM R is chosen to be Toeplitz8 , which has the form:

(R) ij = ρ |i-j| (35)
Parameter ρ, which controls the correlation structure of the CM, is set to 0.7. Scales {τ i } n i=1 are drawn from a Gamma distribution with shape parameter α and scale parameter β = 1 α . In all experiments, we fix α = 1. To generate a covariance matrix admitting the structure (4), we compute:

Σ = I p + σ 2 U U (36)
where U ∈ St p,r is the underlying signal subspace basis obtained from the eigenvalue decomposition of R and σ is a free parameter corresponding to the signal to noise ratio.

In our experiments, the data dimension is p = 15 and n = {63, 109, 190, 331, 575, 1000}.

Sets {y i } n i=1 are drawn from the MSG distribution with covariance Σ. As the aim is to estimate the structured covariance matrix Σ, the estimation is performed on 500 sets {y i } n i=1 simulated for each value of n. Indeed, for the sake of the experiment, missing data are also simulated, which allows a full control on their ratio and pattern. Importantly, as n increases from 63 to 1000, the missing data ratio decreases as {44, 22, 11, 5, 2, 1}%.

The following covariance estimators, which are reported in Table ??, are considered for comparison:

i) The covariance matrix for MSG distributions obtained from Algorithm 1 (full-rank) and from Algorithm 2 (low-rank), named Σ EM-Tyl and Σ EM-Tyl-r , respectively.

ii) The covariance matrix for the Gaussian distribution obtained from Algorithm 1 (full-rank) and from Algorithm 2 (low-rank) with the solution given by ( 28), named Σ EM-SCM and Σ EM-SCM-r , respectively.

iii) The sample covariance matrix (SCM) from the clairvoyant data (without missing data):

Σ SCM-clair = 1 n n i=1 y i y i ( 37
)
and its low-rank version, named Σ SCM-clair-r .

iv) The SCM Σ SCM-obs estimated from the set of fully observed vectors

{y i l } no l=1 ∈ N o where N o = {i l } no l=1 .
v) The fixed point estimator or Tyler's estimator from the clairvoyant data:

Σ Tyl-clair = p n n i=1 y i y i y i Σ -1 y i ( 38
)
and its low-rank version, named Σ Tyl-clair-r .

vi) Tyler's estimator Σ Tyl-obs from the set of fully observed vectors

{y i l } no l=1 ∈ N o with N o = {i l } no l=1 .
vii) Multiple imputation (MI) [START_REF] Royston | Multiple imputation of missing values[END_REF][START_REF] Van Buuren | Flexible imputation of missing data[END_REF] from which a robust version is proposed (RMI): for each y i with missing values y m i , q vectors with imputed missing entries (referred to as imputed vectors hereafter) ( y i1 , y i2 , . . . , y iq ) are generated with missing values drawn from a MSG distribution:

y m ij ∼ N (µ o ij , √ τ ij σ o ij ), j = 1, . . . , q (39)
where µ o ij and σ o ij are the mean and variance of the observed components of y ij , and τ ij are the scales parameters drawn from the Gamma distribution with shape parameter α = 1. The estimated covariance is the mean of the q Tyler's estimators from the q imputed vectors:

Σ RMI = 1 q q j=1 Σ Tyl,j = p nq q j=1 n i=1 y ij y ij y ij Σ -1 Tyl,j y ij (40)
The low-rank version is also estimated, named Σ RMI-r .

viii) Robust stochastic imputation (RSI): this procedure is tantamount to MI, but with q = 1.

viiii) Mean imputation [START_REF] Van Buuren | Flexible imputation of missing data[END_REF]: missing components of each vector y i are imputed by the mean of the observed components:

y m i = µ o i (41)
Then, the covariance Σ Mean-Tyl is estimated using Tyler's estimator. The low-rank version is also estimated, named Σ Mean-Tyl-r .

Type Estimator Description

Low In the experiments, low-rank estimators keep their names with an additional "-r", e.g., EM-SCM-r.

Covariance estimation with missing data

This experiment shows the performances of covariance estimation as functions of the missing data ratio and missing data patterns. To compare the estimated covariance matrix to the true data covariance, the following geodesic distance is used [START_REF] Bhatia | Positive definite matrices[END_REF]:

δ 2 S p ++ (Σ, Σ) = || log(Σ -1 2 ΣΣ -1 2 )|| 2 2 (42)
This distance, which is sufficient to measure estimation errors, emanates from the Fisher metric for the Gaussian distribution on S p ++ [START_REF] Bouchard | Riemannian geometry for compound Gaussian distributions: Application to recursive change detection[END_REF]. Three missing data patterns are examinated, which correspond to the following configurations (see Fig. 1):

• One block of missing data of size (7 × 20). This is a special case of the monotone missing data pattern as studied in [START_REF] Liu | Efficient ML estimation of the multivariate normal distribution from incomplete data[END_REF][START_REF] Liu | Regularized robust estimation of mean and covariance matrix for incomplete data[END_REF]. This configuration corresponds to the case where one group of variables is missing at the same observations, e.g., a group of sensors that are equally interrupted in time. This case is referred to as monotone pattern in the experiments. • Multiple blocks of missing data with various size and randomly scattered into the data set. This configuration corresponds to the general missing data pattern. It is most likely to happen in real-life applications (see, e.g., [START_REF] Musial | Technical note: Comparing the effectiveness of recent algorithms to fill and smooth incomplete and noisy time series[END_REF][START_REF] Shen | Missing information reconstruction of remote sensing data: A technical review[END_REF]). This case is referred to as general pattern.

• Randomly distributed missing data. This configuration is also a general missing data pattern except that multiple values for one variable at one observation are missing across the dataset. This case is referred to as random pattern. configurations. In the full rank case, EM-Tyl shows better performances than all other estimators, but fails when the number of available observations is small, especially when it gets close to p.

For the random pattern configuration, this result is even more notable and when n increases, EM-Tyl exhibits similar errors than mean imputation. Good performances of the latter can be explained by the increasing availability of observations as n increases, which gives a better estimate of µ o i . Unsurprisingly, all estimators (except the ones based on the SCM) reach Tyler's estimator on clairvoyant data as the missing data ratio decreases in large sample size. In the case of Gaussian estimators, the EM-SCM estimator shows very close performances than the clairvoyant SCM. As expected, Tyl-obs and SCM-obs display poor performances when the missing data ratio increases.

Results for the low-rank model are illustrated for r = 5. In all cases, we observe that EM-Tyl-r is a better estimator than EM-SCM-r in terms of distance. For small n, the gap between the EM-Tyl-r estimator and the imputation strategies (RMI-r and Mean-Tyl-r) diminishes compared to the full-rank case, and EM-Tyl-r significantly gives better estimates when the number of samples increases.

In conclusion, depending on the available sample size and the missing data pattern at hand, different estimators are preferable: imputation strategies are more advantageous for small amount of observed samples whereas our estimator performs better as the number of observed samples grows and when the covariance matrix has a low-rank structure.

Should discarded?

Data corrupted by outliers is one of the main motivation of robust estimation [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. In this case, one can be interested to know whether Tyler's estimator or EM-Tyl estimator gives higher performances. In the complete data case, the former can be directly estimated with outliers. To estimate the latter, outliers can be discarded or masked (i.e., set as missing data) if their position are known, which is a common approach in various applications. In the following, a data set {y i } n i=1 is drawn from the MSG distribution with a fixed number of samples n = 200 and with p = 15. Outliers are generated by adding white Gaussian noise (WGN) z i ∼ N (0, σ wgn I p ) with an increasing variance σ wgn on randomly selected observations y i .

Three data sets are considered to compare the CM estimation errors: 1) Corrupted data y i + z i ;

2) Masked data where the outliers are set as missing data.

3) Data without outliers y i ;

Tyler's estimator is computed on data set 1), whereas the EM-Tyl algorithm is applied to data set 2). Finally, the clairvoyant SCM and Tyler's estimators are computed from data set 3). Fig. 5 gives an overlook on what choice of the CM estimation procedure would be preferable. It shows that it mainly depends on the outlier variance: for a small σ wgn , Tyler's estimator should be preferred, whereas EM-Tyl is more suitable when the outlier variance reaches the variance of the signal σ. Note also that EM-Tyl provides better performances than the clairvoyant SCM until the missing data ratio reaches ∼45% of the data.

Robust imputation of missing values with outliers

In this experiment, we propose to apply our estimators to missing data imputation with possible outliers corrupting the data. Missing data imputation [START_REF] Van Buuren | Flexible imputation of missing data[END_REF] concerns a wide range of applications, including remote sensing [START_REF] Shen | Missing information reconstruction of remote sensing data: A technical review[END_REF]. A recent procedure to deal with missing values was developed, namely the EM-EOF method [START_REF] Hippert-Ferrer | EM-EOF: Gap-filling in incomplete SAR displacement time series[END_REF], which uses the EM algorithm and empirical orthogonal functions (EOFs) to iteratively decompose the CM and reconstruct the incomplete data with a few number of selected EOFs k p. The final number of components is the one giving the minimal error between the initial data and the imputed data. As shown in Fig. 6, instead of the SCM (which is used in the aforementioned study), the EM-Tyl-r estimator is used with r = k at the last iteration of the EM-EOF algorithm. This is justified by the fact that this algorithm iteratively updates the SCM with the predicted missing values at each step, whereas the EM-Tyl-r estimator only needs the missing data pattern, which remains the same at each iteration.

The RMI estimator [START_REF] Mian | New robust statistics for change detection in time series of multivariate SAR images[END_REF] is also considered for comparison. Quantity of outliers or gaps (%) Figure 6: Diagram of the EM-EOF method [START_REF] Hippert-Ferrer | EM-EOF: Gap-filling in incomplete SAR displacement time series[END_REF] for missing data imputation. After initializing the missing values, the CM matrix is estimated from the data matrix Y , which is then reconstructed with a number of components k p. Instead of using the SCM, the EM-Tyl-r estimator with r = k is plugged in at the last iteration of the algorithm (in red).

The data is generated using a Haystack-type model [START_REF] Lerman | Robust computation of linear models, or how to find a needle in a haystack[END_REF] which draws samples {y i } n i=1 as inliers y in i and outliers y out i :

{y i } n i=1 = {{y in i } i∈{1,...,nin} , {y out i } i∈{1,...,nout} } (43)
y in i ∼ N (0, I p + σ 2 s U U ) (44)
y out i ∼ N (0, I p + σ 2 o U ⊥ U ⊥ ) (45)
where U ⊂ St p,k is the signal subspace basis, σ 2 s and σ 2 o are respectively the signal to noise ratio (SNR) and the outlier to noise ratio (ONR), and 100 × (n-nin) n is the outlier ratio in the data set in percentage.

In the experiments, we choose σ 2 s = 10, σ 2 o ∈ {0, 1.5, 3, . . . , 30} and the outlier ratio varies from 0 to 50% of the data. Note that inliers and outliers are whole vectors (not just entries) which are randomly chosen in the data set. 30% of the data is discarded as missing values under a general pattern.

To measure the performance of imputation, a set Y = {y l } s l=1 consisting of s random points over {y ij } 1≤i≤n,1≤j≤p is selected. Y is called the cross-validation (CV) data set.

Then, we compute the root-mean-square error (RMSE) between Y and its reconstruction Y = {ŷ l } s l=1 after imputation:

δ imputation = 1 s s l=1 (ŷ l -y l ) 2 1/2 (46)
In total, 1% of the data is chosen for the cross-validation data set. These values are artificially removed and copied before the process, and compared to the new values after the imputation procedure. For model [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF], CV errors illustrated by Fig. 7 show a substantial gain by replacing the SCM by the EM-Tyl-r estimator at the last step of the EM-EOF method, whereas RMI performs poorly. The gain is larger for important signal to outlier ratio σ 2 o , which confirms the robustness of the proposed estimator for robust low-rank imputation under model [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF].

Application to classification and clustering problems

Classification with missing values

In remote sensing, and more particularly in multispectral and hyperspectral imaging, missing data can arise for various reasons including: 1) Clouds when the sensor operates in the visible part of the spectrum [START_REF] Shen | Missing information reconstruction of remote sensing data: A technical review[END_REF];

2) Sensor failure, as the known problem on the scan-line corrector images of the Landsat ETM+ sensor [START_REF] Chen | A simple and effective method for filling gaps in Landsat ETM+ SLC-off images[END_REF] or Aqua MODIS band 6 [START_REF] Rakwatin | Restoration of Aqua MODIS band 6 using histogram matching and local least squares fitting[END_REF];

In these cases, the data (or the entire band itself) is masked, downsampled to avoid temporal gaps or restored using specific gap-filling methods [START_REF] Chen | A simple and effective method for filling gaps in Landsat ETM+ SLC-off images[END_REF][START_REF] Rakwatin | Restoration of Aqua MODIS band 6 using histogram matching and local least squares fitting[END_REF]. When it comes to classification tasks, incomplete data is a challenge. In this scope, existing approaches classify, mask, and interpolate values from cloudy observations in a pre-processing step [START_REF] Rußwurm | Breizhcrops: A time series dataset for crop type mapping[END_REF].

To tackle this challenge, the proposed estimators are used as a set of descriptors 9{θ i } to a classification problem on the BreizhCrops data set [START_REF] Rußwurm | Breizhcrops: A time series dataset for crop type mapping[END_REF]. The data consists of Sentinel-2 multivariate time series of field crop reflectances on the Brittany region over 23 spectral bands. It is divided into four sub-regions called FRH01, FRH02, FRH03 and FRH04 corresponding to the four departments of Brittany. These sub-regions contain field labels which are gathered in 9 selected classes representing crop categories, e.g., barley, nuts, wheat, etc. Following [START_REF] Rußwurm | Breizhcrops: A time series dataset for crop type mapping[END_REF], each band is mean-aggregated over one field parcel to a feature vector y t ∈ R p , where p is the number of features (here the spectral bands) and t a timestamp. Thus, the whole data set is denoted {{y t k } K k=1 } T t=1 , which corresponds to the aggregation of all reflectances at field parcels k ∈

[1, K] and timestamps t ∈ [1, T ].
The experiment is done on the L1C top-of-atmosphere product and on 13 bands selected by the BreizhCrops processing chain [START_REF] Rußwurm | Breizhcrops: A time series dataset for crop type mapping[END_REF]. For the L1C product, the total number of parcels is K = 608263 and the number of observations per parcel is T = 51.

A supervised covariance-based classification is performed using a Minimum distance to Riemannian mean (MDRM) classifier [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF], which works in a test-training form. The training step provides a set of p × p SPD matrices encoding field parcels for the available classes. Assuming that the data is complete, the set of SCMs would be computed as:

Σ k = 1 T T t=1 y t k (y t k ) , k = 1, . . . , K (47)
For each class, a center of mass of the available parcels is estimated. In the test step, a field parcel is also encoded as an SPD matrix and then assigned to the class whose center of mass is the closest according to the distance (42) acting on the manifold S p ++ [START_REF] Congedo | The Riemannian minimum distance to means field classifier[END_REF]. Before the experiment, missing data are artificially created in the data set. For this, successive bands are set as missing and the performances are evaluated as function of the missing data ratio. The classifier is trained first on the FRH01 region and tested on the FRH03 region. A second experiment is conducted, where the train data set is changed to the FRH04 region while the test set remains unchanged. observations in each window denoted X i ∈ R p×n . Thus, we obtain a set of descriptors {θ i } to cluster using a K-means++ [START_REF] Vassilvitskii | k-means++: The advantages of careful seeding[END_REF]. For the descriptors using the low-rank model ( 4 The averaged Overall Accuracy (OA) are reported in Fig. 10. We observe that in general, the descriptors based on RSI and RSI-r estimators give lower accuracies compared to the descriptors using the SCM, SCM-r, EM-Tyl and EM-Tyl-r estimators. As the missing data ratio increases, the performances of the EM-SCM estimator are undermined whereas EM-Tyl estimator accuracies remain stable around 50%. For large missing data ratio, EM-SCM-r performs surprisingly well compared to EM-Tyl and EM-Tyl-r. For lower missing data ratios, clustering achieves its best performances using the EM-Tyl-r estimator.

Conclusion

This article proposes a novel procedure based on the EM algorithm to perform ro-

Figure 1 :

 1 Figure 1: Illustration of three rectangular data sets with different missing data patterns (black=observed,

  The transpose operator is , whereas tr{•} and | • | are respectively the trace and the determinant operators. The eigenvalue decomposition of A into matrices U and Λ is denoted by A EVD = U ΛU . The symbol indicates positive definitiveness, S p ++ is the set of p × p symmetric positive definite (SPD) matrices 2 and St p,r is the real Stiefel manifold of p × r orthogonal matrices 3 . ∝ stands for "proportional to". Q(•) is the score function. Finally, E[•] denotes the expectation operator.

  N o or i ∈ N m accelerates the running time of the EM algorithm.

Figure 3 :

 3 Figure 3: Mean of error measures 42 for r = p and r = (right) of methods EM-Tyl, EM-Tyl-r, EM-SCM, EM-SCM-r, Tyl-clair, Tyl-obs, SCM-clair, SCM-obs, RMI and mean imputation as functions of the number of samples n. The mean are computed over 500 simulated sets {y i } n i=1 for monotone (top), general (middle) and random (bottom) missing data pattern with p = 15.

Fig. 3 20 T 2 Figure 4 :

 32024 Fig.3 and 4show the results in terms of mean distance for two rank settings (the full rank case r = p and the low-rank case with r = 5) and the three aforementioned pattern 20

Figure 5 :

 5 Figure 5: Mean geodesic distances over 200 simulated sets {y i } n i=1 as function of the missing data or outlier ratio (%) for n = 200. Tyl-clair and SCM-clair errors are shown with no outliers. Tyl-clair is displayed on corrupted data with different outlier variance σwgn w.r.t. the data variance σ. EM-Tyl and Tyl-obs errors are shown when outliers are masked as missing values.

Figure 7 :

 7 Figure 7: Mean of cross-validation error (46) over 200 simulated sets {y i } n i=1 as function of the outlier ratio (%) and signal to outlier ratio σ 2 o for the EM-EOF using EM-Tyl-r estimator, EM-EOF using the SCM and RMI procedures. 30 % of missing data under the general pattern are generated. Chosen parameters are n = 200 and p = 15.

1 Figure 8 :

 18 Figure 8: Classification mean overall accuracy (%) versus the number of randomly missing bands (1 band ∼ 7% of the data) over 50 run of the MDRM [62] classifier based on EM-SCM, EM-Tyl-τ i and RSI estimators.

  Fig.9a). As shown by Fig.9b, sensor failure is simulated by adding missing values on random columns of selected bands.

  bust low-rank estimation of the covariance matrix with missing data following a general pattern. The developed tools take advantage of the MSG distribution and of the LR structure of the covariance matrix. Closed-form expressions of the unknown CM and associated scales are derived and integrated to the M-step of the algorithm, which generalizes the Gaussian case. The performance of the proposed estimators are validated on simulated data sets with missing values and corrupted by outliers, as well as realworld incomplete data sets. Compared to the classical Gaussian assumption and to the unstructured model, experiments show the possibility to improve the results in terms of CM estimation and robust imputation, as well as supervised (classification) and unsupervised (image clustering) problems. Some potential extensions of this work include the extension to other classes of compound Gaussian distributions and experimenting classification tasks with temporal gaps rather than spectral gaps, as well as more spectral information by reducing the downsampling rate.

  Missing elements (colored) of each y i are permuted to fall at the bottom of each y i . Notice that a completely observed y i implies P i = Ip.

	y1 y2	yN	y1 y2	yN
	y11		y11
	y12	y i = y i Pi	y13
	y13		y14
	y14		y15
	y15		y12
	Figure 2:
		oo are the block CM of y m i , of y m i and y o i , and of y o i .

Table 1 :

 1 List of estimators used in numerical experiments.

	-rank

It is sometimes named spiked or factor model.

Note that this distribution can be easily transposed to complex-valued vectors[START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF].

This model follows a parametrization of the covariance matrix of the real elliptical model E(0, τ ξ), where ξ is the shape matrix and where τ has only one value. In our case, the scale parameter varies at each observation i = 1, . . . , n. More details can be found in[START_REF] Yao | A representation theorem and its applications to spherically-invariant random processes[END_REF][START_REF] Zoubir | Robust statistics for signal processing[END_REF][START_REF] Maronna | Robust statistics: theory and methods (with R)[END_REF].

Note that P i is invertible and P -1 i = P i .

Indeed, it is not the missing values themselves that are filled, but the function of the missing values (the sufficient statistics) that are computed[START_REF] Little | Statistical analysis with Missing Data[END_REF].

Note that this information on the structure is not used in the estimation procedure.

In machine learning problems, statistical descriptors are classically used as they are often more discriminative than raw data.

This measure can be easily computed by the formula r i=1 λ i / p i=1 λ i .

Acknowledgment

This work was supported by ANR ASTRID project MARGARITA (ANR-17-ASTR-0015).

The authors would like to thank Antoine Collas and Ammar Mian for their precious help on the implementation pipeline of the MDRM and Kmeans + + algorithms for real data, and Florent Bouchard for his deep insights in robust estimation and covariance-based classification.

Results in terms of overall accuracy (OA) for the FRH03 region with two different training sets versus the number of missing bands are displayed in Fig. 8. Three covariance estimation strategies are shown: covariance estimation for Gaussian data (EM-SCM), robust stochastic imputation (RSI) and robust covariance estimator (EM-Tyl) multiplied by the geometric mean of the estimated scales τ i :

Here, the geometric mean acts as an additional power information to the EM-Tyl estimator and better fits the average of textures that might vary on different scales. Unfortunately, experiments have shown that the low-rank structure does not improve the OA, which might be due to the mean-aggregation over each parcel, which essentially acts as a filter.

Results show that, for this data set, the MDRM classifier based on robust covariance matrix estimation (EM-Tyl-Πτ i -based) is generally more robust to missing bands than the EM-SCM one, whereas the RSI-based classifier gives a rapidly decreasing OA when the gaps ratio increases. Furthermore, we observe that the EM-SCM-based classifier is more suited for large missing data ratio when the training phase is performed on the FRH04 region. Interestingly, the classifier based on the EM-Tyl-Πτ i estimator provides an almost stable OA until the gap ratio reaches 28% of the data (4 missing bands) when when FRH01 is the training set and 14% of the data (2 missing bands) when the classifier is trained on FRH04.

Image clustering with missing values

We apply the proposed CM estimators to a hyperspectral image clustering problem using the K-means++ algorithm [START_REF] Vassilvitskii | k-means++: The advantages of careful seeding[END_REF] on the Indian Pines data set [START_REF] Baumgardner | band aviris hyperspectral image data set[END_REF], consisting of a 145 × 145 pixels image with p = 200 bands. As in the classification problem, the proposed estimators EM-SCM and EM-Tyl are used as descriptors {θ i }, as well as EM-SCM-r and EM-Tyl-r. The aim is to partition the descriptors in a number of clusters which correspond to the 16 classes dividing the Indian Pines image.

The image is first centered by subtracting the global mean. Then, a sliding window of size w × w is applied to the image. One descriptor {θ i } is estimated using the n = w 2


Jean-François Silvain
email: jean-francois.silvain@icmcb.cnrs.fr

Helies Boumali

Julie Bourret

Pierre-Marie Geffroy

Sébastien Fourcade

Thomas Batigne

Karim Sinno

Yongfeng Lu

Loic Constantin

Sebastien Fourcade

Porous 3D Cu structures with adaptive heat dissipation properties

Keywords: 3D printing, Cu, debinding, sintering, Additive manufacturing, thermal conductivity

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

With the emerging of new technologies, the field of transport tends to rely more on electrical systems, such as auxiliary controls as well as traction requirements. The aerospace industry is turning to the manufacture of more electric planes, which allow better management of the rationalization / management of energy on board.

The increasing share of electricity in airplanes reduces fuel consumption, cuts costs and maintenance time. As a result, the electrical power used on board in transport systems is increasingly important. During operation, electronic systems undergo thermal gradients, which generate thermomechanical stresses at the interfaces between these different layers of electronic systems. These interfacial stresses, induced by the differences in thermal expansion coefficient between the different layers, then produce thermal fatigue, which leads to irreversible physical degradation of the system. The on-board electronic components are also temperature sensitive and have temperature limits for reliable operations (e.g., ~140 °C maximum for silicon). The heat dissipation in these components is essential to ensure an operating temperature below this maximum temperature. Thus, cooling systems have been then developed by several teams around the world to increase the heat dissipation performance through studying new materials and geometries [1][2][3].

Constantin et. al [START_REF] Constantin | Laser 3D printing of complex Copper structures[END_REF] has shown that increasing the specific surface area of a material improves the heat dissipation properties of the same material. So, two strategies can be envisioned through: 1) producing complex shapes in order to obtain parts with a large specific surface area, and 2) developing materials intrinsically having larger specific surface areas (i.e., porous materials).

The high thermal conductivity (~ 400 W/mK) of pure copper (Cu) makes it an essential material for thermal management in numerous industries, such as microelectronics, power plants, and aerospace [START_REF] Singer | 16th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm)[END_REF][START_REF] Körner | [END_REF][7]. However, the shape (i.e., surface area) of heat sinks, playing a crucial role in heat dissipation, is limited to simple geometries in current fabrication techniques [8][9].

For obtaining exceptional properties, many methods of producing porous metals have been developed around the world, such as sol-gel chemistry, the use of space-holders, and additive manufacturing [10]. The latter is a layer-by-layer approach allowing quick and cost-effective manufacturing of highly complex parts in small quantities without the necessity of specific tools. Using porous structures in manufacturing has a few advantages, including the cost reduction as compared to the powder bed fusion as compared to other technologies such as inkjet or thermoplastic extrusion in which non eco-friendly binders are used.

The main objective of the study is to obtain porous materials by 3D printing, using an extrusion additive manufacturing (EAM) method and eco-friendly paste, to convert pastes to porous structures with ideal heat dissipation properties. Initially, pastes with properties suitable for micro-extrusion are formulated. The pastes must be sufficiently fluid to be extruded, but sufficiently viscous to hold the shapes during the printing process (extrusion and drying) and post treatment (debinding and sintering). For this purpose, the nature of the binder together with the morphology of the fillers and their volume percentage are controlled. The formulation of the Cu pastes was based on the work by J. Bourret et al. [11] using eco-friendly additives. The binder and pastes were characterized macroscopically. The powders are characterized by X-ray diffraction, scanning electron microscopy, and laser granulometry. Secondly, to control porosity, regular parts were printed to study the impact of heat treatments (debinding, sintering) on the structure and microstructure of the parts, which will be characterized. Finally, the physical properties of parts with controlled porosity were studied by heat dissipation and Vickers hardness.

Experimental method

Preparation of the Cu paste

The Cu pastes were prepared using a spherical Cu powder (ECKA ® Granules (Germany, GmbH) with a medium particle size d50 close to 20 µm. Glycerol (Alfa Aesar A16205.AP, Thermo Fisher GmbH, Germany) was used as a plasticizer. The binder was obtained by acoustic mixing during 3 minutes (Labram, Resodyn, USA) of 5 wt.% of a water-soluble polysaccharide powder with distilled water as a solvent. This polyssacharide is eco-friendly and found in a wide range of applications thanks in part to its thickening, stabilizing, gelling and binding properties: food, cosmetic, medical and pharmaceutical industries [12]. Their nontoxicity, life cycles and low cost make them ideal natural additives. Due to its high hydrophilic properties, it is able to form a network due to the swelling of their particles [START_REF] Lapasin | Rheology of Industrial Polysaccharides: Theory and Applications[END_REF]. After preparation, pseudo-gels are obtained and were kept in the refrigerator at 8 °C for around 24 hours.

Cu powders were introduced into a container with the pseudo-gel at room temperature (binder + plasticizer).

The Cu pastes were prepared by mixing all these components in an acoustic mixer for 3 minutes. Once mixing was finished, the Cu paste was then kept in the refrigerator at 8 °C for a minimum of 24 hours of rest before being used in micro-extrusion. In our study, the formulation for the micro-extrusion contains 50 vol.% of Cu particles. This ratio corresponds to a best compromise between suitable rheological properties of the paste (with low viscosity) and the maximum ratio of metallic particles in the extrusion paste [11].

3D printing

The micro-extrusion tests have been performed with a 3D printer (S600D, Lynxter, Bayonne, France) with the PAS11 tool head. Simple 3D cubes (10x10x10 mm 3 ) parts were obtained from a CAD file; sliced using a software (Simplify 3D) to obtain the printing path with the computer G-code language executed by the printer.

The diameter of the nozzle and the distance between 2 strands are equal to 0.86 mm.

Consolidation by heat treatments

After extrusion, the parts were dried at room temperature and then debinded and sintered. Nabertherm RD 30/200/11 tube furnace is use for debinding and sintering. To determine the debinding cycle, thermal analysis was carried out on the binder powder. Several dwell temperatures of debinding (T° = 300, 400, 500, and 600 °C) in air and argon (Ar), with a flow of 200 mL/min, were tested to evaluate the impact of these parameters on the parts' density and microstructure. However, the heating ramp was fixed at 5 °C/min with a dwell time of 1 hour.

Sintering was used to obtain parts with a density lower than 60% to maximize specific surface areas with optimized thermal dissipation. Sintering was performed at 700 °C for 60 min in an Ar/H2(5%) atmosphere with a flow of 200 mL/min. The heating ramp is equal to 6 °C/min with a dwell time of 2 hours.

Structural and micro-structural characterization

The Archimedes method (Mettler Toledo AT201) was used to measure the real density of the samples and the mercury intrusion porosimeter to measure the pore volume and porosity of the Cu materials. The relative density was calculated by dividing the real density by the theoretical density of Cu, CuO and Cu2O (dCu = 8.96 g/cm 3 , dCuO = 6.31 g/cm 3 , dCu2O = 6 g/cm 3 ) knowing the amount of Cu, CuO and Cu2O inside the densified materials. The microstructure of printed parts was analyzed by Scanning Electron Microscopy (SEM, VEGA II SBH, TESCAN). The surface area is the main factor in thermal management. To demonstrate the benefit of additive manufacturing in thermal management, the heat dissipation performance of different structures printed was measured and compared to a dense Cu. In the experimental set-up [START_REF] Constantin | Laser 3D printing of complex Copper structures[END_REF], the heat sink was attached to a memory card chip placed on a hot plate heated to 200 ºC. A fan was located 10 cm from the top of the heat sink to produce an airflow. The temperature of the electronic chip was measured using an IR camera about 0.5 m from the heat sink with a precision of ± 2 °C.

Results and discussion

Debinding of the Cu materials printed

The diffractograms in Figure 1 show that during debinding under argon, whatever the debinding temperature, the structure of the parts is identical to that of the starting powder; it is highly crystallized cubic copper with well-defined peaks. On the other hand, when debinding is carried out in air, the structure evolves. CuO and Cu2O peaks appear from 200 °C. Figure 1: Typical diffraction diagram of Cu materials debinding under argon and ai5 atmosphere at different temperature The predominance of the phases has been plotted as a function of the debinding temperature in air in Figure 2. Figure 2a shows that when the debinding temperature increases, the Cu phase is gradually consumed in favor of the oxide phases; it is completely consumed at 600°C. The proportion of the Cu2O phase stabilizes between 20 and 30% while the CuO phase becomes predominant at 525°C. The Cu2O phase does not increase because it serves as a reaction intermediate for the formation of the CuO phase. In Fig. 2b, it notes the variation of the Cu part volume in mm 3 , in relation with the debinding temperature. As it can be seen in the Fig. 2a, the concentration of Cu oxide increases for a temperature greater then 300 °C and therefore shows the transformation of metal Cu to two different Cu oxides, Cu2O and CuO. Cu oxides have lattice parameters greater than that of pure Cu; this induces volume differences in the unit cell: VCu = 47.24 Å, VCu2O = 77.83 Å and VCuO = 81.22 Å [START_REF] Silvain | [END_REF]. These variations of the cell volume induced the increasing of the volume of the Cu parts during the debinding under air. Fig. 2b shows that for the parts debinded in argon, the volume is quite stable up to 400 °C and slightly increase after this temperature. For parts debinded in air, the volume of the parts is almost stable up to 300 °C, and then increases afterwards. This behavior results from a competition between two phenomena: A first phenomenon of densification caused by pre-sintering of Cu particles is preponderant up to 300 °C. For temperature higher than 300 °C, oxidation of Cu is predominant. This observation agrees with the micrography in Fig. 3a which shows that there is no visible evolution of the Cu powders after a heat treatment under Argon, which compare with the same temperature under air (cf. Fig. 3b and3c). Fig. 3c shows the growth of nanowires on top of the Cu particles. The CuO nanowires start appearing when the temperature reaches 400 °C and grow continuously as the annealing time increases. It is assumed that the formation of these CuO nanowires is driven by the thermal oxidation of Cu during annealing [START_REF] Silvain | [END_REF]. Such a phenomenon has already been observed by Jiang et al. [16]. A step-by-step kinetic study of the growth process shows that the CuO nanowires start appearing when the temperature reaches 400 °C and grow continuously as the annealing time increases, reaching lengths of several microns after 2 h of annealing. It is assumed that the formation of these CuO nanowires is driven by the thermal oxidation of copper during annealing. Cuprous oxide (Cu2O) is known to be the first product of Cu oxidation, while CuO is formed in a second step. Many mechanisms have been proposed in order to explain this phenomenon, such as the catalytic action of Cu2O on CuO formation through a vapor-solid mechanism [16], or the growth stresses in the oxide layers strong enough to extrude CuO wires from the surface [17]. The evolution of the microstructure, for material with a debinded in air, is correlated to the formation of Cu oxide. In contrary, the microstructures of the Cu parts with debinding in Ar, are quite similar at any debinding temperatures due to the absence of the Cu oxidation in Ar atmosphere.

Figure 3: SEM microstructure of the Cu materials debinded at 400 °C under a) Ar, and b, c) air (white arrows shows CuO nanowires) Fig. 4a shows that debinding conditions have a significant impact on the final sample density. For samples debinded in air, as the debinding temperature increases, the density of the parts increases from 200 to 400 °C, then decreases from 400 to 600 ° C. This behavior results from a competition between the densification caused by the sintering (decrease of the porosity inside the material) and the formation of Cu oxide (Cu2O and CuO), which have a density lower than pure Cu. When the samples are debinded in Ar, there is no Cu oxide formation; the density of the parts therefore increases linearly with the sintering temperature. Above 400 °C, the decrease of the measured density for Cu material debinded in air has to be correlated with the increase of the CuO and Cu2O percentage. Fig. 4: Evolution of the porosity level, inside the Cu material, after a) 60 min debinding in air and Ar atmosphere and b) debinding and sintering at 700 °C for 60 min in an Ar/H2(5%) atmosphere.

Sintering of the Cu materials printed

After debinding at different temperatures and atmospheres, Fig. 4b shows the evolution of the density of the sintered at 700 °C for 60 min in an Ar/H2(5%) atmosphere. First, the sintered densities are higher for materials debinded in air than those debinded under Ar at any temperature. Second, the variation of the density of the sintered materials shows the same trend that the variation of just debinded material. For samples debinded in air, an increase in the density is observed first for temperatures up to 400 °C, followed by a decrease at higher temperatures. For samples debinded in Ar, a continuous increase can be observed.

Heat dissipation of the materials 3D printed

Fig. 5 shows the evolution of the temperature difference between the hot plate and the electronic chip as a function of the porosity level of the sintered Cu parts. It has to be mentioned that, for this paper, just selective materials with micrometric porosity, induced by different debinding temperature and sintering conditions, ranging from 5 to 30% have been analyzed. All of these materials have been debinded under air. In a future paper, macro porosities, induced by EAM process, will be associate with micro porosities in order to increase the porosity range and to see the effect of the porosity size and shape with the thermal dissipation behavior. The heat dissipation in the porous Cu part significantly depends on the porosity level inside the Cu part sintered. Ended, when the volume fraction of pores increases, the temperature difference between the hot plate and the electronic chip increases and so does the ability of the porous Cu to dissipate heat. Therefore, it is possible to improve the heat dissipation of a material by increasing its specific surface areas through a higher porosity. It has to be mentioned that, due to the fact that the porosity level inside the sinter Cu materials is much higher than 5%, the structure porosity must be consider has open. This open porosity allows a better dissipation of the heat on the surface and inside the Cu parts. Figure 5: Evolution of the temperature difference between the hot plate and the electronic chip in function of the porosity level of the Cu materials.

Conclusion

Cu pastes with optimal properties for micro-extrusion have been formulated and porous Cu materials have been fabricated by 3D printing, using the EAM method. The influence of the debinding and sintering conditions (atmosphere, temperature, and time) on the porosity level and the heat dissipation behavior have been investigated. The results could be summarized as follows.

1) Debinding conditions have a significant impact on sample density. For samples debinded in air the density of the parts increases from 200 to 400 °C, then decreases from 400 to 600 ° C. This behavior results from a competition between the densification induced by the sintering and the Cu oxidation. When the samples are debinded in Ar, there is no oxide formation. Cu parts debinded in air are more porous than those in argon.

2) The variation of the density of the sintered materials show the same trend that the variation of just debinded material. For the samples debinded in air, an increase in the density is observed up to 400 °C, followed by a decrease at higher temperatures. 3) For Cu materials, the porosity level inside the materials 3D fabricated can be controlled by changing the debinding and sintering conditions. 4) Heat dissipation of the porous Cu part strongly depends on the porosity level inside the sintered Cu parts. For such open porosity structure, the increase of the volume fraction of pores is strictly correlated with the ability of the porous Cu to dissipate heat. 5) Finally, this innovative work shows clearly that the EAM method, using eco-friendly paste, can be a promising method to product the heat sink with a very heat dissipation ability thank to a specific internal porosity of material sintered and also the specific 3D shape of Cu heat sink.

Figure 2 :

 2 Evolution of a) the chemical composition, with air debinding temperature, of the Cu materials and, b) the volume of the materials, in mm 3 of the Cu materials at different air debinding temperatures

The authors would like to thanks Eric Lebraud of the Institut de Chimie de la Matière Condensée de Bordeaux for performing the XRD analyses.

This research did not receive any specific grant from funding from funding agencies in the public, commercial, or not-for-profit sectors.

All data generated or analyzed during this study are included in this published article

None