Prof M Pesavento

Astrid Rubiano

Prof P Larzabal

Dr R Boyer

Dr C Ren

Keywords: MIMO radar DOD Smart control of a soft robotic hand prosthesis Publications: [IC18, IC16, IC15] [J16, JS1, IC33, IC25, IC24, IC22, NC8, NC5, ICS1] [JS5, J14, IC27, IC23, IC19, ICS2, NC6] Adaptive scheme for urban and forest target detection using synthetic aperture radar Contribution to robust adaptive signal processing in MIMO systems without secondary data Performances investigation in the context MIMO radar performance analysis under K-distributed clutter Distributed multi-frequency calibration for sensor array radio interferometers Publications: [IC28, J18] Joint maximum likelihood calibration and DOA estimation with separated arrays Robust covariance matrix estimation in convex structures Publications: [IC30] Distributed processing for very large arrays in radioastron-

• Award of the allocation ministérielle MSER (from the French ministry of education and research) as Ph.D. student (2008Ph.D. student ( -2011)).

• Expert of the French National Research Agency (ANR 2017 and ANR 2018).

• Member of the board "Special Area Team -Signal Processing for Multisensor Systems (SAT-SPMuS)" EURASIP, since Jan. 2018.

• Member of the board "Special Area Team -Theoretical and Methodological Trends in Signal Processing (SAT-TMTSP)" EURASIP, since Aug. 2015.

• Aliated member of the IEEE Big Data SIG (Special Interest Group) since Jan. 2016.

• • Tutorial and plenary talks Tutorial at the international conference DAT-2017, Algeria, Fev. 2017, entitled:

"Lower bounds on the MSE : From Theory to Practice".

Plenary talk in the RADAR-EMP 2015 conference, at Military Polytechnic School, Algiers, May, 7th 2015, entitled: "MIMO Radar: Performance Analysis Investigation".

Invited presentation at the rst SKA-France meeting for signal processsing, Sept.,

8th, 2016 at ENS Paris-Saclay, France. SKA-France (Square Kilometre Array) is a national coordination of industrial, technical and scientic activities preparatory to the SKA project in France, set in place jointly by the Institute for Earth Sciences and Astronomy (CNRS/INSU), Paris Observatory, Côte d'Azur Observatory, Bordeaux University and Orléans University. (Slides are available at https://skafrance.oca.eu/fr/evenements/atelier-scientiques/14-ateliers-scientiques).

• Scientic committee/technical program committee • Co-organization of scientic thematic days GDR-ISIS Entitled: "Array processing: non Gaussian, non circular and non stationary signals",

Dec. 2016 at Telecom ParisTech, France (with Prof. J.-P. Delmas and Prof. P.

Chevalier).

• Elected member of the LEME laboratory board (Paris Nanterre University) since Apr. 2014.

• Head of the ES-group (Equipe Signal ) since Jan. 2018.

• Invited member of the UFR's board (UFR SITEC of Paris Nanterre University) since Sept.

2014.

• Correspondant GDR-ISIS (information, signal, image, vision) for the LEME-EA4416 laboratory of Paris Nanterre University since Avril 2015.

• In charge of the EESC Master of Science (Responsable M1-EESC ) Sep. 2018 at L2S-France).

• Member of 5 assistant professor/PRAG hiring committee since 2015 (3 for Paris Nanterre University (MCF-61, MCF-63 and PRAG-Maths), 2 for Centrale-Supelec (2 LRU equivalent to MCF-61)) .

• Responsible of the work-package integration of calibration techniques developed at ENS Paris-Saclay in the SKA-France (my team is composed of 2 assistant professors, 1 adjoint astronomer, 1 astronomer and 1 Ph.D. student).

• Reviewer for International Journals: ing by deriving lower bounds on the mean square error as well the statistical resolution limit.

Specically,

• The rst part of my Ph.D. had been dedicated to the calculus of the Cramér-Rao bound (CRB) for dierent data models under Gaussian noise assumption adapted for the asymptotic area (i.e., for a high signal to noise ratio and/or with innite number of observations.).

More precisely, I derived and analyzed the so-called conditional and unconditional CRBs for a single time-varying near-eld source. In each case, I obtained non-matrix closed-form expressions. This calculus has two advantages: i) due to the fact that one has to inverse the Fisher information matrix (FIM), the computational cost for a large number of snapshots (in the case of the conditional CRB) and/or for a large number of sensors (in the case of the unconditional CRB), of a matrix-based CRB can be high while the proposed approach is low and ii) some useful information had been deduced from the behavior of the bound.

In particular, an explicit relationship between the conditional and the unconditional CRBs was provided.

In the same vein, in order to characterize the optimal performances in the non-asymptotic region, I derived the McAulay-Seidman, the Hammersley-Chapman-Robbins, the McAulay-Hofstetter bounds and the so-called Todros-Tabrikian bound.

• In the second part of my Ph.D., I focused on the concept of the statistical resolution limit (i.e., the minimum distance between two closely spaced signals embedded in an additive noise that allows a correct resolvability/parameter estimation.) More precisely, I dened and derived the statistical resolution limit using the CRB and the hypothesis test approaches for the mono-dimensional case, i.e., one parameter of interest per source/signal.

Then, I extended this concept to the multidimensional case, i.e., multiple parameters of interest per source/signal. Applications that I considered in this context go from the polarized sources (sensed using the so-called COLD array; co-centered orthogonal loop and dipole array) to the study of the resolvability for source localization with known clutter interference in a MIMO radar context. First, I summarizes my research activity in the estimation context:

• Since January 2014, I am a member of the ANR MAGELLAN (Machine learning methods for the very large arrays in radio astronomy), and more specically, I am in charge of the task 2.2 related to the advanced calibration methods for the next radio interferometers generation. In this context, I am working, with my colleagues from Lagrange-University of Nice and SATIE-ENS Paris-Saclay, in proposing new robust and fast self-calibration algorithms. The radio interferometers context, imposes us to i) work with correlations a.k.a. visibilities (and not directly with the sensed measurements them self as it is usually the case in most array processing applications), ii) solve the unknown antenna gains and phases as well as the unknown atmospheric and ionospheric disturbances, iii) to handle a large number of elements and a large eld of view and nally, iv) robustify the calibration process since radio interferometer observations are often aected by the presence of outliers which are due to several causes, e.g., weak non-calibrator sources or man made radio frequency interferences. This makes calibration a daunting parameter estimation task for which the existing methods are ineective.

The aforementioned topics have been considered and/or are still currently been studied • An other important research track was initiated during my guest visits in Germany with a collaboration with Prof. M. Pesavento from Technische Universitat Darmstadt: designing novel estimation procedure in the context of parameterized mean and variance under compound Gaussian (CG) distributed clutter without secondary data (i.e., without assuming the existing of target-free signals). This rises, for example, in the recent high resolution radars for which the central limit theorem is not valid any more and thus the radar clutter cannot be correctly modeled as a Gaussian process. In this context, we recently devised in [JS3,J13,IC23,IC15] the conditional (i.e., the relaxed iterative maximum likelihood estimator), the joint (i.e., the iterative maximum a posteriori estimator) and the marginal (i.e., the exact iterative maximum likelihood estimator) maximum-likelihood-based iterative estimator for parameterized mean (for direction-of-departure and direction-of-arrival estimation in the MIMO radar context as an application). The proposed estimators employ a stepwise numerical concentration approach w.r.t. the objective function related to the conditional, joint and marginal likelihoods of the observation data. Our estimators leads to superior performance w.r.t. the existing state-of-the art. Interconnections, pros and cons of these three proposed estimators are discussed in [JS3, J13].

The aforementioned studies have been considered in the Ph.D. thesis of X. Zhang.

• In the same vein, a fructus discussions with my colleagues of SATIE laboratory (Prof. P.

Forster), SONDRA-CentraleSupelec (Dr. C. Ren) and LEME-Paris Nanterre University (Dr. A. Breloy) lead us to consider the problem of estimating covariance matrices in convex structure. More precisely, in [JS4, [START_REF] Mériaux | Ecient estimation of scatter matrix with convex structure under T-distribution[END_REF][START_REF] Mériaux | Robust-COMET for Covariance Estimation in Convex Structures: Algorithm and Statistical Properties[END_REF], we are considering the estimation of structured covariance matrices of a complex elliptically symmetric (CES) distributed observations (the complex elliptically symmetric distributions are a generalization of the CG distribution). We take into account the specic structure of the covariance matrix (e.g., Toeplitz structure appears in array processing with uniform linear array) and the non gaussianity of the data in order to improve substantially the estimation accuracy. Specically, we are tackling this problem by proposing a novel estimator, named StructurEd ScAtter Matrix Estimator (SESAME), which is based on the EXtended Invariance Principle (EXIP). In addition, we are conducting theoretical analysis on the unbiasedness and the asymptotic eciency of the proposed SESAME. Finally, we are proposing an iterative procedure of the proposed SESAME, called Iterative-SESAME (I-SESAME), reaching faster the CRB.

This study is currently addressed in the Ph.D. thesis of B. Mériaux.

• One of my current interest is the problem of low dimensional signal subspace estimation in a Bayesian context [START_REF] Ben Abdallah | Bayesian robust subspace estimation in presence of compound Gaussian sources[END_REF][START_REF] Ben Abdellah | Bayesian Low-Rank Signal Subspace Estimation for Compound Gaussian Sources[END_REF][START_REF] Ben Abdallah | Minimum Mean Square Distance Estimation of Subspaces in presence of Gaussian sources with application to STAP detection[END_REF]. This has been initiated with my colleague of Paris Nanterre University, Dr. A. Breloy and the with the help of Pr. D. Lautru, with whom, we obtained a scholarship (allocation ministérielle MSER) for R. Ben Abdellah as Ph.D.

candidate. The essence of this work is to devise new estimators using the so-called minimum mean square distance (MMSD), which minimizes the expected natural distance between the true range space of interest and its estimate taking into account the presence of a subspace basis priori distribution. Such approach has shown its interest for signal subspace estimation for small sample support and/or low signal to noise ratio contexts. Following this framework, we are deriving the MMSD of the signal subspace estimators in the context of CG sources embedded in white Gaussian noise, which is an useful model for various array processing applications. As a byproduct, we also introduce a generalized version of the complex Bingham Langevin distribution in order to model the prior on the subspace unitary basis.

As mentionned above, this work is under consideration in the Ph.D. thesis of R. Ben Abdellah.

• Miscellaneous: For sake of brevity some of my recent works are not detailed in this manuscript. As an example, the reader is referred to [J19,J12] for topics related to autoregressive based methods for the generalized likelihood ratio test. Specically, in [J12],

we focus on the adaptive detection of range-spread target in CG clutter without secondary data, whereas in [J19], we set up a sequential detection test in non-Gaussian correlated clutter using bootstrap framework.

This works have been done in collaboration with Prof. A. Zoubir from Technische Universitat Darmstadt, Germany.

Second, as mentioned above, I maintains parallel research activities related to new problems in the context of array processing performance analysis. Nevertheless, it is worth mentioning that these topics dier from my previous Ph.D. works that had been essentially concentrated around the calculus of the CRB and the statistical resolution limit under Gaussian noise for dierent array processing applications. Specically, my recent works related to the performance analysis context, can be categorized as follows

• Our recent paper [J17] is a result of a collaboration with Prof. A. Zoubir and his Ph.D. student, A. Mennad. In this latter, we derived Slepian-Bangs-type formulas for CES distributed data vectors in the presence of model misspecication. The basic Slepian-Bangs (SB) formula has been introduced in the array processing literature as a convenient and compact representation of the FIM for parameter estimation under (parametric) Gaussian data model. In [J17], we provided a new generalization of the classical SB formula to parametric estimation problems involving i) non-Gaussian, heavy-tailed and CES distributed data in ii) the presence of model misspecication. Moreover, we showed that our proposed formulas encompass the special cases of the SB formula for CES distributions under perfect model specication, the SB formulas in the presence of misspecied Gaussian models, and the SB formula for the estimation of the scatter matrix of a set of CES distributed data under misspecication of the density generator.

The aforementioned study has been considered in the Ph.D. thesis of A. Mennad.

• much attention as it plays a signicant role in several signal and array processing applications. Nevertheless, the study of the optimal estimation performances in such context is a dicult task since the unknown parameter vector of interest may contain both continuous (observations distribution's parameters) and discrete (change-point locations) parameters.

Our idea is to handle this by deriving a lower bound on the mean square error. More precisely, a Hybrid Cramér-Rao-Weiss-Weinstein bound (HCRWWB) and its associated closed-form expressions, whatever the considered distribution of the data, is proposed.

Furthermore, contrary to several works about performance analysis in the change-point literature, our study is adapted to multiple changes but also to unknown observations distribution's parameters [JS5,J14]. In addition, a semidenite programming formulation of the minimization procedure is given in order to compute the tightest proposed HCRWWB in an ecient way. This latter consists of nding the unique minimum volume covering the set constituted by hyper-ellipsoid elements which are generated using the derived candidate HCRWWB matrices w.r.t. the so-called Loewner partial ordering.

The aforementioned work is currently being studied in the Ph.D. thesis of L. Bacharch (funded by an allocation ministérielle MSER).

• Miscellaneous: Again, for sake of concisely some of my recent works are not detailed in this manuscript. The reader can refer to [J15,J8,IC19] for more details. As an example, and briey speaking, on the one hand, in [J15,J8] and in collaboration with Prof. E. Chaumette, we addressed, respectively i) the problem of fundamental limitations on resolution in deterministic parameters estimation by introducing a new rigourous denition of resolvability based on a probabilistic framework and incorporating a requirement for accuracy unlike most existing denitions; and ii) we proposed a new class of Weiss-Weinstein bounds (essentially free from regularity conditions on the probability density functions support) and we discussed its relationship with the Bobrovsky-Mayer-Wolf-Zakai bounds in [J15].

On the other hand, in [JS6,J11] with collaboration with Prof. J.-P. Delmas, we focused on the geometry design of planar antenna array in the context of near-eld source localization with and without the presence of a variable power prole. Specically, we studied the class of square and cross-based centro-symmetric arrays and highlight some of their attractive features. In particular, we identied key geometric parameters that control the near-eld array performance. Opportunistically, these geometric parameters were used in order to design non-uniform square and cross-based centro-symmetric arrays that achieve better near-eld localization accuracy. Such design is handled by minimizing the relative peak sidelobe level ratio derived from the conventional array beampattern. This latter, is a max-min problem under constraints, which can be transformed into a nite sequence of convex linear matrix inequality problems solved using the Matlab GloptiPoly utility.

After this brief introduction, in the following, I details some of my current activities. It is worth mentioning that the following research activities description has been partially presented and/or extracted from my publications which are listed in section.0.6. Consequently, some slight changes in notations from chapter to chapter may exist.

Notations: Chapter 2

Robust scatter matrix estimation and subspace estimation with application to radar

In most of the array processing literature, the additive noise or clutter (in the radar applications)

is simply assumed to be a Gaussian stochastic process. Such assumption is generally a good approximation in many cases and has its theoretical basis in the central limit theorem. However, in certain specic scenarios, the radar clutter cannot be correctly described by the Gaussian model anymore. As an example, experimental measurements reveal that the ground clutter data heavily deviate from the Gaussian model [START_REF] Billingsley | Ground clutter measurements for surface-sited radar[END_REF]. This is also true, e.g., for the sea clutter in a high-resolution and low-grazing-angle radar context, where the scatter number is random and the clutter shows nonstationarity [START_REF] Gini | Performance analysis of two adaptive radar detectors against non-Gaussian real sea clutter data[END_REF]. To account for such problems, where the noise and/or the clutter is a non-Gaussian process, numerous clutter models have been developed. Among them,

• The so-called complex elliptically symmetric (CES) distributions (cf. Section. 2.1.1.1).

Its main advantage lies in its generality to encompass a wide variety of non-Gaussian distributions (Generalized Gaussian, compound Gaussian, t-distribution, W -distribution and K-distribution, etc.) which turns out to model accurately impulsive noise, spiky radar clutter measurements or other heavy-tailed observations [START_REF] Sangston | Adaptive detection of radar targets in compound-Gaussian clutter[END_REF][START_REF] Ward | Maritime surveillance radar. part 1: Radar scattering from the ocean surface[END_REF].

• An other alternative is the compound-Gaussian (CG) distribution (referred as the spherically invariant random process (SIRP) in the radar community [START_REF] Ward | Maritime surveillance radar. part 1: Radar scattering from the ocean surface[END_REF]). The CG is a subclass of the CES and is a good alternative; while parametric estimator based on the CG modeling are less complex than those based on CES modeling, the CG own the feasibility to describe dierent scales of the clutter roughness, as well as to encompass some important heavy-tailed distribution as K-distribution, t-distribution, Laplace, Cauchy and Weibull distribution, etc (cf. Section. 2.3.1.2). The CG is a two-scale, complex process with random power, structured as the product of two independent components: a complex Gaussian process with zero mean and unknown covariance matrix, and the square root of a positive scalar random process. As an example, in the radar context, the former describes the local scattering and is usually referred to as speckle, while the latter, modeling the local power changing, is called texture. advantage of being density generator free, but at the cost of a loss of information and thus loss of optimal estimation accuracy. The CAE modeling is not discussed in this Chapter due to space limitation. Nevertheless, the reader is referred to our paper [START_REF] Meriaux | Robust-COMET for Covariance Estimation in Convex Structures: Algorithm and Statistical Properties[END_REF] in which we focused on structured covariance matrix estimation in a robust statistical framework using the CAE modeling.

In the following, I present some of our recent works related to the scatter, covariance and signal subspace robust estimation as well robust source localization without the use of secondary data. This Chapter is based on the following articles the Toeplitz structure appears in array processing with Uniform Linear Array (ULA) or in time series analysis [START_REF] Fuhrmann | On the existence of positive-denite maximum-likelihood estimates of structured covariance matrices[END_REF]. In MIMO communications or spatio-temporal noise processes in MEG/EEG data, the CM exhibits a Kronecker structure, where the factor matrices could be themselves structured [START_REF] Wirfält | On kronecker and linearly structured covariance matrix estimation[END_REF][START_REF] Werner | On estimation of covariance matrices with kronecker product structure[END_REF]. In some applications, the CM lies in a small-dimension subspace [START_REF] Forster | Generalized rectication of cross spectral matrices for arrays of arbitrary geometry[END_REF][START_REF] Combernoux | Generalized rectication in the l1-norm with application to robust array processing[END_REF].

In the afore-mentioned works, the CM estimation is unequivocally improved, when the prior structure is considered. However, they usually assume complex Gaussian distributed samples.

Then, the structured CM estimation is realized either by projecting onto a subset describing the structure of the Sample Covariance Matrix (SCM), which is the unstructured Maximum Likelihood (ML) estimator for the CM [START_REF] Greenewald | Robust sar stap via kronecker decomposition[END_REF][START_REF] Ottersten | Covariance matching estimation techniques for array signal processing applications[END_REF] or by deriving the ML with constraints. In some practical applications, performance is degraded, because the assumption of Gaussian distribution is inappropriate and the previous algorithms are not robust to outliers. In order to overcome this issue, a wide class of distribution free methods based on the unstructured Tyler's estimate has been proposed [START_REF] Meriaux | Robust-COMET for Covariance Estimation in Convex Structures: Algorithm and Statistical Properties[END_REF]1618]. Those methods begin by normalizing the zero mean observations to get rid of the texture. Specically, in [START_REF] Meriaux | Robust-COMET for Covariance Estimation in Convex Structures: Algorithm and Statistical Properties[END_REF], we proposed a robust extension of the COvariance Matching Estimation Technique (RCOMET). In [START_REF] Sun | Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions[END_REF][START_REF] Breloy | Robust rank constrained kronecker covariance matrix estimation[END_REF], estimators have been proposed which minimize a constrained version of Tyler's cost function using iterative Majorization-Minimization algorithms. In [START_REF] Soloveychik | Tyler's covariance matrix estimator in elliptical models with convex structure[END_REF], a COnvexly ConstrAined (COCA) CM estimator is presented, based on the generalized Method of Moments (MoM) for the Tyler's estimate subject to convex constraints.

However, normalizing the observation causes a loss of information, notably the scaling factor, which is the reason why we chose to work directly with a larger class of distributions, namely, the CES distribution [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF].

In this section, we introduce a StructurEd ScAtter Matrix Estimator (SESAME) for any given CES distribution, whose scatter matrix (which is proportional to the covariance matrix) owns a convex structure. This method is carried out in two steps. In the same vein as COMET but generalized to CES distributions, SESAME combines the unstructured ML estimation of the scatter matrix and the EXtended Invariance Principle (EXIP) [START_REF] Söderström | System Identication[END_REF]. A theoretical analysis of DISTRIBUTION 27 the proposed SESAME's asymptotic performance (consistency, bias, eciency and asymptotic distribution) is conducted. Also, we propose an iterative implementation of SESAME that leads to a second method, called Iterative-SESAME (I-SESAME), which performs better than SESAME in the sense that the asymptotic behavior is reached faster. It is worth noting that the potential extra-parameters characterizing the distribution are assumed to be known. In practice, they could be jointly estimated with the unstructured ML estimator of the scatter matrix [START_REF] Fortunati | Matched, mismatched, and robust scatter matrix estimation and hypothesis testing in complex t-distributed data[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized gaussian distributions[END_REF] during the rst step of the proposed algorithms.

Background and problem setup

In this section, we present briey the CES distribution model, the associated scatter matrix estimators and the corresponding compact form of the FIM. A fuller survey on CES distribution can be found in [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF].

CES distribution

A m-dimensional random vector (r.v.), y ∈ C m follows a CES distribution if and only if it admits the following stochastic representation [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]:

y d = m y + QAu (2.1)
where m y ∈ C m is the location parameter and the non-negative real random variable Q, called the 2nd-order modular variate, is independent of the complex r.v. u. The latter r.v. is uniformly distributed on the unit k-sphere

CS k z ∈ C k | z = 1 , denoted by u ∼ U CS k . The matrix A ∈ C m×k has rank(A) = k. Furthermore, if it
exists, the mean of y is equal to m y and the covariance matrix of y is proportional to the scatter matrix, R = AA H , more precisely

E (y -m y )(y -m y ) H = 1 k E [Q] R.
In the following, we assume m y is known, or equivalently assumed equal to zero vector, without loss of generality and k = m, so rank(R) = m to belong to the absolutely continuous case. In this case, the probability density function (pdf ) of such a vector exists and can be written as [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]:

p Y (y; R, g) = C m,g |R| -1 g y H R -1 y (2.2)
in which the function g : R + → R + is called the density generator and satises

δ m,g +∞ 0 t m-1 g(t)dt < ∞ and C m,g = Γ(m) π m δ -1
m,g is the normalizing constant. In this case, we denote y ∼ CES (0, R, g) in short. Furthermore, thanks to (2.1), the quadratic form y H R -1 y d = Q has the following pdf:

p Q (q) = δ -1 m,g q m-1 g (q) . (2.3)
As already indicated, the CES distributions embrace various non-Gaussian distributions, like Generalized Gaussian, t-distribution, K-distribution. The expression of the density generator function is given in the Table 1 for some commonly used CES distributions.

M-estimators

Let us consider N i.i.d. zero mean CES distributed observations, y n ∼ CES m (0, R, g) , n = 1, . . . , N , with N > m.

M = 1 N N n=1 u y H n M -1 y n y n y H n , and R ∝ M (2.4)
In the case of the ML estimator, the function u(•) is given by u ML (s) = - g (s) g(s)

, where g (•) refers to the derivative of g(•). In the case of unknown density generator function g(•), the previous function u ML (•) is replaced by another, which is described later. Existence and uniqueness of the solution of (2.4) have been rst studied in the real case by [START_REF] Maronna | Robust M-estimators of multivariate location and scatter[END_REF][START_REF] Kent | Redescending M-estimates of multivariate location and scatter[END_REF], where the function u(•)

satises a set of general assumptions. The extension for the complex case has been proposed

in [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]. Furthermore, we consider the following equation, which is the limit towards which the equation (2.4) converges almost surely by the law of large number

M = E u y H M -1 y yy H (2.5)
where y ∼ CES m (0, R, g). Under the previous assumptions for the function u(•), the equation (2.5) (respectively (2.4)) admits a unique solution M (respectively M) and

M = σ -1 R
where σ is the solution of E ψ(σ|t| 2 ) = m with t ∼ CES m (0, I, g) and ψ(s) = su(s). In addition, the estimator M can be easily obtained by an iterative procedure and is consistent. Finally, the asymptotic distribution of M has been established for the complex case in [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF][START_REF] Mahot | Asymptotic properties of robust complex covariance matrix estimates[END_REF]:

√ N vec M -M d → CN (0, Σ, Ω) (2.6)
with

   Σ = σ 1 M T ⊗ M + σ 2 vec (M) vec (M) H Ω = σ 1 M T ⊗ M K + σ 2 vec (M) vec (M) T (2.7) and            σ 1 = a 1 (m + 1) 2 (a 2 + m) 2 σ 2 = 1 a 2 2 (a 1 -1) -a 1 (a 2 -1) m + (m + 2)a 2 (a 2 + m) 2 (2.8)
where K is the commutation matrix, which satises Kvec (A) = vec A T [START_REF] Magnus | The commutation matrix: some properties and applications[END_REF]. The matrix Σ (respectively Ω) denotes the covariance matrix (respectively pseudo-covariance matrix). The coecients a 1 and a 2 are dened by:

     a 1 = 1 m(m + 1) E ψ 2 σ|t| 2 a 2 = 1 m E σ|t| 2 ψ σ|t| 2 (2.9)
We remark that |t| 2 = t H t d = Q, thus in the following, we note A = E ψ 2 σ|t| 2 = E ψ 2 (σQ)

and B = E σ|t| 2 ψ σ|t| 2 = E σQψ (σQ) . The subscript ML is used when A and B are com- puted with ψ ML (s) = su ML (s).

EFFICIENT ESTIMATION OF COVARIANCE/SCATTER MATRICES WITH CONVEX STRUCTURE UNDER CES DISTRIBUTION 29

Regarding the ML estimator, we obtain σ = 1, since we prove that E [ψ ML (Q)] = m. Simi- larly, we can show that B ML = A ML -m 2 , thus the coecients σ 1 and σ 2 , fully describing the asymptotic (pseudo)-covariance matrices of the scatter ML estimator, can be simplied by

σ 1 = m(m + 1) A ML and σ 2 = -σ 1 (1 -σ 1 ) 1 + m(1 -σ 1
) .

(2.10)

Fisher information matrix

Applying [START_REF] Besson | On the Fisher information matrix for multivariate elliptically contoured distributions[END_REF] to the particular case of a zero mean CES distribution and a scatter matrix R parameterized by the real vector, µ, the (k, ) element of the FIM for a single vector of observation is given by:

[F (µ)] (k, ) = κ 1 Tr R -1 Rk R -1 R + κ 2 Tr R -1 Rk Tr R -1 R (2.11) in which Rk = ∂R ∂µ k , κ 1 = A ML m(m + 1)
and κ 2 = κ 1 -1. We recall that A ML = E ψ 2 ML (Q) . We rewrite (2.11) as

F (µ) = ∂r(µ) ∂µ H Y ∂r(µ) ∂µ (2.12) where Y = κ 1 R -T ⊗ R -1 + κ 2 vec R -1 vec R -1 H
and ∂r(µ) ∂µ refers to the Jacobian matrix of r(µ) = vec (R(µ)).

Problem setup

Let us consider N i.i.d. zero mean CES distributed observations, y n ∼ CES m (0, R e , g) , n = 1, . . . , N . We assume that the scatter matrix belongs to a convex subset S of Hermitian matrices, for which there exists a one-to-one dierentiable mapping µ → R(µ) from R P to S. The unknown parameter of interest is the vector µ with exact value µ e , and R e = R(µ e ) corresponds to the true scatter matrix. The log-likelihood function is given, up to an additive constant, by

L (y 1 , . . . , y N ; µ) = -N log |R(µ)| + N n=1 log g y H n R(µ) -1 y n (2.13)
The above function is generally non-convex w.r.t. R, its minimization w.r.t. µ is therefore a laborious and computationally prohibitive problem. To overcome this issue, we propose in the next section a new estimation method that gives unique, consistent and asymptotically ecient estimates. Furthermore, for linear structures, we obtain closed form expressions of these estimates.

SESAME : StructurEd ScAtter Matrix Estimator

In this section, we propose a two-step estimation procedure of µ. The rst step consists in computing the unstructured ML estimator of R e . The estimation of µ is then obtained by solving a weighted least squares problem derived from the so-called EXIP approach [START_REF] Söderström | System Identication[END_REF]. For notational convenience, we omit the dependence on N for the estimators based on N observations when there is no ambiguity.

R FP = 1 N N n=1 u ML y H n R -1 FP y n y n y H n H N ( R FP ) (2.14)
As already mentioned above, the iterative algorithm R k+1 = H N (R k ) converges to R FP for any initialization point [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]. Moreover, the consistency and the asymptotic Gaussianity, verifying (2.6), of this estimator are established in [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF][START_REF] Mahot | Asymptotic properties of robust complex covariance matrix estimates[END_REF].

Step 2

For the second step, guided in part by the EXIP principle [START_REF] Söderström | System Identication[END_REF], we estimate µ by minimizing the following criterion J R FP , R (µ):

µ = arg min µ J R FP , R (µ) with J R FP , R (µ) = κ 1 Tr R -1 R FP -R(µ) R -1 R FP -R(µ) + κ 2 Tr R -1 R FP -R(µ) 2 (2.15)
where R refers to any consistent estimator of R e , such as for instance R FP . In the same way to get (2.12), we obtain the following reformulation for (2.15) as

µ = arg min µ ( r FP -r(µ)) H Y ( r FP -r(µ)) = arg min µ Y 1/2 ( r FP -r(µ)) 2 (2.16) with Y = κ 1 W -1 + κ 2 vec R -1 vec R -1 H , W = R T ⊗ R and r FP = vec R FP .
In practice, we choose R = R FP and minimize the cost function J R FP , R FP (µ). Another consistent estimator of R will be used in Section 2.1.4.

Given R FP and R, the function J R FP , R (µ) is convex w.r.t R(µ). Therefore, for R ∈ S convex set, the minimization of (2.16) w.r.t. R(µ) is a convex problem that admits a unique solution. Consequently, the one-to-one mapping ensures the uniqueness of µ.

For a discussion on the practical implementation for holding the PSD constraint, the reader is referred to Appendix 6.1.1.

Application of SESAME to some CES distributions

As already mentioned, the CES distributions encompass a large number of standard distributions like the multivariate complex Generalized Gaussian, multivariate complex t-, W -and K-distributions. In this subsection, we provide explicit expressions for the constant A ML , ap- pearing in the SESAME criterion for these common CES distributions. The coecients σ 1 and σ 2 related to the scatter matrix unstructured ML estimator are also given. All the results are recapped in Table 1 for a centered m-dimensional complex random vector. The guidelines for the above calculations are given in Appendix 6.1.2.

In the following, µ will be referred to the SESAME estimate of µ, for which we address the study of consistency and eciency in section 2.1.3.

EFFICIENT ESTIMATION OF COVARIANCE/SCATTER MATRICES WITH CONVEX STRUCTURE UNDER CES DISTRIBUTION 31

1

Gaussian Generalized Gaussian Student W -distribution K-distribution 1 CNm CGN m,s,b Ct m,d CW m,s,b CKm,ν g(t) exp (-t) exp (-t s /b) s, b > 0 (1 + t/d) -(d+m) d > 0 t s-1 exp (-t s /b) s, b > 0 √ t ν-m K ν-m 2 √ νt ν > 0 Cm,g π -m sΓ(m)b -m/s π m Γ(m/s) Γ(m + d) π m d m Γ(d) sΓ(m)b -(m+s-1)/s π m Γ((m + s -1)/s) 2 ν (ν+m)/2 π m Γ(ν) A ML m(m + 1) m(m + s) m(m + 1)(m + d) d + m + 1 s(m + s -1) + m 2 2 -(m+ν) Γ(ν)Γ(m) R + x m+ν+1 K 2 ν-m-1 (x) K ν-m (x) dx σ 1 1 m + 1 m + s d + m + 1 d + m m(m + 1) s(m + s -1) + m 2
No closed form, numerical evaluation from (2.10)

σ 2 0 1 -s s(m + s) d + m + 1 d(d + m) m(1 -s)(m + s) s(m + s -1)(s(m + s -1) + m 2 )
Table 1: Density generator and their corresponding coecients for some commonly used CES distributions.

Asymptotic analysis

This section provides a statistical analysis of the proposed estimator SESAME, µ, which is the unique solution minimizing the criterion (2. is a consistent estimator of R(µ e ).

Proof. Using the consistency of R FP and R [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF] and for large N, we obtain Y P → Y e , r FP P → r e .

Consequently in (2.16), we obtain µ

P → µ ∞ where µ ∞ is the solution of the asymptotic criterion (2.16) for N → ∞ µ ∞ = arg min µ Y 1/2 e (r e -r (µ)) 2 .
(2.17)

Since Y 1/2
e is a full-rank matrix and the mapping is one-to-one, the unique solution of the above problem is µ ∞ = µ e with probability one, which establishes the consistency of µ. Finally, the continuous mapping implies R( µ) P → R(µ e ). Theorem 2 Let µ N be the SESAME estimate of µ e from N i.i.d. observations, y n ∼ CES m (0, R(µ e ), g). Therefore µ N is asymptotically unbiased, ecient and Gaussian distributed. Specically, Proof. The above expression of the CRB follows straightforwardly from (2.12). The estimate µ N is given by minimizing the function J R FP , R (µ). The consistency of µ N allows us to write the following Taylor expansion around µ e :

√ N ( µ N -µ e ) d → N (0, CRB) (2.
0 = ∂J R FP , R (µ) ∂µ µ= µ N = ∂J R FP , R (µ) ∂µ µ=µ e +   ∂ 2 J R FP , R (µ) ∂µ∂µ T µ=ξ N   ( µ N -µ e )
with ξ N on the line segment connecting µ e and µ N , i.e., ∃ c ∈ ]0, 1

[ such that ξ N = cµ e +(1-c) µ N [28, Theorem 5.4.8], leading to √ N ( µ N -µ e ) = -   ∂ 2 J R FP , R (µ) ∂µ∂µ T µ=ξ N   -1 √ N g N (µ
(µ) = ∂J R FP , R (µ) ∂µ .
• The consistency of µ N implies ξ N P → µ e . Moreover, by consistency of R FP and R, the continuous mapping yields to • Furthermore, the gradient g N (µ) is expressed by

∂ 2 J R FP , R (µ) ∂µ∂µ T µ=ξ N P → ∂ 2 J Re,Re (µ)
g N (µ) = -∂r(µ) ∂µ H Y( rFP -r(µ)) -( rFP -r(µ)) H Y ∂r(µ) ∂µ (2.20) = -2 ∂r(µ) ∂µ H Y( rFP -r(µ))
The gradient is so real-valued. Thus, using the asymptotic distribution of r FP given by (2.6), we can shown that (see Appendix 6.1.3 for details)

- √ N g N (µ e ) d → N 0, 4CRB -1 (2.
√ N ( µ N -µ e ) d
→ N (0, CRB), which concludes the proof on the unbiasedness, the asymptotic eciency and Gaussianity of SESAME estimator.

Iterative SESAME

In this section, we propose an iterative procedure for SESAME, which possesses the same asymptotic performance than SESAME and turns out to be a numerical improvement in most cases.

The SESAME algorithm provides an estimate µ of µ e by minimizing J R FP , R (µ), where J R FP , R (µ)

µ SESAME = arg min µ J R FP , R (µ)
with R FP is the unstructured ML estimator of R e obtained with (2.14) and R is any consistent estimator of R e . According to Theorem 1., R ( µ SESAME ) is a consistent estimator of R e .

Intuitively, the better the estimator R is, the better the solution µ should be. This points naturally lead to an iterative procedure, where the minimized norm is rened at each step by updating R with R ( µ SESAME ). For a nite number of steps, we obtain the Iterative SESAME (I-SESAME) for µ, denoted µ I-SESAME and achieved at the k-th stage by solving

µ (k+1) = arg min µ J R FP , R (k) (µ) with R (k) = R µ (k) (2.22)
The I-SESAME algorithm is recapped as By applying Theorems 1. and 2. at each iteration, the below theorem follows immediately Step 2 Initialize µ (1) by minimizing J R FP , R FP (µ) Step

3 for k = 1 to N it , Compute µ (k+1) from (2.22) Return µ I-SESAME = µ (N it +1)
Theorem 3 Let µ I-SESAME N be the I-SESAME estimate of µ e based on N i.i.d. observations,

y n ∼ CES m (0, R(µ e ), g). µ I-SESAME N is a consistent estimator of µ e . Likewise, R µ I-SESAME N is a consistent estimator of R(µ e ).
Moreover, it is asymptotically unbiased, ecient and Gaussian distributed:

√ N µ I-SESAME N -µ e d
→ N (0, CRB) .

(2.23)

The improvement of I-SESAME compared to SESAME lies in the fact that, as noticed by numerical simulations, I-SESAME reaches empirically faster the CRB than SESAME, i.e., with a smaller number of samples.

Applications and numerical results

This section presents the SESAME in the particular case of linear structure for the scatter matrix.

SESAME with linear parameterization

In the case of linear parameterization for the scatter matrix, which is indeed a convex structure, there exists a matrix, P ∈ C m 2 ×2m-1 , which relates the vectorized matrix R(µ) to µ according to:

r(µ) = vec (R(µ)) = Pµ (2.24)
In this case, the PSD constraint is not taking into account. Considering linear structure, the SESAME criterion (2.16) reads

µ = arg min µ Y 1/2 r FP -Y 1/2 Pµ 2 (2.25)
The well known analytical solution gives

µ = P H YP -1 P H Y r FP (2.26)
In this particular case, the CRB on µ for a single observation is easily shown to be We compare the performance of the proposed algorithms to the state-of-the-art and the CRB.

CRB = P H Y e P -1 (
(µ) =       R 1 R 2 • • • R m R * 2 . . . . . . . . . . . . . . . . . . R 2 R * m • • • R * 2 R 1       and µ =          R 1 (R 2 ) (R 2 ) . . . (R m ) (R m )          ∈ R 2m-1 (2.
Furthermore, we display the performance of SESAME estimation scheme by replacing the rst step by the joint-algorithm proposed in [START_REF] Fortunati | Matched, mismatched, and robust scatter matrix estimation and hypothesis testing in complex t-distributed data[END_REF] to deal with the possibility of unknown parameter d.

For the performance of I-SESAME, we consider only the unstructured ML estimator as rst step at each iteration. Our algorithms are compared to RCOMET from [START_REF] Meriaux | Robust-COMET for Covariance Estimation in Convex Structures: Algorithm and Statistical Properties[END_REF], COCA from [START_REF] Soloveychik | Tyler's covariance matrix estimator in elliptical models with convex structure[END_REF] and Constrained Tyler from [START_REF] Sun | Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions[END_REF]. The three methods are based on the Tyler's scatter estimator [START_REF] Tyler | A distribution-free M-estimator of multivariate scatter[END_REF] using normalized observations z n = y n / y n . It should be noted that, for Constrained Tyler, the Algorithm 3 in [START_REF] Sun | Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions[END_REF] derived for real-valued PSD Toeplitz matrices can not be applied. However, the Vandermonde factorization of PSD Toeplitz matrices [START_REF] Bäckström | Vandermonde factorization of toeplitz matrices and applications in ltering and warping[END_REF] allows us to use the Algorithm 2 of [START_REF] Sun | Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions[END_REF]. Thus, the structure set for Constrained Tyler becomes S = R | R = APA H with the diagonal matrix P 0 and A = [a(-90

• ), a(-88 • ), • • • , a(86 • ), a(88 • )],
where a (θ) = 1, e -jπ sin(θ) , • • • , e -jπ(m-1) sin(θ) T

.

Finally, we compare to the intuitive estimate µ obtained by averaging the real and imaginary parts of diagonals of the unstructured ML estimator, which corresponds to the Euclidean projection onto the Toeplitz set.

The asymptotic eciency of our estimator is checked on Fig. 2.1: its Mean Square Error (MSE) reaches the CRB as N increases. RCOMET, Constrained Tyler and COCA do not reach this bound since they do not take into account the underlying distribution of the data. Despite the absence of convergence proof for the joint-algorithm, we notice that optimal asymptotic performances for µ may be approached with unknown d. Performance of the proposed estimation scheme with the joint-algorithm as rst step are not displayed for small N , since the jointalgorithm does not converge for part of the 5000 runs. In addition, the asymptotic unbiasedness of SESAME as well as those of the other algorithms can be indirectly observed on the Fig. shows its limits outside these standards regimes. This estimator is also known to be sensitive to missmodeling, e.g., outliers or non-Gaussian observations. A possible solution to ensure better performance in these contexts is to incorporate a prior knowledge into the estimation process. In a Bayesian context, a prior distribution of the subspace unitary basis can be assumed in order to overcome the afromentionned drawbacks. This approach yield estimators such as the maximum a posterior (MAP) which maximizes the posterior probability [START_REF] Besson | Joint bayesian estimation of close subspaces from noisy measurements[END_REF], the minimum mean square error (MMSE), which minimizes the average Euclidean distance between the true unitary basis of subspace of interest and its estimate, and the minimum mean square distance (MMSD), that minimizes the distance between the true projection matrix and its estimate [START_REF] Srivastava | A bayesian approach to geometric subspace estimation[END_REF][START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF]. The latter is an intuitively appealing method, as it is based on a natural metric in the complex Grassmann manifold [START_REF] Chikuse | Statistics on Special Manifold[END_REF], i.e., the set of P -dimensional subspaces in C N (where P is the rank of subspace and N is the dimension of the observation space). In the context of subspace estimation, the MMSD has been addressed in [START_REF] Srivastava | A bayesian approach to geometric subspace estimation[END_REF][START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF]. More specically, [START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF] derives a practical formulation of the MMSD estimators when the subspace of interest is parameterized by its unitary basis. This formulation is then used in [START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF] to propose MMSD estimators involving an uniform prior for the sources distribution.

In this section, we focus on the context of sources following a CG distribution [START_REF] Ollila | Complex elliptically symmetric distributions survey, new results and applications[END_REF] embedded in white Gaussian noise. Indeed, the CG distributions have been considered in many modern array processing applications, as they can account for local power uctuations and presents good agreement to several real data set [START_REF] Ollila | Compound-gaussian clutter modeling with an inverse gaussian texture distribution[END_REF][START_REF] Greco | Statistical analysis of measured polarimetric clutter data at dierent range resolutions[END_REF]. Note that this familly covers a large panel of well known distributions, notably heavy-tailed distributions such as Weibull, Student's t-distribution, and K-distribution (see [START_REF] Ollila | Compound-gaussian clutter modeling with an inverse gaussian texture distribution[END_REF] and references therein). Hence, the considered mixture distribution can accurately model clutter (or power-uctuating sources) plus thermal noise observations, which are common in array processing. For this reason, we present herein, a new MMSD estimators in the context of CG distributed sources embedded in white Gaussian noise. The derivations are built upon the framework of [START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF] (which considered an uniform prior on the sources), and two Gibbs-sampler based algorithms are established to compute the proposed estimators. In addition, our development requires to extend the framework of [START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF] to the case of data with complex entries. To this aim, we consider a generalization of complex Bingham-Langevin (also U ∈ U N P the unknown unitary (orthonormal) basis of subspace of interest, S ∈ C P ×K the signal of interest and N ∈ C N ×K the additive noise, in which, N denotes the size of data, K the number of samples, P the dimension of subspace of interest (P N ) and p(Y|U) the conditional probability of Y given U. In the following subsection, we present the main background theory on which our derivations are based.

referred

Minimum mean square distance estimator

The MMSD estimator minimizes the average Euclidean distance between the true range space R(U) = UU H and its estimate R( U) = U U H , which corresponds to the natural distance between the subspace spanned by U and U in the complex Grassmann space. Straightforward, extending the formulation of [START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF] to the complex case, the MMSD estimator is expressed as:

U MMSD = arg min U E U U H -UU H 2 F = arg max U E{Tr{ U H UU H U}} (2.29)
Adopting the proof in [START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF] for the complex case, the MMSD estimator is obtained as

U MMSD = P P {M(p(U|Y))} (2.30)
with the matrix

M(p(U|Y)) = UU H p(U|Y)dU (2.31)
and where the operator P P {.} extracts the P strongest eigenvectors of a given matrix. This generic expression of the MMSD depends on p(U|Y). The latter has to be specied based on both the data model, the parameters S and U, and the prior distribution assigned onto the parameters. In the following sections, we detail the distributions on the source signal S and the unitary basis U that we consider. Specically, in order to be robust to various signal distributions, the CG with unknown deterministic textures [START_REF] Ollila | Compound-gaussian clutter modeling with an inverse gaussian texture distribution[END_REF] is chosen as a prior distribution for the sources.

Meanwhile, for the unitary basis U, we propose the use of the Complex Generalized Bingham Langevin (CGBL) distribution, as it extends several directional distributions [4042] to the case of complex data.

Compound Gaussian distribution

Let us assume that s follows a zero mean CG distribution, denoted s ∼ CG(0, Σ, f τ ), then, it has the following stochastic representation

s d = √ τ d , (2.32)
where i) τ is an positive random scalar, called texture, of probability density function f τ . This parameter is statistically independent of d. Depending on f τ , we can obtain various standard multivariate distribution for s (see [START_REF] Ollila | Complex elliptically symmetric distributions survey, new results and applications[END_REF] for more details). In order to propose algorithms that are robust to these distributions, we consider here this parameter as unknown deterministic for each realization. ii) d follows a zero-mean multivariate complex Gaussian distribution of covariance matrix Σ, denoted, d ∼ CN (0, Σ). The parameter Σ is referred to as the scatter matrix. Notice that if E{τ } < ∞, the covariance matrix of s exists and is proportional to the scatter matrix. ) is a diagonal matrix whose elements are the i.i.d. the texture realizations, aggregated in a vector τ .

Complex generalized Bingham Langevin distribution

The CGBL is adopted as prior distribution to gather the pre-established knowledge on the random basis U∈ U N P . The proposed CGBL is a probability distribution on the set of unitary matrices which combines linear and quadratic terms that is parametrized by a set of matrices {A p } ∈ S + R and the matrix C. We denote U ∼ CGBL(C, {A p }) when the pdf of U reads:

p CGBL (U) ∝ exp    P p=1 Re{c H p u p } + u H p A p u p    (2.33)
where c p and u p stands for the p-th column vector of respectively C and U.

{y k } k∈[[1,K]
] are assumed to be drawn as:

y k = s k + n k (2.34)
where

• s k ∼ CG(0, Σ, τ k ) are the CG distributed sources. τ ∈ R K + denotes the vector of textures τ k that are assumed to be unknown deterministic. Moreover, the source scatter matrix will be parameterized by its low-rank SVD as Σ = UΛU H . The eigenvectors U span the signal subspace basis, and we set a CGBL prior (of pdf p CGBL ) for this unitary matrix. The eigenvalues matrix Λ = diag(λ), where λ = [λ 1 , . . . , λ P ] T , is positive diagonal with deterministic unknown entries.

• n k ∼ CN (0, σ 2 I N ) is an additive white Gaussian noise of known or pre-estimated variance σ 2 . APPLICATION TO RADAR Hence, we have for each sample the representation

(y k |U, Λ, τ k ) ∼ CN (0, Σ k ) with Σ k = τ k UΛU H + σ 2 I N ∀k ∈ [[1 . . . K]] .
(2.35)

Thanks to the Sherman Morrison Woodbury lemma, the expression of Σ -1 k is simplied as

Σ -1 k = σ -2 I -UΓ k U H , where Γ k = σ -2 I P -(τ k Λ + σ 2 I P ) -1 is a diagonal matrix of entries [Γ k ] p,p = γ p k = τ k λ p σ 2 (τ k λ p + σ 2 )
.

(2.36)

We nally have the expression the probability density function of the sample set Y as

p(Y|U, λ, τ ) = K k=1 p(y k |U, λ, τ k ) ∝ K k=1 exp -y H k Σ -1 k y k det(Σ k )
.

(2.37)

Now that the model is established, the MMSD estimator of U is expressed as the solution of the following optimization problem: min

U,τ ,λ E U U H -UU H 2 F s.t. τ k ≥ 0 ∀k, λ p ≥ 0 ∀p Σ k = τ k U diag(λ) U H + σ 2 I N U H U = I P . (2.38)
In order to compute this estimator eciently, we derive in the following section an iterative algorithm that sequentially updates the variables U, τ and λ. The update of U requires a Gibbs sampling scheme, while for updating both the texture τ and the eigenvalues λ, we use a Majorization-Minimization (MM) procedure [START_REF] Sun | Majorization-minimization algorithms in signal processing communications and machine learning[END_REF], that allows to obtain closed form updates.

The overall algorithm is summed up in the box Algorithm 1, and details are given below.

Algorithm derivation

In the following, the superscript t denotes the value of the variable at iteration t of the algorithm. First, an initialization of the variables U 0 ,λ 0 , and τ 0 is set. This initialization can, for example, be taken from the P strongest eigenvectors and eigenvalues of the SCM for U and λ, and the norm of samples for τ . Now, we derive an algorithm that cyclically updates the blocks U, λ and τ by improving the value of the objective function.

Update the basis U

For xed blocks τ t and λ t . The update U t+1 is obtained by solving the problem:

max

U E U U H -UU H 2 F s.t. Σ t k = τ t k U diag(λ t ) U H + σ 2 I N U H U = I P (2.39)
Which corresponds to a MMSD estimation process. Thanks to the expression given in section 2.2.1.1, the update is obtained by

U t+1 = P P M(p(U|Y, τ t , λ t )) (2.
p(U|Y, τ t , λ t ) ∝ p(Y|U, τ t , λ t )p CGBL (U) ∝ K k=1 exp -y H k Σ t k -1 y k det(Σ t k ) p CGBL (U) ∝ K k=1 exp -y H k (-UΓ t k U H + σ -2 I N )y k p CGBL (U) ∝ exp K k=1 y H k UΓ t k U H y k p CGBL (U) ∝ exp    P p=1 u H p M t p u p    p CGBL (U) ∝ exp    P p=1 Re{c H p u p } + u H p A p + M t p u p    (2.42)
with matrix M p in (2.42) dened as

M t p = K k=1 γ t k,p y k y H k (2.43)
The posterior probability is therefore recognized as (U|Y,

τ t , λ t ) ∼ CGBL(C, {G t p }) with G t p = A p + M t p , ∀p ∈ [[1, P ]].
With this general distribution, there is no closed form for computing the expression of M(p(U|Y)) in (2.41). However, the update can still be evaluated using the so-called induced arithmetic mean (IAM) [START_REF] Sarlette | Consensus optimization on manifolds[END_REF] of the unitary matrix, as

U t+1 = P P    1 N r N bi +Nr n=N bi +1 U t (n) U t (n) H    (2.44)
where U t (n) are sampled as U t (n) ∼ CGBL(C, {G t p }), N bi stands for the burn-in samples (number of thrown samples from the Markov chain), and N r is the number of samples used to evaluate the update. In order to do so, an ecient Gibbs sampling procedure to draw the CGBL distribution as detailed in Appendix 6.2.

Update the eigenvalues λ and textures τ

For xed U t+1 the update of λ t+1 and τ t+1 boils down a maximum likelihood estimation process, since this parameter is considered unknown deterministic. Hence, updating λ requires solving the following problem max λ,τ p(Y| U t+1 , τ , λ)

s.t. τ k ≥ 0 ∀k, λ p ≥ 0 ∀p Σ t+1 k = τ k U t+1 diag(λ) U t+1 H + σ 2 I N (2.45)
The above maximization problem is equivalent to the minimization of the negative log-likelihood which is obtained via (2.37) as:

min λ,τ K k=1 ln det(Σ t+1 k ) + y H k Σ t+1 k -1 y k s.t. τ k ≥ 0 ∀k, λ p ≥ 0 ∀p Σ k = τ k U t+1 diag(λ) U t+1 H + σ 2 I N (2.46) APPLICATION TO RADAR input : Y, C, {A p }, σ 2 , P , K, N , N bi , N r output : U, τ , λ initialize: U (0) ← U init , τ (0) ← τ init , λ (0)
←

end

This problem has no closed form solution but can be solved by following the MM approach of [START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF].

The MM algorithm [START_REF] Sun | Majorization-minimization algorithms in signal processing communications and machine learning[END_REF] performs, at each iteration, an update of the variables by minimizing a surrogate (majorizing function) of the objective, therefore improving the value of the latter.

Using proposition 1 and 2 of [START_REF] Sun | Low-complexity algorithms for low rank clutter parameters estimation in radar systems[END_REF], we obtain the following updates:

• Update {τ k }: at the point τ t k , the updates of the MM algorithm are given for all k ∈ [[

Numerical simulations

To illustrate the performance of proposed estimators, we evaluate their average fraction of energy (AFE) through Monte Carlo simulation. The AFE is considered as a criteria of performance, since, it evaluates the closeness of the true range space UU H towards its estimate U U H . The AFE of a given estimator U is expressed as:

AFE( U) = E{Tr{U H U U H U}} (2.50)
The number of burn-in iterations in the Gibbs sampler N bi =10, N r =200 and N =20, P =5, unless otherwise stated. The matrix S is generated from a compound Gaussian distribution, i.e., its columns are independant, in which, s k d

=

√ τ k d. The texture parameter {τ k } follows the Gamma distribution with unitary mean, i.e., τ k ∼ CG(ν, 1 ν ), ∀k and λ = 1. The matrix N is generated from a Gaussian distribution with zero-mean and covariance matrix σ 2 I. We denote the signal to noise ratio SNR=log( 1 σ 2 ). Finally, for sake of clarity, we assign to the basis of interest U the CIB distribution [START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF][START_REF] Mardia | Directional Statistics[END_REF], i.e., p B (U) ∝ etr{κU H Ū ŪH U} where κ (the concentration parameter) and Ū (randomly initiated) are the prior knowledge which we have about U.

Specically, we generate the matrix data according to Y = US + N, i.e., s k ∼ CG(0, I, τ k ), y k ∼ CN (0, τ k UU H + σ 2 I), ∀k. We compare the following estimators:

• U CG MMSD : the proposed MMSD estimator.

• U Benchmark MMSD : the proposed MMSD estimator with known texture realization τ .

• U U MMSD : the MMSD estimator for uniformly distributed sources of [START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF], dened as U MMSD-CV = P P {κ Ū ŪH + 1 2σ 2 YY H }.

• U SCM : the estimator built from SVD of the SCM, i.e., U SCM = P P {YY H }.

• U MLE : the estimator from [START_REF] Breloy | Maximum likelihood estimation of clutter subspace in non homogeneous noise context[END_REF], corresponding to the MLE estimator of the considered context, while assuming no prior distribution on the subspace basis U. From the analysis of these gures, the following conclusions are drawn:

• The commly used SCM estimator U SCM shows poor performances for low SNR and for small number of samples. As expected, Bayesian estimators perform better in these severe contexts since they include prior knowledge about the basis subspace U in the estimation process. mators.

• The MMSD estimator U U MMSD [START_REF] Besson | Minimum mean square distance estimation of a subspace[END_REF] shows its interest outside the asymptotic regimes. However, the proposed MMSD estimator U CG MMSD outperforms the U U MMSD in the critical cases thanks to the inclusion of the CG model for the signal sources.

Finally, we exemplify briey the use of the MMSD in the context of STAP low rank detection (the reader is referred to our paper [START_REF] Ben Abdallah | Minimum Mean Square Distance Estimation of Subspaces in presence of Gaussian sources with application to STAP detection[END_REF] for more details.) Fig. 2.7 displays the output of the dierent detectors (i.e., the value of the ANMF, adaptive normalized matched lter, for a grid of target parameters) for a conguration where the tested cell contains 10 targets and the number of secondary data is limited.

It is noticed that, the detector which is based on the MMSD provides a trade o that allows for interference rejection, reduced false alarm and reliable target detection. This detector even outperforms the state-of-the-art SFPE (the shrinkage xed point estimator for which the clutter subspace projector is estimated from the SVD of the shrinkage xed point estimator with an optimal regularization parameter) [START_REF] Pascal | Generalized robust shrinkage estimator and its application to stap detection problem[END_REF].

Robust parameterized mean estimation without secondary data

In array processing, there exists in the literature an abundant works to investigate algorithms for parameterized mean estimation and their associated performances, as DoA and target localization

[5355], mostly under the umbrella of Gaussian noise, clutter or environment, depending on the application. For example, in the radar context, the validity of the Gaussian clutter assumption is rooted in the central limit theorem and is realistic in the case of suciently large number of independent and identically distributed (i.i.d.) elementary scatterers. In applications of highresolution radars, the radar clutter exhibits non-stationarity, and a Gaussian modeling of the clutter, be it white or colored, deviates heavily from the real data and thus is inadequate [START_REF] Gini | Performance analysis of two adaptive radar detectors against non-Gaussian real sea clutter data[END_REF][START_REF] Breloy | Robust subspace clustering paradigm for stap[END_REF].

While the discussion, analysis and derivation remains valid whatever the chosen application, in the following, and for sake of clarity and space limitation, we focus on the MIMO radar applications. Such derivation and analysis can be easily applied to any array processing problem under non-Gaussian environment.

Back to the MIMO radar context, out of the many models proposed to deal with these non-Gaussian clutter cases, the CG distribution has, thanks to its ability to describe dierent scales of the clutter roughness and to incorporate various non-Gaussian distributions, become most widely used one [START_REF] Gini | Performance analysis of two adaptive radar detectors against non-Gaussian real sea clutter data[END_REF]5759]. The texture and speckle, of a CG distributed r.v., account for the local temporal power variation of the observations and the scattering received in the antennas, respectively. Though abounding works have been dedicated to the estimation algorithms in the CG clutter context with zero mean observations (the so-called secondary data in the radar community) [6063], there are, to the best of our knowledge, extremely few works dealing with unknown parameterized mean and unknown covariance matrix in the context of CG clutter

In [START_REF] Wang | Maximum likelihood estimation of compound Gaussian clutter and taget parameters[END_REF] and [START_REF] Akcakaya | Adaptive MIMO radar design and detection in compound-Gaussian clutter[END_REF], the authors devised parameter-expanded expectation-maximization (PX-EM) algorithms to estimate the signal as well as clutter parameters for the traditional phased-array radar and MIMO radar, respectively. Nevertheless, the algorithms proposed in [START_REF] Wang | Maximum likelihood estimation of compound Gaussian clutter and taget parameters[END_REF] and [START_REF] Akcakaya | Adaptive MIMO radar design and detection in compound-Gaussian clutter[END_REF] are restricted to a special, linear signal model, called the generalized multivariate analysis of variance (GMANOVA) model [START_REF] Dogandºi¢ | Generalized multivariate analysis of variance: A unied framework for signal processing in correlated noise[END_REF], under which category our context does not fall. To the best of our knowledge, no available algorithm in the current literature addresses the target estimation problem, or the problem of the direction-of-departure/arrival (DoD/DoA) estimation [START_REF] Haardt | Subspace methods and exploitation of special array structures[END_REF] (a highly non-linear problem) in general, under the CG clutter in a comprehensive manner.

In this section, we propose the Iterative Maximum Likelihood Estimator (IMLE) and its Bayesian version, the Iterative MAP Estimator (IMAPE) which are based, respectively, on the conditional likelihood of the observations on the texture realizations, and the joint likelihood between the two. As a consequence, these two estimators are both eo ipso suboptimal. To overcome these algorithms suboptimality, we propose also an iterative ML estimator that is based on the marginal (exact ) observation likelihood, named the iterative marginal ML estimator (IMMLE). Finally, interconnections and comparison between the IMMLE, IMAPE and IMLE is provided.

Model setup 2.3.1.1 Observation model

Consider a MIMO radar system with linear and possibly non-uniform arrays both at the transmitter and the receiver. Further assume that K targets are illuminated by the MIMO radar, all modeled as far-eld, narrowband, point sources [START_REF] Li | MIMO radar Signal Processing[END_REF]. The radar output for the l th pulse in a CPI, and after matched ltering, of the observation Y , in the case of transmission of orthogonal APPLICATION TO RADAR waveforms [START_REF] Haimovich | MIMO radar with widely separated antennas[END_REF], reads

Z(l) = 1 √ T Y (l)S H = K k=1 √ T α k e 2jπf k l a (R) θ (R) k a T (T ) θ (T ) k + N (l), for l = 0, . . . , L -1 (2.51)
where L denotes the number of radar pulses per CPI; α k and f k denote a complex coecient proportional to the radar cross section (RCS) and the normalized Doppler frequency of the k th target, respectively; T is the number of snapshots per pulse, θ

a (T ) (θ (T ) k ) = [e j 2π sin θ (T ) k λ d (T ) 1 , . . . , e j 2π sin θ (T ) k λ d (T ) M ] T and a (R) (θ (R) k ) = [e j 2π sin θ (R) k λ d (R) 1 , . . . , e j 2π sin θ (R) k λ d (R) N ] T ,
in which M and N represent the number of sensors at the transmitter and the receiver, respectively; d (T ) i and d (R) i denote the distance between the ith sensor and the reference sensor for the transmitter and the receiver, respectively; λ stands for the wavelength; N (l) denotes the received clutter matrix at pulse l; and (•) T denotes the transpose of a matrix.

By stacking the output in (2.51) into an M N × 1 vector denoted by z(l), we further have:

z(l) = vec {Z(l)} = A (θ) v(l) + n(l), l = 0, . . . , L -1, (2.52)
in which vec{•} stands for the vectorization of a matrix,

A (θ) = a θ (T ) 1 , θ (R) 1 , . . . , a θ (T ) K , θ (R) K
denotes the steering matrix after matched ltering, where θ = θ

(T ) 1 , θ (R) 1 , θ (T ) 2 , . . . , θ (R) K
is a parameter vector introduced to incorporate all the unknown DoDs and DoAs of the targets, and, a θ

(T ) k , θ (R) k = vec a (R) θ (R) k a T (T ) θ (T ) k = I M ⊗ a (R) θ (R) k a (T ) θ (T ) k
, in which I M stands for the identity matrix of size M , and we recall that ⊗ denotes the Kronecker product, v(l) =

√

T α 1 e 2jπf 1 l , . . . , √ T α K e 2jπf K l T and n(l) = vec {N (l)} denotes the clutter vector after matched ltering at pulse l.

Observation statistics

We model the clutter vectors m(l), l = 0, . . . , L -1 as i.i.d. CG random variable, which can be formulated as the product of two components statistically independent of each other: n(l) = τ (l)x(l), l = 0, . . . , L -1, in which the texture terms τ (l), are i.i.d. positive random variables; the speckle terms x(l) are i.i.d. M N -dimensional circular complex Gaussian vectors with zero mean and second-order moments E x(i)x H (j) = δ ij Σ where Σ denotes the speckle covariance matrix, δ ij is the Kronecker delta. To avoid the ambiguity in the model arising from scaling the texture or the speckle, we assume that tr{Σ} = M N , in which tr{•} denotes the trace.

In this section, we mainly focus on two kinds of CG clutters that are prevalent in the literature, namely, the K-distributed and the t-distributed clutters. Extension to other distribution of the family of CG is straightforward. We recall that, in both cases the texture is characterized by two parameters, the shape parameter a and the scale parameter b:

• K-distributed clutter, in which τ (l) follows a gamma distribution (denoted by τ (l) ∼ Gamma(a, b)), namely, p(τ (l); a, b) = 1 Γ(a)b a τ (l) a-1 e -τ (l)
b , in which Γ(•) denotes the gamma function.

• t-distributed clutter, in which τ (l) follows an inverse-gamma distribution (denoted by

τ (l) ∼ Inv-Gamma(a, b)), thus, p(τ (l); a, b) = b a Γ(a) τ (l) -a-1 e -b τ (l) .

Unknown parameter vector and likelihood function

Under the assumptions above, the unknown parameter vector of our problem is given by:

ξ = θ T , {α} T , {α} T , f T , ζ T , a, b T , (2.53)
in which α = [α 1 , . . . , α K ] T is a complex vector parameter including the RCS coecients of all K targets, f = [f 1 , . . . , f K ] T contains the normalized Doppler frequencies of the targets, ζ is a M 2 N 2 -element vector containing the real and imaginary parts of the entries of the lower triangular part of Σ, {•} and {•} denote the real and the imaginary part, respectively. Let z = z T (1), ..., z T (L -1)

T denote the full observation vector after matched ltering, and τ = [τ (0), . . . , τ (L -1)] T represents the vector of texture realizations at all pulses. The full observation likelihood conditioned on τ can be written as:

p z|τ ; ξ = L-1 l=0 exp -ρ(l) 2 τ (l)
|πΣ| τ M N (l) ;

(2.54)

in which ξ = θ T , {α} T , {α} T , f T , ζ T T
is the unknown parameter vector that does not contain the texture parameters a and b and

ρ(l) = Σ -1/2 (z(l) -A (θ) v(l)) , (2.55)
which represents the clutter realization at pulse l with its speckle spatially whitened. The conditional likelihood in (2.68), multiplied by p(τ ; a, b), leads to the joint likelihood between z and τ :

p (z, τ ; ξ) = p z|τ ; ξ p(τ ; a, b) = L-1 l=0 exp -ρ(l) 2 τ (l) |πΣ| τ M N (l) p(τ (l); a, b).
(2.56)

Finally, the full observation marginal (exact) likelihood, w.r.t. ξ, is obtained by integrating out τ from the joint likelihood in (2.70), as:

p (z; ξ) = +∞ 0 p (z, τ ; ξ) dτ = L-1 l=0 +∞ 0 exp -ρ(l) 2 τ (l) |πΣ| τ M N (l) p(τ (l); a, b)dτ (l).
(2.57)

Iterative marginal maximum likelihood estimator

The derivation procedure of the IMMLE is presented in this section. To begin with, let Λ denote the marginal log-likelihood (LL) function, which is obtained from (2.57), as:

Λ = ln p (z; ξ) = -LM N ln π -L ln |Σ| + L-1 l=0 ln g M N ρ(l) 2 , a, b , (2.58) APPLICATION TO RADAR in which g M N ρ(l) 2 , a, b = +∞ 0 exp -ρ(l) 2 τ (l) τ M N (l) p(τ (l); a, b)dτ (l) =                2 ρ(l) a-M N K a-M N 2 ρ(l) /b 1 2 b M N +a 2 Γ(a) , K-distributed clutter, b a Γ(M N + a) Γ(a) ρ(l) 2 + b M N +a , t-distributed clutter, (2.59)
where K n (•) is the modied Bessel function of the second kind and order n.

We rst devise estimates for the clutter parameters, i.e., of the speckle covariance matrix Σ, and the texture parameters a and b using an iterative optimization approach. Let Σ denote the estimate of Σ when all the remaining unknown parameters are xed, which can be obtained by solving the equation ∂Λ/∂Σ = 0, as

Σ = 1 L L-1 l=0 h M N ρ(l) 2 , a, b • (z(l) -A (θ) v(l)) (z(l) -A (θ) v(l)) H , (2.60)
in which

h M N ρ(l) 2 , a, b = - ∂g M N ( ρ(l) 2 ,a,b) ∂ ρ(l) 2 g M N ρ(l) 2 , a, b =              K a-M N -1 2 ρ(l) /b 1 2 b 1 2 ρ(l) K a-M N 2 ρ(l) /b 1 2 , K-distributed clutter, M N + a ρ(l) 2 + b , t-distributed clutter. (2.61)
We further need to normalize Σ to fulll the assumption that tr{Σ} = M N . Let Σn denote the normalized estimate Σ, which is:

Σn = M N Σ tr Σ .
(2.62)

Similarly, the estimates of a and b when remaining unknown parameters in ξ are xed, denoted by â and b, can be found by equating ∂Λ/∂a and ∂Λ/∂b to zero, respectively, i.e., by solving numerically:

∂Λ ∂a = L-1 l=0 w M N ρ(l) 2 , a, b g M N ρ(l) 2 , a, b = 0 and ∂Λ ∂b = L-1 l=0 k M N ρ(l) 2 , a, b g M N ρ(l) 2 , a, b = 0, (2.63)
w.r.t. a and b, respectively, in which

w M N ρ(l) 2 , a, b = ∂g M N ρ(l) 2 , a, b ∂a =                          - 1 b a Γ(a) +∞ 0 exp - ρ(l) 2 τ (l) - τ (l) b τ (l) -M N +a-1 • ln b τ (l) + Ψ(a) dτ (l), K-distributed clutter, - b a Γ(M N + a) ln ρ(l) 2 b + 1 -Ψ(M N + a) + Ψ(a) Γ(M N ) ρ(l) 2 + b M N +a , t-distributed clutter, ( 2
k M N ρ(l) 2 , a, b = ∂g M N ρ(l) 2 , a, b ∂b =                      1 b a+2 Γ(a) +∞ 0 exp - ρ(l) 2 τ (l) - τ (l) b • τ (l) -M N +a-1 • (τ (l) -ab) dτ (l), K-distributed clutter, - ab a-1 Γ(M N + a) -a ρ(l) 2 + M N b Γ(a + 1) ρ(l) 2 + b M N +a+1
, t-distributed clutter.

(2.65)

Next, we consider the estimate v(l), by solving ∂Λ/∂v(l) = 0, which reads

v(l) = ÃH (θ) Ã (θ) -1 ÃH (θ) z(l), (2.66)
for à (θ) = Σ -1 2 A (θ) , and z(l) = Σ -1/2 z(l), representing the steering matrix and the observation at pulse l, both pre-whitened by the speckle covariance matrix Σ, respectively.

As the expressions in Eqs. (2.60)-(2.66) suggest, the estimation of each of the individual parameters a, b, Σ and v(l) requires the knowledge of all the others parameters, and furthermore the knowledge of the parameter vector θ. This mutual dependence between the unknown parameters makes it impossible to concentrate the LL function (2.58) analytically, i.e., to obtain a closed-form expression for the LL function concentrated w.r.t. each of the aforementioned parameters that are independent of the other ones. Instead, we resort to the so-called stepwise numerical concentration approach [START_REF] Pesavento | Maximum-likelihood direction-of-arrival estimation in the presence of unknown nonuniform noise[END_REF]. This approach consists in concentrating the LL function iteratively, by assuming that certain parameters are known from the previous iteration. For the task under consideration, we assume, at each iteration, that Σ, â and b are known and use them to compute v(l), which is then used in turn to update the values of Σ and â and b to be used in the next iteration. This sequential updating procedure is repeated until convergence criteria is met or a maximum iteration number is reached.

Next, we turn to the estimation of θ. The iterative update approach explained above allows us to drop all the constant terms in the LL function (2.58) (including those terms that contain only Σ, a and b as unknown parameters, as these are assumed to be xed at each iteration). Furthermore, by inserting the expression of v(l) in (2.66) into what remains in the LL function, we obtain the estimate of θ, denoted by θ, as:

θ =              arg min θ L-1 l=0 (M N -a) ln P ⊥ Ã(θ) z(l) -ln K a-M N 2 P ⊥ Ã(θ) z(l) b 1 2 , K-distributed clutter, arg min θ L-1 l=0 ln P ⊥ Ã(θ) z(l) 2 + b , t-distributed clutter.
(2.67)

in which P ⊥ Ã(θ) = I M N -Ã(θ) ÃH (θ) Ã(θ) -1 ÃH (θ)
is the orthogonal projection matrix onto the null space of Ã(θ). Finally, the whole procedure of the IMMLE is summarized in Table 2.4.

Iterative maximum likelihood estimator & its Bayesian variant

In this section, we present the procedure algorithm for the IMLE and the IMAPE schemes. The design and derivation of these estimators is obtained using the methodology as those present in Section. 2.3.2. For more details the reader is invited to consul our papers [START_REF] Zhang | Maximum likelihood and maximum a posteriori direction-of-arrival estimation in the presence of SIRP noise[END_REF][START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF]. APPLICATION TO RADAR

We recall that the IMLE is based on the full observation likelihood conditioned on τ , i.e.,

p z|τ (z|τ ; ψ) = T t=1 exp -1 τ (t) ρ H (t)ρ(t) | πτ (t)Σ | ; (2.68)
in which ψ contains all the unknown parameters except τ The latter pdf, leads to the following

IMLE procedures 1. Find θ(i)
by numerically computing:

θ(i) = arg min θ L-1 l=0 1 τ (i) (l) P ⊥ Â(i) (θ) (l) ẑ(i) (l) 2 (2.69) in which • ẑ(i) (l) = Σ(i) n -1 2 z(l), Â(i) (θ) = Σ(i) n -1 2 A(θ) • P ⊥ Â(i) (θ) (l) = I M N -Â(i) (θ) Â(i) (θ) H Â(i) (θ) -1 Â(i) (θ) H • τ (i) (l)
is considered as a weighting factor which removes bad pulses.

Iterative expression of Σ(i+1)

:

Σ(i+1) = M N L L-1 l=0 n(i) (l) n(i) (l) H n(i) (l) H Σ(i) n -1 n(i) (l)
in which the estimated clutter realization reads n(i

) (l) = z(l) -A θ(i) v(i) (l).
3. Σ(i+1) (after normalization) to update:

τ (i+1) (l) = 1 M N n(i) (l) H Σ(i+1) n -1 n(i) (l)
On the other hand, the IMAPE (its Bayesian variant) is based on the joint likelihood between z and τ , viz.:

p z,τ (z, τ ; ξ) = p z|τ (z|τ ; ψ) p τ (τ ; a, b) = T t=1 exp -1 τ (t) ρ H (t)ρ(t) | πτ (t)Σ | p τ (t) (τ (t); a, b).
Σ(i+1) =                                          2 L L-1 l=0 n(i) (l) n(i) (l) H â(i) -M N -1 b(i) + 4 b(i) n(i) (l) H Σ(i) -1 n(i) (l) + â(i) -M N -1 2 b(i) 2 1 2 , K-distributed clutter â(i) + M N + 1 L L-1 l=0 n(i) (l) n(i) (l) H b(i) + n(i) (l) H Σ(i) -1 n(i) (l) , Student's t-distributed clutter in which τ (t) =                            1 2 (a -N -1) b + (a -N -1) 2 b 2 + 4b (z(t) -v(t)) H Σ -1 (z(t) -v(t)) 1 2
,

K-distributed clutter, (z(t) -v(t)) H Σ -1 (z(t) -v(t)) + b a + N + 1 , t-distributed clutter. (2.71) and b =          T t=1 τ (t) T a , K-distributed clutter, T a T t=1 1 τ (t)
, t-distributed clutter.

(2.72)

On the other hand, calculating ∂Λ J /∂a yields:

∂Λ J ∂a =              -T Ψ(a) -T ln b + T t=1 ln τ (t), K-distributed clutter, -T Ψ(a) + T ln b - T t=1 ln τ (t), t-distributed clutter, (2.73)
From (2.73) it turns out that ∂Λ J /∂a = 0 does not allow an analytical expression of the root, thus â, unlike b in (2.72), can only be calculated numerically. Eqs. (2.71) -(2.73) reveal that the estimates of τ (t), a and b are mutually dependent.

Discussions

Remark 1: Let us recall the expression of θ for the conventional ML estimator (CMLE) [START_REF] Vantrees | Detection, Estimation and Modulation Theory: Radar-Sonar Signal Processing and Gaussian Signals in Noise[END_REF],

which treats the clutter as uniform white Gaussian distributed, θCMLE = arg min θ

L-1 l=0 P ⊥ A(θ) z(l) 2 ,
This expression shows that the CMLE considers simply the mean of P ⊥ A(θ) z(l)

2

(the square of the norm of the projection of the observation at pulse l onto the null space of the steering matrix), while the IMLE and IMAPE, as the expression of θIMLE/IMAPE shows, consider the modied arithmetic mean (pre-whitened by the speckle covariance matrix, and weighted by the inverse of the texture realization at each pulse). It is precisely because of this modication that the IMLE and IMAPE gain their advantages in performance over the CMLE. On the other hand, we can see from (2.67) that the proposed IMMLE considers, instead of direct or modied arithmetic mean of the projections, the sum of their logarithms (modied by some algebraic operations), which is equivalent to a geometric mean. This dierence leads the CMLE, IMLE and IMAPE to treat all the pulses equally, whereas the IMMLE focus only on the pulses corresponding to small values of the texture realization. In this sense, we note that geometric mean provide more robustness to heavy tail noise than the arithmetic mean. Due to space limitation, the reader is referred to Table . 2.5 for a concise comparison between the IMMLE, IMAPE and IMLE.

Remark 2: We remark, that the proposed algorithms do not entail the estimation of the RCS coecients α k , and the normalized Doppler frequencies f k , of the targets, but rather only involves estimating the vectors v(l), which are functions of them. Indeed, in applications where the estimation of those parameters are of interest, one can naturally nd the ML or LS estimates of them by respectively equating an adequate cost function to zero, and then complement our algorithm accordingly. This, however, it is not to be discussed in this section due to the space limitation.

Remark 3: Remark 5: In our estimation procedure we have made the assumption that tr{Σ} = M N . In real-life applications, however, it is common that the true value of tr{Σ} is unknown, or that there exists a discrepancy between the assumed and the true tr{Σ}. It can nevertheless be demonstrated that our proposed IMMLE (like the IMLE and IMAPE) is robust against such tr{Σ} misassumptions, in terms that its performance is independent of the assumed value of tr{Σ}. A theoretical proof of this is provided in the Ph.D. thesis of X. Zhang.

Performance analysis

To evaluate the performance of our algorithms, we further derived, in [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF], expressions for the CRB and for its variants, including the extended Miller-Chang bound (EMLB), the modied CRB (mCRB) and the hybrid CRB (HCRB), w.r.t. the target's spacing parameter. We then provided an extended examination of their relationships, and studied the eect of the texture parameters. It is worth mentioning that our derived CRB expressions dier from those derived by Besson & Abramovich in [START_REF] Besson | On the Fisher information matrix for multivariate elliptically contoured distributions[END_REF], since, i) the parameter's distribution are considered unknown in our work (i.e., the shape and scale parameters), furthermore, ii) for the model considered herein, where the size of the unknown signal parameter vector (hence the dimension of the resulting FIM) is large, a block-wise expression for the computation of the CRB w.r.t. the signal DoDs and DoAs is proposed in order to alleviate the cost complexity of the numerical evaluation of Finally, in order to fully characterize the performance analysis, we have investigated the resolvability problem of two closely spaced targets. In the MIMO radar context, a few recent works, e.g., [7278], have addressed this problem. The clutter in these works, however, is unexceptionally modeled as a Gaussian process. In our paper [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF], we have focused on the resolvability problem concerning two (colocated) MIMO radar targets under non-Gaussian clutter (modeled as CG). To be more specic, we have set as its principal aim the solution to the following question: What is, in a colocated MIMO radar context under non-Gaussian clutter, the minimum angular separation (between two closely spaced targets) required, under which these two targets can still be correctly resolved? " No work in the current literature, to the best of our knowledge, has been dedicated to this question, except our preliminary work [START_REF] Zhang | MIMO radar performance analysis under K-distributed clutter[END_REF], in which we approached this problem by numerical means. In [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF], we carried on with what was set out in [START_REF] Zhang | MIMO radar performance analysis under K-distributed clutter[END_REF] and bring it to completion, by proposing an analytical expression as the solution to the question under discussion, and by considering a wider range of clutter distributions.

To approach this question we resort, in a similar way to [7274], to the concept of the resolution limit (RL), which provides the theoretical foothold of our work to characterize the resolvability of two targets. The RL is dened as the minimum distance w.r.t. the parameter of interest (e.g., the DoDs/DoAs or the electrical angles, etc.) that allows distinguishing between two closely spaced sources [START_REF]Statistical resolution limit for multiple signals and parameters of interest[END_REF][START_REF] Smith | Statistical resolution limits and the complexied Cramér Rao bound[END_REF][START_REF] Shahram | On the resolvability of sinusoids with nearby frequencies in the presence of noise[END_REF]. Various approaches have been devised to account for the RL, generally categorized, in view of the respective theories they rest on, into three families:

those based on the mean null spectrum analysis [START_REF] Cox | Resolving power and sensitivity to mismatch of optimum array processors[END_REF], those capitalizing on the detection theory [START_REF] Shahram | On the resolvability of sinusoids with nearby frequencies in the presence of noise[END_REF]8386], and nally, those concerning the estimation theory and exploiting the CRB [START_REF] Smith | Statistical resolution limits and the complexied Cramér Rao bound[END_REF]8790].

Belonging to the family of the third approach, a widely recognized criterion is proposed by Smith [START_REF] Smith | Statistical resolution limits and the complexied Cramér Rao bound[END_REF], according to which two targets are resolvable if the distance between the targets (w.r.t. the parameter of interest) is greater than the standard deviation of the distance estimation. The prevalence of Smith's criterion, over other criteria derived from the estimation theory, e.g., the one proposed in [START_REF] Lee | The Cramér-Rao bound on frequency estimates of signals closely spaced in frequency[END_REF][START_REF] Lee | The Cramér-Rao bound on frequency estimates of signals closely spaced in frequency (unconditional case)[END_REF][START_REF] Dilaveroglu | Nonmatrix Cramér-Rao bound expressions for high-resolution frequency estimators[END_REF], is largely attributable to its merit of taking the coupling between the parameters into account. Moreover, it enjoys generality in contrast to the mean null spectrum approach, as the latter is designed for certain specic high-resolution algorithms and not for a specic signal model itself [START_REF] Kaveh | The statistical performance of the MUSIC and the minimumnorm algorithms in resolving plane waves in noise[END_REF]. Finally, the RL yielded by Smith's criterion is closely related, as recently revealed in [START_REF]Statistical resolution limit for multiple signals and parameters of interest[END_REF], to the class of the detection theory based approach, meaning that these two approaches can in fact be unied. In view of these merits, we focused on the RL in Smith's sense in [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF]. First, we proposed an analytical expression for the angular resolution limit (ARL 1 )

between two closely-spaced targets in a colocated MIMO radar system under CG clutter. As a byproduct, closed-form expressions of the standard CRB w.r.t. the angular spacing have been derived. Furthermore, we provided numerical illustrations to vindicate our expression, as well as to inspect the properties revealed by it, as the impact of the shape and/or scale parameter variation and the type of the considered CG distribution.

Due to space limitation, we refer the reader to our paper [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF] in which the derivation and analysis are presented in details.

Numerical simulations

We consider a MIMO radar comprising M = 3 sensors at the transmitter and N = 4 at the receiver, both with half-wave length inter-element spacing. The DoD and DoA of the rst source are respectively 18 • and 20 [Σ] m,n = σ 2 0.9 |m-n| e j π 2 (m-n) , m, n = 1, . . . , M N , in which σ 2 is a factor to adjust speckle power. Each point of the MSE in the gures is generated by averaging the results of 100 Monte-Carlo trials.

In Figs. 2.8 and 2.9 we plot the SCR = 1

L L-1 l=0 (A(θ)v(l)) H (A(θ)v(l)) E{τ (l)}tr{Σ}
, in which E{τ (l)} is equal to ab for a K-distributed clutter and b/(a -1) for a t-distributed clutter (for a > 1). The IMMLE procedures Initialization i = 0, set â(0) , b(0) to be two arbitrary positive numbers and Σ(0)

n = I M N Step 1 iteration i, calculate θ(i) from (2.67) using â(i) , b(i) and Σ(i) n calculate v(i) (l) from (2.66) using θ(i) , â(i) , b(i
θ(i) , v(i) (l), â(i+1) and b(i+1) Set i ← i + 1
Step 3 Chapter 3

Robust and scalable parametric calibration with application to radio astronomical arrays Advanced radio interferometers, as the existing low frequency array (LOFAR) [START_REF] Van Haarlem | LOFAR: The LOw-Frequency ARray[END_REF] and the future square kilometre array (SKA) [START_REF] Dewdney | The Square Kilometre Array[END_REF], form large sensor arrays, which are constituted of many small antenna elements. As an example, the LOFAR consists of 50 stations, mainly located across the Netherlands. Each station is a closed packed sensor array, composed of at least 96 low-band antennas (30-90MHz) and 48 high-band antennas . Such interferometers oer a large aperture size and deliver large amounts of data in order to reach high performance in terms of resolution, sensitivity and survey speed [START_REF] Dewdney | The Square Kilometre Array[END_REF]. Nevertheless, to meet the theoretical optimal performances of such next generation radio interferometers, a plethora of signal processing challenges must be conducted, among them, calibration, data reduction and image one or more cooperative sources, named calibrator sources. In the radio astronomy context, calibration is commonly treated using the rst approach as we have access to prior knowledge thanks to tables describing accurately the position and ux of the brightest sources [START_REF] Baars | The absolute spectrum of CAS A -an accurate ux density scale and a set of secondary calibrators[END_REF][START_REF] Kimball | Ivezi¢, A unied catalog of radio objects detected by NVSS, FIRST, WENSS, GB6, and SDSS[END_REF]. Meaning that, the characteristics of the calibration sources, i.e., their true/nominal directions and their powers without the eects of the ionosphere nor antenna imperfections, are a priori knowledge which is required to solve such calibration problems [119125]. Based on this knowledge, stateof-the-art calibration algorithms operate mostly in an iterative manner in a mono-wavelength scenario [START_REF] Van Der Veen | Signal processing tools for radio astronomy[END_REF]126130]. For instance, the (weighted) alternating least squares (WALS) approach has been adapted for LOFAR station calibration in [START_REF] Wijnholds | Fish-eye observing with phased array radio telescopes[END_REF][START_REF] Wijnholds | Multisource self-calibration for sensor arrays[END_REF]. Nevertheless, the major drawback of these schemes can be devised into two categories 2. Second, the Gaussianity assumption which is not realistic in the radio astronomy context.

Specically, the presence of outliers has multiple causes, among which i) the radio frequency interferers, which corrupt the observations and are not always perfectly ltered in practice [START_REF] Leshem | Multichannel interference mitigation techniques in radio astronomy[END_REF], ii) the presence of unknown weak sources in the background [START_REF] Kazemi | Robust radio interferometric calibration using the tdistribution[END_REF], iii) the presence of some punctual events as interference due to the Sun or due to strong sources in the sidelobes which can also randomly create outliers [START_REF] Boonstra | Radio frequency interference mitigation in radio astronomy[END_REF]. To the best of our knowledge, the proposed scheme in [START_REF] Kazemi | Robust radio interferometric calibration using the tdistribution[END_REF][START_REF] Yatawatta | Robust radio interferometric calibration[END_REF], represents the only alternative to the existing calibration algorithms based on a Gaussian noise model using a particular noise modeling, namely, a student-t modeling which may not incorporate all these undesirable perturbation phenomena.

In the following, I present two methodologies to overcome theses drawbacks. The rst one is based on a classical array signal processing model and the second one is based on the socalled Jones parametric model which is commonly used in the astronomical community. In the rst part we base our algorithm design on the complexity ecient and the second part is more devoted to the robustness aspect. Note that this Chapter is based on our following articles Furthermore, the sensor noise power estimation considers the presence of non-calibration sources.

For sake of clarity, let us focus on the so-called regime 3 (a.k.a. 3DC for direction dependent distortion regime with a compact set of antennas) where all lines of sight toward a source in the sky cross the same ionospheric layer and where the thickness of the ionosphere can be direction dependent [START_REF] Lonsdale | Calibration approaches[END_REF], which is represented in Fig. 3.1 and well adapted for the calibration of a LOFAR station and the future SKA stations as well as the core of these arrays. Consequently, in this regime, the ionospheric phase delays modify the geometric delays and introduce angular-shifts for the source directions [START_REF] Van Der Veen | Signal processing tools for radio astronomy[END_REF][START_REF] Cotton | Beyond the isoplanatic patch in the VLA Low-frequency Sky Survey[END_REF], which are direction and wavelength dependent [START_REF] Thompson | Interferometry and Synthesis in Radio Astronomy[END_REF][START_REF] Cohen | Probing Fine-Scale Ionospheric Structure with the Very Large Array Radio Telescope[END_REF]. By estimating calibrator shifts (i.e., the dierence between the true calibrator directions, known from tables [START_REF] Baars | The absolute spectrum of CAS A -an accurate ux density scale and a set of secondary calibrators[END_REF][START_REF] Bennett | The revised 3C catalog of radio sources[END_REF][START_REF] Kimball | Ivezi¢, A unied catalog of radio objects detected by NVSS, FIRST, WENSS, GB6, and SDSS[END_REF], and their estimated apparent directions), interpolation methods can be eciently applied in order to obtain a phase screen model, that captures the ionospheric delays over the entire Field-of-View [START_REF] Cotton | Beyond the isoplanatic patch in the VLA Low-frequency Sky Survey[END_REF]. We emphasize that in addition to the phase screen reconstruction step, the calibration usually involves the estimation of the complex undirectional gains of the antennas, their directional gains toward each calibrator and their noise powers [START_REF] Wijnholds | Fish-eye observing with phased array radio telescopes[END_REF],

for the whole available range of wavelength range, i.e., processing bandwidth.

Model setup a

Consider an array comprised of P elements, with known locations, each referred by its Cartesian coordinates ξ p = [x p , y p , z p ] T for p = 1, . . . , P , that we stack in Ξ = [ξ 1 , . . . , ξ P ] T ∈ R P ×3 . This array is exposed to Q known strong calibration sources and Q U unknown weak non-calibration sources.

Let D K = d K 1 , . . . , d K Q ∈ R 3×Q and D U = d U 1 , . . . , d U Q U ∈ R 3×Q U
denote the known (true/nominal ) calibrator direction cosines and unknown non-calibrator direction cosines, respectively, in which each source direction d = [l, m, n] T can be uniquely described by a couple (l, m),

since n = √ 1 -l 2 -m 2 .
The ionosphere introduces an unknown angular-shift for each source direction, depending on the wavelength λ, which is related to the frequency f = c λ , with c denoting the light speed. Consequently, we distinguish between the unknown apparent directions w.r.t. the calibrators, denoted by D λ = [d λ,1 , . . . , d λ,Q ], and their true/nominal known directions D K , i.e., without the propagation disturbances.

Under the narrowband assumption, the steering vector a λ (d) toward the direction d at wavelength λ is given by

a λ (d) = a λ (l, m) = 1 √ P exp -j 2π λ Ξd , (3.1)
Figure 3.1: The so-called regime 3, which is considered in this section, assumes that V S and A S. This leads to ionospheric perturbations which are direction dependent (after [START_REF] Van Der Tol | Self-calibration for the LOFAR radio astronomical array[END_REF][START_REF] Lonsdale | Calibration approaches[END_REF]).

that we gather for multiple directions in the steering matrix

A D λ = 1 √ P exp -j 2π λ ΞD λ . (3.2)
We assume that all antennas have identical directional responses. Their directional gain responses (and propagation losses) are modeled by two diagonal matrices,

Γ λ ∈ C Q×Q and Γ U λ ∈ C Q U ×Q U
, toward the calibration and non-calibration sources, respectively.

The received signals from each antenna are split into narrow subband and stacked in a vector, leading to

x λ (n) = G λ A D λ Γ λ s λ (n) + A D U λ Γ U λ s U λ (n) + n λ (n), (3.3)
for the n-th observation and wavelength λ, where G λ = diag(g λ ) ∈ C P ×P models the undirec-

tional antenna gains, s λ (n) ∈ C Q and s U λ (n) ∈ C Q U represent,
Σ λ = diag (σ λ ) ∈ R Q×Q , Σ U λ = diag (σ U λ ) ∈ R Q U ×Q U and Σ n λ = diag (σ n λ ) ∈ R P ×P be the diago-
nal covariance matrices for the calibrators, non-calibration sources and sensor noises, respectively.

Consequently, the covariance matrix R λ = E x λ x H λ of the observations corresponding to model

(3.3) is given by R λ = E D λ M λ E H D λ + R U λ + Σ n λ , (3.4)
in which

E D λ = G λ A D λ Σ 1 2 λ , (3.5)
M λ = Γ λ Γ H λ = diag (m λ ) , (3.6)
3.1. SPARSE AND PARALLEL MULTI-WAVELENGTH CALIBRATION ALGORITHM 63 and where we have dened the unknown covariance matrix for the non-calibration sources as

R U λ = G λ A D U λ Γ U λ Σ U λ G λ A D U λ Γ U λ H . (3.7)
To the best of our knowledge, all algorithms present in the literature, ignore the contribution of the non-calibrator sources and thus, use the following model

R λ = E D λ M λ E H D λ + R U λ + Σ n λ = = E D λ M λ E H D λ + Σ n λ (3.8)
In order to overcome the scaling ambiguities in the observation model (3.8), we consider the following commonly used assumptions in radio astronomy: i) to resolve the phase ambiguity of g λ , we take its rst element as the phase reference; and ii) m λ shares a common scalar factor with g λ and when solving for the calibrator directions, a common rotation of all steering vectors can be compensated by the undirectional gain phase solution. We therefore consider an additional know source toward the zenith to remove both ambiguities, by xing its direction d 0 and its directional gain.

3.1.2 Model eects of the wavelength on antenna gains, source direction shifts and source powers a

In the radio astronomy context, it is commonly adopted that the antenna and source parameters are wavelength dependent. Consequently, we assume smooth or/and known variation of the parameters g λ , Γ λ , Σ λ , Σ U λ and Σ n λ over λ, as commonly used in recent astronomy applications.

We summarize the particular behavior of the underlying parameters as follows:

C1) The undirectional gains, g λ , vary smoothly over λ. Common models for characterizing these behaviors consist of classical polynomials of power law over λ.

C2) The directional gains, Γ λ , are inversely proportional to λ, i.e., Γ λ ∝ λ -1 , as observed in practice. Note that the proposed algorithm can be straightforwardly adapted with another given behavior (including the extreme cases: a) a constant behavior across the frequency range and b) statistical independence between frequency bins).

C3) As a consequence of the ionospheric delays, that are at the origin of the directional shifts, the shifts are proportional to λ 2 .

C4) The source powers, Σ λ and Σ U λ , vary commonly with a power law with dierent spectral indexes. We consider the calibrator powers, Σ λ , to be known from tables.

C5)

p λ = g T λ , d T λ,1 , . . . , d T λ,Q , m T λ , σ nT λ T from J sample covariance matrices Rλ = 1 N N n=1 x λ (n)x H λ (n) λ∈Λ , (3.10)
where Λ = {λ 1 , . . . , λ J } represents the set of the J available wavelengths for the whole network. In this subsection, we dene the main steps of our proposed algorithm, then, in the following subsections, we describe each step.

Algorithm 1: Parallel Multi-wavelength Calibration Algorithm Input:

Rλ , m [0] λ λ∈Λ , D K , η p ; 1 Init: set i = 0, g λ = g [0] λ λ∈Λ , D λ = D K , m λ = m [0] λ , Ω λ = 1 P ×P λ∈Λ ; 2 repeat 4 4 i = i + 1; 6 6 Estimate in parallel g [i] λ λ∈Λ with Algorithm. 1.1; 8 8 Estimate in parallel D [i] λ , m [i] λ , σ n[i] λ λ∈Λ
with Algorithm. 1.2; 10 10 Update locally

Ω [i] λ = σ n[i] λ σ n[i]T λ -1 2 λ∈Λ ; 11 until p [i-1] -p [i] 2 ≤ p [i] 2 η p ; Output: p = p [i]T λ 1 , . . . , p [i]T λ J T ;
Let us note that a large number of samples, statistically ecient estimators can be devised using the WLS approach. In this sense, we dene, for each λ ∈ Λ, the local cost function [START_REF] Wijnholds | Multisource self-calibration for sensor arrays[END_REF] to be minimized as

κ λ (p λ ) = W -1 2 λ R λ (p λ ) -Rλ W -1 2 λ 2 F , (3.11)
in which

R λ (p λ ) = E D λ M λ E H D λ + Σ n λ (3.12)
denotes the covariance matrix when the contribution of the (weakest) non-calibrators is ignored, and W λ is the weighting matrix. The optimal weighting matrix for Gaussian noise is the inverse of the covariance of the residuals [START_REF] Ottersten | Covariance matching estimation techniques for array signal processing applications[END_REF], which is generally unknown. Justied by physical reasons, we consider W λ = I in our alternating algorithm as an initialization and rene it as W λ = Σ n λ once we obtain an estimate of Σ n λ . Since Σ n λ is diagonal, we rewrite the local cost function (3.11) as

κ λ (p λ ) = R λ (p λ ) -Rλ Ω λ 2 F (3.13) Ω λ = σ n λ σ nT λ - 1 2 .
( (3.15)

Our aim is to estimate p by minimizing κ(p) in an alternating and parallel manner. Note that the overall problem is non-convex, so we can not claim nding the global minimum when the algorithm converges (neverthless, in practice we might have good initialisation that converges to the global minimum as demonstrated by our simulations and application to real data). We rst estimate locally in each agent the parameter vector {g λ } λ∈Λ , with the remaining parameters in p xed as described in Section 3.1.4.2, by reformulating the problem as a consensus problem, in which the updates of the consensus parameter α, that is explained below, provide coherence among wavelength. In a second step, we estimate the variables {m λ , d λ,1 , . . . , d λ,Q , σ n λ } λ∈Λ for xed {g λ } λ∈Λ , by using a sparse representation approach as described in Section 3.1.4.3. Finally, we update the weighting matrices {Ω λ } λ∈Λ . During these procedures, the amount of information that needs to be exchanged between the fusion center and the compute agents is much less than the volume of data being calibrated, making this scheme computationally feasible. The overall procedure, referred to as PMCA, is presented in the box Algorithm 1. The algorithm is carefully initialized with the true/nominal calibrator parameters and an initial guess for the antenna gains, e.g. from preceding calibration, or by default by the unit sensor gain, g

[0] λ = 1. The stopping criterion η p has to be suciently small to assure convergence. In the following subsections, we present the two major alternating optimization steps of the proposed PMCA.

Direction independent antenna gain estimation

In this subsection, we describe Algorithm. Then, each agent transfers its parameter estimates to the fusion center. At the fusion center, smoothness of the parameters across wavelength is enforced (line 3 of Algorithm.1.1), using the coherence parameter vector α, that is passed to each agent to provide coherent processing across the whole wavelength range, and thus improving the local calibration procedure 1 .

With this network setup, we formulate the parallel calibration procedure as

minimize α,{g λ } λ∈Λ λ∈Λ κ[i] λ (g λ ) s.t. g λ = B λ α, ∀λ ∈ Λ, (3.16)
in which

κ[i] λ (g λ ) = κ λ (p λ |m [i] λ , D [i] λ , σ n[i] λ ), i is the [i]-th iteration of Algorithm 1, B λ = b λ ⊗ I represents the coecient basis with b λ = [b 1,λ , . . . , b K,λ ] T , α = [α 1,1 , . . . , α 1,P , α 2,1 , . . . , α K,P ] T
denotes the augmented vector of hidden variables and where the cost function consists of a sum of independent cost functions, one for each subband, that are coupled through the coherence constraints which however are independent across sensors. A commonly way to solve (3.16)

L [i] {g λ } λ∈Λ , α, {y λ } λ∈Λ = λ∈Λ κ[i] λ (g λ ) + y H λ (g λ -B λ α) + ρ 2 g λ -B λ α 2 2 = λ∈Λ L [i] λ (g λ , α, y λ ) , (3.17)
where {y λ } λ∈Λ are the J Lagrange multipliers and ρ is the regularization term, chosen by the practitioner. Thus, denoting t the local iteration counter of Algorithm.1.1, the updates for the [t]-th iteration are given by g

[t]
λ = arg min

g λ L [i] λ g λ , α [t-1] , y [t-1] λ , λ ∈ Λ, (3.18
)

α [t] = arg min α λ∈Λ L [i] λ g [t]
λ , α, y

[t-1] λ , ( 3.19
) Minimization of (3.18) To minimize (3.18), we follow an iterative optimization approach based on [START_REF] Salvini | Fast gain calibration in radio astronomy using alternating direction implicit methods: Analysis and applications[END_REF][START_REF] Goldstein | Fast alternating direction optimization methods[END_REF], that we adapt to our parallel minimization. Let us assume that g λ and g * λ are two independent variables.

y [t] λ = y [t-1] λ + ρ g [t] λ -B λ α [t] , λ ∈ Λ, (3.20
We then regard g * λ as xed and minimize L

[i] λ g λ , g * λ , α [t] , y

[t] λ w.r.t. g λ only, and without considering the diagonal elements in the cost functions in (3.16) that contain the unknown noise variances σ n λ . Let us dene the operator S p (.), that converts the p-th row of a matrix to a vector and removes the p-th element of this selected vector. Further, dene the vector rp λ = S p Rλ , the weighting vector ω = S p Ω

[i] λ , z = S p R K λ in which the estimated calibrator sky model is R K λ = A D λ M λ A H D λ
. Consequently, after some manipulations, this leads to

[ĝ

[t+1] λ ] p = 2z H ω rp λ ω + ρ B λ α [t] p -[y [t] λ ] p 2 (z H ω z ω ) + ρ , ( 3.21)
where

z ω = z S p Ω [i] λ . Then, we directly update [g * [t+1] λ ] p = [ĝ [t+1] λ ] p *
and proceed in the same manner with the remaining parameters in g λ . This procedure is summarized in Algorithm 1.1a and is repeated until convergence.

Minimization of (3.19)

After gathering the estimates {g By taking into account the aforementioned assumptions, i.e., the resolution grid are ∆l q λ ∝ λ 2 , ∆m q λ ∝ λ 2 , and by scaling them around λ 0 , our grid resolution reads

∆l q λ = λ λ 0 2 ∆l q λ 0 , ( 3.25
)

∆m q λ = λ λ 0 2 ∆m q λ 0 . (3.26)
Thus, we dene JQ dictionaries of steering vectors, Ãq,λ ∈ C P ×Nq , for q = 1, . . . , Q, λ ∈ Λ, which contain N q steering vectors, centered around the true/nominal direction of the q-th calibrator, namely d K q , with resolution ∆l q λ , ∆m q λ and N q 1. These dictionary steering matrices are gathered in

Ãλ = Ã1,λ , . . . , ÃQ,λ ∈ C P ×Ng , λ ∈ Λ, (3.27)
with N g = Q q=1 N q denoting the total number of directions on the grid. We dene then J vectors, { mλ } λ∈Λ , as mλ = mT 1,λ , . . . , mT

Q,λ T ∈ R Ng , λ ∈ Λ, (3.28)
which contain the squared DD gains towards all calibrators, where mq,λ is the sparse vector associated with Ãq,λ . Due to the assumption of non-overlapping displacement sectors, each mq,λ is exactly 1-sparse, i.e., mq,λ 0 = 1, for q = 1, . . . , Q, ∀λ ∈ Λ. Since the shift resolution in the dictionaries is made proportional to λ 2 , the support of mq,λ is independent of λ. To proceed, we exploit that Γ λ ∝ λ -1 (see Section 3.1.2) in order to estimate a unique sparse vector for all wavelengths, namely m. More precisely, under this assumption, we dene m = mT 1 , . . . , mT

Q T as mλ = λ 0 λ 2 m, ∀λ ∈ Λ, (3.29)
which can be also adjusted for other existing models of Γ λ . Consequently, we formulate the initial minimization problem as

minimize m,{σ n λ } λ∈Λ λ∈Λ rλ -V λ m -N λ σ n
V λ = (Σ n λ ) -1 2 Ẽ * λ ⊗ (Σ n λ ) -1 2 Ẽλ ,N λ = Σ n-1 2 λ • Σ n-1 2 λ ,r λ = vec Rλ Ω λ , in which Σλ = blkdiag I N 1 ×N 1 [σ λ ] 1 , . . . , I N Q ×N Q [σ λ ] Q and Ẽλ = λ 0 λ G λ Ãλ Σ 1 2 λ .
To handle the above constraints, which are non-convex, we choose the Distributed Iterative Hard Thresholding method [145147], which is based on Iterative Hard Thresholding [START_REF] Blumensath | Iterative Hard Thresholding for compressed sensing[END_REF].

This greedy algorithm consists of a projected gradient descend direction algorithm and oers strong theoretical guarantees that have been successfully employed in the DoA estimation context [START_REF] Ollila | Robust iterative hard thresholding for compressed sensing[END_REF][START_REF] Ollila | Multichannel sparse recovery of complex-valued signals using huber's criterion[END_REF]. Particularly, when the grid is ne and the columns of Ãq,λ are strongly coherent, we can guarantee that each estimated mq is exactly 1-sparse. Thus, using the Coordinate Descent algorithm [START_REF] Friedman | Pathwise coordinate optimization[END_REF], we obtain an analytic solution for each sub-problem and the sparsity of the desired minimizer m reduces the computational complexity. Each step involves the hard thresholding operator H s (.), that keeps the s-largest components of a vector and sets the remaining entries equal to zero, thus, it automatically satises both constraints of sparsity and positivity. Let us denote rλ and V λ , that refer, respectively, to rλ and V λ , where the elements corresponding to the diagonal of R λ are discarded. Then, we obtain the update for the [k]-th iteration as

m[k] q = H 1 m[k-1] q + τ [k] q λ∈Λ V qT λ rq λ -V q λ m[k-1] q , ( 3.31)
with V q λ the q-th column of V λ , where we can allow a step size τ

[k] q that depends on m[k-1] q and the [k]-th iteration, by the use of the Normalized Iterative Hard Thresholding [START_REF] Blumensath | Normalized Iterative Hard Thresholding: Guaranteed stability and performance, Selected Topics in Signal Processing[END_REF]. In addition, the residual is dened as rq λ = rλ -Q q =1,q =q V q λ mq . Afterward, the estimation of {σ n λ } λ∈Λ is performed locally, without the need of transmitting the estimated values. Firstly, note that without considering non-calibration sources, i.e., R U λ ≈ 0, the estimation of σ n λ is given by

σ n λ = vecdiag Rλ -Ĝλ RK λ ĜH λ , (3.32)
since we assume independence of σ n λ across wavelength. Secondly, we remove the bias introduced by the non-calibration sources as follows: we calculate the power

σr λ = a H λ (d λ,r ) Rλ -Ĝλ RK λ ĜH λ a λ (d λ,r ) a λ (d λ,r ) 2 2 (3.33)
of the residual sample covariance matrix for a random direction d λ,r , where no source is supposed to be present. We then approximate a H λ (d λ,r )a λ (d λ,q ) ≈ 0 for any d λ,r = d λ,q , which yields σr λ as the sum of the sensor noise powers [START_REF] Van Der Veen | Signal processing tools for radio astronomy[END_REF][START_REF] Leshem | Radio-astronomical imaging in the presence of strong radio interference[END_REF]. By imposing

P p=1 [σ n λ ] p = σr
λ to the minimization of (3.13) w.r.t. σ n λ , the new unbiased solution is given by

σn λ = σ n λ + 1 P σr λ -1 T P ×1 σ n λ 1 P ×1 , (3.34)
which concludes the description of Algorithm.1.2.

Simulations

The proposed method is evaluated in realistic situations, with typical parameter values commonly used in radio astronomy applications (with consisting of Q = 2 strong calibration sources and Q U = 8 weak non-calibration sources taken from the ten strongest sources in the table of [START_REF] Bennett | The revised 3C catalog of radio sources[END_REF]). 1 Init: set t = 0, g [t] λ = g

[i-1] λ λ∈Λ , R K λ = A D [i-1] λ M [i-1] λ A H D [i-
Input: Rλ , R K λ , g [t-1] λ , y [t-1] λ , Ω [i] λ λ∈Λz , α [t-1] , η g ; 1 Init: set t λ = 0, g [t λ ] λ = g [t-1] λ λ∈Λz ; 2 foreach λ ∈ Λ z do 3 repeat 5 5
t λ = t λ + 1 ; 6 for p = 1, . . . , P do 8 8

Estimate [g

[t λ ] λ ] p using (3.21) ; 9 end 10 until g [t λ -1] λ -g [t λ ] λ 2 ≤ g [t λ ] λ 2 η g ; 11 end Output: ĝλ = g [t λ ]
λ λ∈Λz

;

The reader is referred to [START_REF] Brossard | Parallel multi-wavelength calibration algorithm for radio astronomical arraysk[END_REF] for a full and detailed description of the parameter set up which are adapted for the LOFAR's Initial Test Station [START_REF] Wijnholds | Sky noise limited snapshot imaging in the presence of RFI with LOFAR's Initial Test Station[END_REF], see Fig. 3.2.

The statistical performance of the proposed algorithm is compared with mono-calibration scheme, the unconstrained-CRB C λ and the multi-constrained-CRB, C. In Fig. 3.3, we plot the Root Mean Square Error (RMSE) for the estimates of g λ . Results are averaged for 500 Monte-Carlo runs, for each chosen value of K and ρ. We approach the multi-constrained-CRB for K = K K (which denotes the exact order polynomial value) and even with both K = 2 and K = 4 (i.e., with a misspecication in the order polynomial value), we signicantly improve mono-calibration and are close to the CRB. Moreover, we also have errors due to polynomial interpolation, which is clearly seen at the edge wavelengths. In addition, numerical simulations show that the convergence of the proposed iterative algorithm is attained for less then 5 iterations.

On the other hand, during the DoA estimation, we choose initially a coarse grid, with the same resolution for each coordinate of each calibrator. We apply grid renements [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF] until we avoid o-grid mismatch.

In order to investigate the statistical performance, we perform 200 Monte-Carlo runs for different sample sizes N , and setting ρ = 5 • 10 3 P and K to its true value. We plot the RMSE of dierent parameters in Fig. 3.4 and Fig. 3.5, as function of the number of samples N and compare them to their corresponding multi-constrained-CRB. As expected, the method approaches RADIO ASTRONOMICAL ARRAYS

Algorithm 1.2: estimation of D λ , m λ , σ n λ λ∈Λ Input: Rλ , g [i] λ λ∈Λ , p [i-1] , η m; 1 Init: set k = 0, g [k] λ = g [i] λ , M [k] λ = M [i-1] λ λ∈Λ , D [k] λ = D [i-1] λ , σ n[k] λ = σ n[i-1] λ λ∈Λ ; 2 repeat 4 4 k = k + 1; 5 for q = 1, . . . , Q do 6 foreach A z , z = 1, . . . , Z do 7 foreach λ ∈ Λ z do 9 9
Calculate locally rq λ = rλ -Q q =1 q =q V q λ mq ; q with (3.31); 17 17

The fusion center communicates the non-zero element of

m[k]

q and its associated direction d the multi-constrained-CRB. This clearly shows the good robustness of the method in low SNR scenarios to the presence of non-calibration sources. This is mainly due to i) considering simultaneously multi-wavelength observations (thus, more observations) and ii) forcing smoothness constraints which can attenuate the eect of such unmodelled sources considered as outliers.

[k] λ 0 ,q ; 18 end 19 until m[k-1] -m[k] 2 ≤ m[k] 2 η m; 21 21 Deduce locally mλ λ∈Λ from m[k] = m[k]T 1 , . . . , m[k]T Q T , with (3.29)
λ 0 = d [k] λ 0 ,1 , . . . , d [k] λ 0 ,Q with (3.
However, it should be noted that this robustness has its own limits depending on the number and the power of the unmodelled sources.

Let us now turn to the application to LOFAR data. For this test, we use data from a single polarization measurement with a 48-element low band antenna array of a single core station, We select data for 31 frequencies in the range 47.6-53.5 MHz, with central frequency (respectively, wavelength) f 0 = 50.5 MHz (λ 0 = 5.9 m) and choose a polynomial order K = 3. We plot results of the images produced from the calibrated (using the PMCA algorithm) covariance matrix for the frequency channel centered at 49.4 MHz in Fig. 3.8, compared with the state-of-the art [START_REF] Wijnholds | Multisource self-calibration for sensor arrays[END_REF], which considers the apparent direction of calibrators as xed w.r.t. λ. closer look at the two images in Fig. 3.8 reveals that the diuse emission on the South of the sky, which is a bright section of the Galactic plane (which can be seen by eye in the night sky as the Milky Way), appears slightly brighter in the image calibrated using PMCA. This indicates that PMCA is more robust to the presence of unmodelled sources in the image thereby reducing the likelihood that power in the unmodelled sources is tted to the calibrator sources by the calibration algorithm, which causes a bias in the calibration solutions.

Another indication that PMCA provides robustness to unmodelled ux is that the DI gain solutions, as shown in Fig. 3.9 for a few representative examples, are similar to those found with the conventional approach despite the fact that the diuse emission was not spatially ltered as is done in the conventional approach. This comparison also shows that the gain solutions produced by the conventional approach vary quite signicantly from one frequency channel to the next while the PMCA provides smoothed solutions, which are physically more realistic. for a few representative examples between the solutions found using PMCA (solid lines) and the conventional per-channel calibration approach (dots).

Robust calibration of radio interferometers in non-Gaussian environment

First, unlike the previous section, we note that the parametric model used here to describe the perturbation eects is based on the so-called Jones matrices [START_REF] Hamaker | Understanding radio polarimetry[END_REF][START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF]. Such formalism describes in a exible way the conversion of the incident electric eld into voltages and captures some eects that cannot be described by the classical array processing model. Indeed, along its propagation path, the signal is aected by various eects and transformations which correspond to matrix multiplications in the mathematical Jones framework. Multiple distortion eects caused by the environment and/or the instruments can be easily incorporated into the model using an adequate parametrization of the Jones matrices. Such eects can represent, for example, the ionospheric phase delay resulting in angular shifts, the atmospheric distortions, the typical phase delay due to geometric path-length dierence, the voltage primary beam, the Faraday rotation, the crossleakage or also the electronic gains [START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF][START_REF] Noordam | The MeqTrees software system and its use for thirdgeneration calibration of radio interferometers[END_REF]. For the above reasons and due to its exibility [157159, [START_REF] Thompson | Interferometry and synthesis in radio astronomy[END_REF][START_REF] Yatawatta | Radio interferometric calibration using the SAGE algorithm[END_REF], we adopt this parametric model herein. Second, we make a distinction between the non-structured and the structured cases: in the rst case, one total Jones matrix stands for all the eects along the full signal path while in the second case, we regard each physical eect separately thanks to individual Jones terms in a cumulative product. Thus, dierent corruptions are described by dierent kinds of Jones matrices. We emphasize that the proposed algorithm, entitled relaxed concentrated ML estimator, is a generic algorithm as it is based on a nonstructured Jones matrices formulation as a rst step. However, it can be adapted to various regimes describing distinct calibration scenarios in which an array can operate [START_REF] Lonsdale | Calibration approaches[END_REF]. Finally, we consider mainly the specic example of the direction dependent distortion regime with a compact set of antennas, which is referred as regime 3 or 3DC regime (direction dependent distortion regime with a compact set of antennas). The array is therefore considered as a closely packed group of antennas but the array elements have a wide eld of view. This is particularly well-adapted for calibration of compact arrays, typically a LOFAR station.

Data model

Case of non-structured Jones matrices

Let us consider M antennas with known locations that receive D signals emitted by calibrator sources. Each antenna is dual polarized and composed of two receptors, in order to provide sensitivity to the two orthogonal polarization directions (x, y) of the incident electromagnetic plane wave. Consequently, the relation between the i-th source emission and the measured voltage at the p-th antenna is given by [START_REF] Hamaker | Understanding radio polarimetry[END_REF][START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF][START_REF] Ker | Radio AGN evolution with low frequency radio surveys[END_REF] v i,p (θ) = J i,p (θ)s i

(3.35)
where

s i = [s ix , s iy ] T is the incoming signal, v i,p (θ) = [v i,px (θ), v i,py (θ)]
T is the generated voltage with one output for each polarization direction and J i,p (θ) denotes the so-called 2 × 2 Jones matrix, parametrized by the unknown vector of interest θ. The Jones matrix models the array response and all the perturbations introduced along the path from the i-th source to the p-th sensor. Since each propagation path is particular, we can associate a dierent Jones matrix with each source-antenna pair (i, p), leading to a total number of DM Jones matrices. In this section, we consider the non-structured case where no specic perturbation model is used to describe the physical mechanism behind each perturbation eect and the unknown elements correspond to the entries of all Jones matrices [START_REF] Yatawatta | Radio interferometric calibration using the SAGE algorithm[END_REF][START_REF] Nunhokee | Link between ghost artefacts, source suppression and incomplete calibration sky models[END_REF] (a structured example is given for 3DC calibration regime, in Section 3.2.1.

2). RADIO ASTRONOMICAL ARRAYS

For each antenna pair, we compute the correlation of the output signals, resulting in the typical observations recorded by a radio interferometer. The correlation between voltages is given, in the case of noise free measurements, for the (p, q) antenna pair, by

V pq (θ) = D i=1 J i,p (θ)C i J H i,q (
θ) for p < q, p, q ∈ {1, . . . , M },

where the signals emitted by the sources are assumed uncorrelated and the 2 × 2 matrix C i = E{s i s H i } is known from prior knowledge. Let us remark that autocorrelations are not considered as shown by the condition p < q in (3.36) (this is a typical situation in the radio astronomy context where radio interferometric systems automatically ag the autocorrelations [START_REF] Van Haarlem | LOFAR: The LOw-Frequency ARray[END_REF]). We rewrite (3.36) as a 4 × 1 vector

ṽpq (θ) = vec V pq (θ) = D i=1 u i,pq (θ) (3.36) in which u i,pq (θ) = J * i,q (θ) ⊗ J i,p (θ) c i with c i = vec(C i ).
We stack all the noisy measurements within a full vector

x = v T 12 , v T 13 , . . . , v T (M -1)M T ∈ C 4B×1 , where B = M (M -1)
2 denotes the total number of antenna pairs and v pq = ṽpq (θ) + n pq with n pq the noise sample at a specic antenna pair. Specically, x reads

x = D i=1 u i (θ) + n (3.37) in which u i (θ) = u T i,12 (θ), u T i,13 (θ), . . . , u T i,(M -1)M (θ) T and n = n T 12 , n T 13 , . . . , n T (M -1)M
T is the full noise vector which accounts for thermal noise, but also the presence of outliers in our data. Therefore, the noise can no longer be considered Gaussian and a robust calibration method is required. To investigate non-Gaussian noise modeling and encompass a broad range of noise distributions, we propose to adopt the CG noise model [START_REF] Jay | Détection en environnement non gaussien[END_REF][START_REF] Yao | Spherically invariant random processes: Theory and applications[END_REF]. Specically, the noise at each antenna pair is assumed to be generated as

n pq = √ τ pq g pq , (3.38)
where the positive real random variable τ pq is referred to as texture, whereas the complex speckle component g pq follows a zero-mean Gaussian distribution, i.e., g pq ∼ CN (0, Ω).

(3.39)

In order to remove scaling ambiguities, we impose tr {Ω} = 1. Note that the choice of this constraint is arbitrary and does not aect the estimates of interest as argued in [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF].

In this section, we adopted the non-structured Jones matrices formulation which is relevant in the radio astronomical context [START_REF] Yatawatta | Radio interferometric calibration using the SAGE algorithm[END_REF][START_REF] Nunhokee | Link between ghost artefacts, source suppression and incomplete calibration sky models[END_REF]. In this case, there is no need to specify the full propagation path, avoiding misspecication in the model. Besides, it is highly exible and can be adapted to dierent scenarios [START_REF] Lonsdale | Calibration approaches[END_REF]. In the following, we present the direction dependent distortion regime with a compact set of antennas, named 3DC regime or regime 3.

Specic case of the 3DC calibration regime

In this regime, direction dependent distortions play a signicant role since individual elements in the array have a wide eld of view. Indeed, this implies dierent propagation conditions towards

ROBUST CALIBRATION OF RADIO INTERFEROMETERS IN NON-GAUSSIAN ENVIRONMENT 77

distinct sources in the eld of view. However, the array being relatively compact, made of similar elements, some eects might be the same for all antennas. In the following, we introduce a particular sequence of Jones matrices with specic parametrizations, in the context of regime 3 [START_REF] Yatawatta | Reduced ambiguity calibration for LOFAR[END_REF] J

i,p (θ 3DC i,p ) = G p (g p )H i,p Z i,p (α i )F i (ϑ i ) (3.40)
for i ∈ {1, . . . , D}, p ∈ {1, . . . , M } and θ 3DC i,p = [ϑ i , g T p , α T i ] T . We note H i,p the only assumed known matrix thanks to electromagnetic simulations in terms of antenna response and a priori knowledge given by calibrator source and antenna positions [START_REF] Kazemi | Blind calibration for radio interferometry using convex optimization[END_REF][START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF][START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF][START_REF] Yatawatta | Reduced ambiguity calibration for LOFAR[END_REF], whereas the remaining matrices are explained in the following items.

• Ionospheric eects :

Propagation through the ionosphere, the outer layer of the earth's atmosphere, introduces propagation delay on the signal which is aected by spatially variable uctuations. In the case of a compact array, the ionospheric delay matrix is in fact a scalar direction-dependent phase given by Z i,p (α i ) = exp jϕ i,p I 2

(3.41) in which ϕ i,p = η i u p + ζ i v p where α i = [η i , ζ i ]
T is the vector of unknown osets resulting in a shift of the i-th source direction and r p = [u p , v p ] T is the vector of known antenna position in units of wavelength.

On top of that, passing through the ionosphere is associated with a rotation of the polarisation plane of each signal source around the line of sight. We call it the ionospheric Faraday rotation matrix F i (ϑ i ) and write it as

F i (ϑ i ) = cos(ϑ i ) -sin(ϑ i ) sin(ϑ i ) cos(ϑ i ) (3.42)
where ϑ i is the unknown Faraday rotation angle, assumed identical for all antennas, since the array has a limited spatial extent [START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF].

• Instrumental eects : Individual antennas are described by electronic complex gains which appear in G p (g p ) = diag{g p } with g p the unknown electronic gain vector. Therefore, in this specic structured case, the physical model parameters in (3.40) are collected in the vector ε

3DC = P[θ 3DC T 1,1 , θ 3DC T 1,2 , . . . , θ 3DC T D,M ]
T where P is an appropriate rearrangement matrix such that ε 3DC = [ϑ 1 , . . . , ϑ D , g T 1 , . . . , g T M , α T 1 , . . . , α T D ] T .

Robust calibration estimator

This section is devoted to the design of a robust calibration estimator based on the model (3.37). As it can be seen from (3.38), one has to specify the pdf of each texture parameter τ pq in order to obtain the exact ML estimates. Nevertheless, in pratical scenarios, such prior knowledge is not available. Consequently, our idea is to make use of a relaxed version of the exact model, i.e., we assume deterministic but unknown texture realizations in the estimation process [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 maximum-likelihood estimates in compound Gaussian noise: Existence and algorithm analysis[END_REF]. This ensures more exibility in our algorithm as the texture distribution is not precisely described and avoids any possible model misspecication, which is consistent with our motivation to design a broad robust estimator w.r.t. the presence of outliers. On the other hand, we adopt here an iterative procedure in order to reduce the computational cost. In doing so, the proposed algorithm sequentially updates each block of unknown parameters while xing the remaining parameters. This leads to the iterative relaxed concentrated ML based calibration In the following, we aim to reduce the computational cost of the minimization procedure in (3.46) by use of the Expectation-Maximization (EM) algorithm. For generality, we rst adopt the nonstructured Jones matrix formulation, which can also be specied depending on the scenario, as shown in Section 3.2.2.2.

Estimation in the case of non-structured Jones matrices

Due to the possible large size of θ, a multi-dimensional parameter search needs to be carried out to solve the optimization problem in (3.46) which requires signicant computation time. To reduce this complexity, we make use of the EM algorithm. As mentioned above, the parameters of interest θ represent the entries of all Jones matrices. Consequently, it is natural to consider the following partition

θ = [θ T 1 , . . . , θ T D ] T = [θ T 1,1 , . . . , θ T 1,M , . . . , θ T D,1 , . . . , θ T D,M ] T , (3.47)
for which the vector θ i,p ∈ R 8×1 is the parametrization of the path from the i-th calibrator source to the p-th sensor, i.e., J i,p (θ) = J i,p (θ i,p ).

Use of the EM algorithm to solve (3.46) The EM algorithm [173175] enables to compute the ML estimates and reduce the computational cost, via the iteration of two steps. The rst one is the E-step which reduces, in our scenario, to the computation of the conditional expectation of the complete data given the observed data and the current estimate of parameters [START_REF] Yatawatta | Radio interferometric calibration using the SAGE algorithm[END_REF][START_REF] Feder | Parameter estimation of superimposed signals using the EM algorithm[END_REF].

Afterwards, the log-likelihood function of this conditional distribution is maximized in the Mstep. As we show in the following, the optimization step is carried out w.r.t. to θ i ∈ C 4M ×1 instead of θ ∈ C 4DM ×1 . Therefore, the global multiple source estimation problem is reduced to multiple single source sub-problems.

1) E-step: For the i-th source, we introduce the so-called complete data vector w i such that

x = D i=1 w i (3.48)

ROBUST CALIBRATION OF RADIO INTERFEROMETERS IN NON-GAUSSIAN ENVIRONMENT 79

with w i = u i (θ i ) + n i and n = D i=1 n i , in which n i ∼ CN (0, β i Ψ). We have D i=1 √ β i = 1 and Ψ is the covariance matrix of n. Since n pq ∼ CN (0, τ pq Ω) and with the independence property, we obtain the following block-diagonal expression for Ψ Ψ = bdiag{τ 12 Ω, . . . , τ (M -1)M Ω}. With [177, p. 36], and after some calculus, the expression of the conditional expectation is given by

ŵi = E{w i |x; θ, τ , Ω} = u i (θ i ) + β i x - D l=1
u l (θ l ) .

(3.51)

2) M-step: The goal of this step is to estimate θ i . Once ŵi are computed for i ∈ {1, . . . , D}, the estimated complete data vector ŵ can be evaluated. The M-step reads as the minimization w.r.t. θ i of the following cost function

φ i (θ i ) = ŵi -u i (θ i ) H (β i Ψ) -1 ŵi -u i (θ i ) . (3.52)
To decrease even more the complexity cost of the proposed robust calibration scheme, we use the Block Coordinate Descent (BCD) algorithm [START_REF] Friedman | Pathwise coordinate optimization[END_REF][START_REF] Hong | A unied algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing[END_REF] in the M-step. Consequently, we obtain analytical solutions for each single source sub-problems in (3.52), as shown below. Notation and calculations being introduced briey in Appendix 6.3, we only present here the results obtained, i.e., the expression of the estimated entries of the Jones matrix associated with the path from the i-th calibrator source to the p-th sensor which is given by θi,p =

   (Σ H i A i,p Σ i + Υ H i Ãi,p Υ i ) -1 (Σ H i A i,p w i,p + Υ H i Ãi,p wi,p ) for 1 < p < M (Σ H i A i,p Σ i ) -1 Σ H i A i,p w i,p for p = 1 (Υ H i Ãi,p Υ i ) -1 Υ H i Ãi,p wi,p for p = M (3.53)
Variables in (3.53) are dened in Appendix 6.3. Therefore, θ i,p for p ∈ {1, . . . , M } are estimated in an iterative loop. With (3.52) and (3.53), it can be proven that the BCD algorithm leads to unique solutions and thus, convergence to at least a local maximizer, is ensured [START_REF] Bertsekas | Nonlinear programming[END_REF]. If the Mstep is performed exactly (i.e., the BCD gives the exact minimizer of (25) and consequently, the M step is exactly solved), convergence of the EM algorithm to a stationary point is ensured (to avoid the unusual case of convergence to a saddle point, a proper initialization is required) [START_REF] Mclachlan | The EM algorithm and extensions[END_REF], with a theoretical innite number of iterations. Finally, in this case, convergence of the concentrated MLE is guaranteed for an innite number of iterations since the value of the cost function at each step can either improve or maintain but cannot worsen [START_REF] Vorobyov | Maximum likelihood direction-of-arrival estimation in unknown noise elds using sparse sensor arrays[END_REF]. In practice, only a nite number of iterations is considered in each loop, so we might not attain local convergence. However, we show in the simulation section the relatively good numerical stability of the algorithm.

The scheme of the proposed algorithm is described in the box Algorithm 3.

Structured Jones matrices

We recall that the output of Algorithm 3 is the estimate of each Jones matrix denoted by Ĵi,p for i ∈ {1, . . . , D} and p ∈ {1, . . . , M }. In the following, we consider the data model in (3.40) for the specic 3DC calibration regime, and intend to estimate the unknown parameter vector RADIO ASTRONOMICAL ARRAYS

p ) = D i=1 || Ĵi,p -G p (g p )H i,p Z i,p F i || 2
F . We rewrite the cost function as

κ(g p ) = D i=1
[X i,p ] k,k = D i=1 [W i,p ĜH p ] k,k = D i=1 [W i,p ] k,k [ Ĝ * p ] k,k (3.57)
for k ∈ {1, 2}, since G p is a diagonal matrix. Therefore, each complex gain element is estimated as

[ĝ p ] k = D i=1 [W * i,p ] k,k -1 D i=1 [X * i,p ] k,k . ( 3
where M i,p = Ĵi,p

F H i H H i,p G H p .
In the case of a compact array, we can write for the i-th source

ϕ T i = αT i Λ (3.62)
where

ϕ i = [φ i,1 , . . . , φi,M ] T and Λ = u 1 , . . . , u M v 1 , . . . , v M
. Therefore, estimation of the directional shifts due to propagation in the ionosphere is given by

αT i = ϕ T i Λ H M p=1 v 2 p -M p=1 u p v p -M p=1 v p u p M p=1 u 2 p M p=1 u 2 p M p=1 v 2 p -( M p=1 u p v p ) 2 . (3.63)
3) Estimation of ϑ i : We consider the following minimization problem θi = arg min

ϑ i M p=1 || Ĵi,p -G p H i,p Z i,p F i (ϑ i )|| 2 F . (3.64)
We assume a large number of antennas M while the number of calibrator sources D is relatively reduced, such that observations outnumber unknown parameters. For each source, the 1D optimization in (3.64) can be computed in a reasonable computational time through a classical data grid search or a Newton type algorithm.

Finally, the proposed algorithm for the structured Jones matrices case regarding 3DC calibration regime is given in the box Algorithm 4.

Extension to the multi-frequency case

In this section, we present the multi-frequency version of our proposed algorithm above in the case of structured Jones matrices. The case of non structured Jones matrices can be deduced easily. For sake of clarity, the index 3DC is omitted in the following Let us recall, that the aim of calibration in a multi-frequency scenario is to estimate the pa-

rameter vector of interest = [ [f 1 ] T , . . . , [f F ] T , g T ] T where [f ] = [ϑ [f ] 1 , . . . , ϑ [f ] D , α [f ] T 1 , . . . , α [f ] T D ] T and g = [g T 1 , . . . , g T M ]
T . To do so, we introduce the following cost function

l [f ] ( [f ] ) = D i=1 l [f ] i ( [f ] i ) (3.65) in which l [f ] i ( [f ] i ) = M p=1 || Ĵ[f] i,p -G p H [f ] i,p Z [f ] ip (α [f ] i )F [f ] i (ϑ [f ] i )|| 2 F with [f ] i = [ϑ [f ] i , α [f ] T i ] T . Prior information on Ĵ[f]
i,p for i ∈ {1, . . . , D}, p ∈ {1, . . . , M } and f ∈ F is provided by the output of the Algorithm 3.

We wish to distributedly solve the following constrained optimization problem thanks to a network of agents

{ˆ [f ] } f ∈F , ẑ = arg min [f 1 ] ,..., [f F ] ,z f ∈F l [f ] ( [f ] ) (3.66) s.t. [f ] i = B [f ] z i , i ∈ {1, . . . , D}, f ∈ F (3.67)
where B [f ] = 1 f 2 I 3 is the known frequency model, z i is an unknown associated global variable, independent w.r.t. frequency and shared by all agents, and z = [z T 1 , . . . , z T D ] T . To solve (3.66), we use a consensus optimization scheme as in the ADMM procedure [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. Instead of considering the original objective function (3.65), we study the following augmented Lagrangian

L( [f 1 ] , . . . , [f F ] , z, y [f 1 ] , . . . , y [f F ] ) = f ∈F D i=1 L [f ] i [f ] i , z i , y [f ] i (3.68)
where

L [f ] i [f ] i , z i , y [f ] i = l [f ] i [f ] i + h [f ] i [f ] i , z i , y [f ] i and h [f ] i [f ] i , z i , y [f ] i = y [f ] T i [f ] i -B [f ] z i + ρ 2 || [f ] i -B [f ] z i || 2 2 .
(3.69)

We note y

[f ] = [y [f ] T 1 , . . . , y [f ] T D ]
T the associated Lagrange parameters (or dual variables) and ρ > 0 a penalty factor. We notice separability of the Lagrangian w.r.t. source direction but above all, separability w.r.t.frequency, meaning that each agent solves a subproblem locally at a given frequency. The ADMM consists in updating sequentially the three following quantities:

• ˆ [f ] i t+1 = arg min [f ] i L [f ] i [f ] i , (ẑ i ) t , ŷ[f] i t (3.70)
performed locally by each agent for i ∈ {1, ..., D}

• (ẑ i ) t+1 = arg min z i f ∈F L [f ] i ˆ [f ] i t+1 , z i , ŷ[f] i t (3.71)
performed globally for i ∈ {1, ..., D}

ROBUST CALIBRATION OF RADIO INTERFEROMETERS IN NON-GAUSSIAN ENVIRONMENT 83

• ŷ[f] i t+1 = ŷ[f] i t + ρ ˆ [f ] i t+1 -B [f ] (ẑ i ) t+1 (3.72)
performed locally by each agent for i ∈ {1, ..., D} where t is the iteration counter. Minimization (3.71) needs access to local solutions from all agents, i.e., at all frequencies, and leads to the following closed-form expression [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF][START_REF] Yatawatta | Distributed radio interferometric calibration[END_REF],

ẑi =   f ∈F ρB [f ] T B [f ]   -1   f ∈F B [f ] T (y [f ] i + ρ [f ] i )   .
(3.73) Minimization (3.70) is addressed iteratively. To this end, we compute the gradient of L

[f ] i ( [f ] i , z i , y [f ] i ) w.r.t. [f ] i , which induces ∂L [f ] i ( [f ] i , z i , y [f ] i ) ∂ϑ [f ] i = [y [f ] i ] 1 + ρ(ϑ [f ] i - 1 f 2 [z i ] 1 ) + M p=1 tr S [f ] i,p + S [f ] H i,p (3.74) where S [f ] i,p = -G p H [f ] i,p Z [f ] i,p ∂F [f ] i (ϑ [f ] i ) ∂ϑ [f ] i Ĵ[f] H i,p , and ∂L [f ] i ( [f ] i , z i , y [f ] i ) ∂η [f ] i = [y [f ] i ] 2 + ρ(η [f ] i - 1 f 2 [z i ] 2 ) + M p=1 tr D [f ] i,p + D [f ] H i,p (3.75) where D [f ] i,p = ju [f ] p Z [f ] * i,p M [f ] i,p and M [f ] i,p = Ĵ[f] i,p F [f ] T i H [f ] H i,p G H p . Likewise, we have ∂L [f ] i ( [f ] i , z i , y [f ] i ) ∂ζ [f ] i = [y [f ] i ] 3 + ρ(ζ [f ] i - 1 f 2 [z i ] 3 ) + M p=1 tr V [f ] i,p + V [f ] H i,p (3.76) where V [f ] i,p = jv [f ] p Z [f ] * i,p M [f ]
i,p . Using (3.74-3.76), we obtain ˆ

[f ] i with a root-nding algorithm or thanks to standard numerical optimization tools as Newton or gradient descent-type algorithm [START_REF] Nocedal | Numerical optimization[END_REF].

Estimation of the gains g p is done as minimization of the following least squares cost function

κ(g p ) = f ∈F D i=1 || Ĵ[f] i,p -G p (g p )H [f ] i,p Z [f ] i,p F [f ] i || 2 F (3.77)
leading to the following estimate of each complex gain element for k ∈ {1, 2}

[ĝ p ] k = f ∈F D i=1 [W [f ] * i,p ] k,k -1 f ∈F D i=1 [X [f ] * i,p ] k,k (3.78)
where

X [f ] i,p = R [f ] i,p Ĵ[f] H i,p and W [f ] i,p = R [f ] i,p R [f ] H i,p in which R [f ] i,p = H [f ] i,p Z [f ] i,p F [f ] i . The
initialize: ˆ ← init , ẑ ← z init , {ŷ [f ] ← y [f ] init } f

Numerical simulations

In this section, we compare the performance of the proposed MSCA with the mono-frequency case, i.e., SCA in which frequency diversity is not taken into account. We recall that radio astronomy observations are aected by the presence of outliers. Thus, we also compare our robust approach with an algorithm based on a classical Gaussian noise assumption [START_REF] Yatawatta | Radio interferometric calibration using the SAGE algorithm[END_REF], which amounts to solve a non-linear least squares problem. First, in Fig. 3.10, we plot the mean square error of η

[f 1 ] 1
as a function of the signal-to-noise ratio, the behavior being the same for any other parameter of . We compare the estimation performance for dierent number of frequencies F and notice better statistical performance when multi-frequency robust calibration is performed.

Robust calibration based on the Student's t [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF] is not exposed in the simulations due to a dierent model which is not adapted to the 3DC regime.

In the next gures, we use Meqtrees [START_REF] Noordam | The MeqTrees software system and its use for thirdgeneration calibration of radio interferometers[END_REF] to generate the data model, the observations and compare its least squares solver to MSCA. Here, we choose to correct for Faraday rotation matrices, which are the only introduced perturbations in the observations. We consider M = 7 antennas (KAT-7 instrument), D = 1, D = 16 weak realistic background sources taken from the SUMSS survey using a spectral index of 0.7. The full duration of the observation is 12 hours, for 60 seconds integration time per data sample. After calibration and subtraction of the bright calibrator source, a dirty image, namely the corrected residual, is constructed with Meqtrees using lwimager. Fig. 3.11 gives the corrected residual image at 895 MHz in a small area surrounding the calibrator, whose position corresponds to the red cross. Fig. 3.12 gives the recovered ux for one of the D weak outlier sources. Therefore, we notice better ux estimation of weak background sources and better calibrator removal using joint frequency dependent calibration with MSCA compared to a frequency independent calibration.

We consider an other case in which M = 7 receptors, D = 1 and D = 16 taken from the SUMSS survey with a spectral index of 0.7. An image area of 3.5 by 3.5 degrees is shown in The results presented here are technical, the reader can refer to its application in the radar context in our paper [START_REF] Ren | Performance bounds under misspecication model for MIMO radar application[END_REF].

• The second non standard condition presented in this Chapter is the presence of discrete and continuous parameter in the unknown parameter vector. It is worth mentioning that among the plethora of lower bounds on the MSE, the CRB is certainly the most famous one in the signal processing community. Indeed, the attractiveness of the CRB comes from the fact that various closed-form expressions are available for a large class of observation models. However, an important point in our context, is that the unknown vector parameter contains discrete parameters, consequently, the regularity conditions to apply the CRB are not fullled since the likelihood involved in the Fisher information computation is not dierentiable w.r.t. these parameters. Our alternative is to propose a Hybrid Cramér-Rao-Weiss-Weinstein bound (HCRWWB) and derive its associated closed-form expressions whatever the considered distribution of the data. Furthermore, we consider the change point application for which contrary to existing works about performance analysis in the change-point literature, our study is adapted to multiple changes but also the unknown observations distribution's parameters. Finally, a strategy in order to pick the The asymptotic performance analysis of an estimation algorithm mostly relies on two simplied assumptions: i) the data are assumed to be Gaussian distributed and ii) the data model used to derive the estimation algorithm is supposed to be correctly specied, that is the pdf assumed to derive an estimator of the parameters of interest and the true pdf that statistically characterizes the data are exactly the same.

Although these assumptions guarantee the possibility to perform elegant performance assessment, e.g., by evaluating the CRB for the estimation problem at hand and/or by obtaining a closed form expression for the Mean Square Error (MSE) of a given estimator, the everyday engineering practice clearly calls the hypotheses i) and ii) into question. Regarding the Gaussian model assumption, large-scale measurement campaigns and the subsequent statistical analysis of the data gathered from a plethora of engineering applications have highlighted the impulsive, heavy-tailed behaviour of the observations [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. These experimental evidences have motivated the need to go beyond the Gaussian model and develop new statistical models able to better characterize the data. One of the more exible and general non-Gaussian model is represented by the set of the CES distributions [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]. Let us recall that the pdf of a CES distributed Ndimensional random vector x l ∈ C N is completely characterized by the mean value γ, the scatter (or shape) matrix Π and by a real valued function w(t) : R + → R, called the density generator, i.e., x l ∼ CES N (γ, Π, w) [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF][START_REF] Richmond | Adaptive array signal processing and performance analysis in non-Gaussian environments[END_REF].

Other experimental evidences reveal recurring violations of the matched model assumption, that is the claim of a perfect match between the assumed and the true data model. The mathematical bases of a formal theory of the parameter estimation under model misspecication has been rstly developed by statistician as Huber [START_REF] Huber | The behavior of maximum likelihood estimates under nonstandard conditions[END_REF], White [START_REF] White | Maximum likelihood estimation of misspecied models[END_REF] and Vuong [START_REF] Vuong | Cramér-rao bounds for misspecied models[END_REF] and recently rediscovered by the Signal Processing (SP) community [193195] and applied to a variety of well-known engineering problems [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF][START_REF] Ren | Performance bounds under misspecication model for mimo radar application[END_REF].

This brief discussion clearly highlights the need to overtake both the Gaussian and the matched model assumptions while assessing the (asymptotic) performance of an estimator. Under the matched model and Gaussian assumptions, the CRB can be evaluated thanks to the Slepian-Bangs (SB) formula. The rst generalization of the SB formula to a non-Gaussian, but still perfectly matched, data model has been proposed by Besson and Abramovich in [START_REF] Besson | On the sher information matrix for multivariate elliptically contoured distributions[END_REF]. The second important step ahead has been made by Richmond and Horowitz in [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF] and then by Parker and Richmond in [START_REF] Parker | Methods and bounds for waveform parameter estimation with a misspecied model[END_REF] where the classical, Gaussian-based, SB formula has been extended to estimation problems under model misspecication, i.e., when the assumed parametric Gaussian model, say CN (γ(θ), Π(θ)), could dier from the true (possibly non parametric) one, indicated as CN (µ, Σ). In other words, we allow the assumed parametric mean value γ(θ) and the assumed parametric covariance matrix Π(θ) to dier from the true µ and Σ for every possible value of the parameter vector θ ∈ Θ, i.e., CN (γ(θ), Π(θ)) = CN (µ, Σ), ∀θ ∈ Θ. It is worth to underline that in the estimation framework under model misspecication, the FIM loses its classical statistical sense and it has to be substituted by the matrices A(θ) and B(θ) dened in [START_REF] Fortunati | The Misspecied Cramér-Rao bound and its application to scatter matrix estimation in complex elliptically symmetric distributions[END_REF], eqs. ( 1) and ( 7), respectively. Consequently, in this context, SB-type formulas could be exploited to obtain A(θ) and B(θ) needed to evaluate the counterpart of the CRB in the presence of model misspecication, i.e., the Misspecied CRB (MCRB) [START_REF] Huber | The behavior of maximum likelihood estimates under nonstandard conditions[END_REF]192194,[START_REF] Greco | Maximum likelihood covariance matrix estimation for complex elliptically symmetric distributions under mismatched conditions[END_REF]. In particular, in [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF] the authors derived SB-type formulas for the decoupled scenario in which the unknown parameter vector θ ∈ Θ can be partitioned in two sub-vectors η and ν, named deterministic and stochastic parameter sub-vectors respectively, such that θ = [η T , ν T ] T and CN (γ(θ), Π(θ)) CN (γ(η), Π(ν)) = CN (µ, Σ), ∀θ ∈ Θ. The ndings presented in [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF] have 90 CHAPTER 4. PERFORMANCE ANALYSIS WITH APPLICATIONS TO ARRAY PROCESSING been extended in [START_REF] Parker | Methods and bounds for waveform parameter estimation with a misspecied model[END_REF] to include the coupling of deterministic and stochastic parameters. More formally, in [START_REF] Parker | Methods and bounds for waveform parameter estimation with a misspecied model[END_REF], SB-type formulas have been derived for the following misspecied scenario

CN (γ(θ), Π(θ)) CN (γ(η, ω), Π(ν, ω)) = CN (µ, Σ), ∀θ ∈ Θ where the unknown parameter vector θ ∈ Θ is partitioned as θ = [η T , ν T , ω T ] T .
The natural extension of the works [START_REF] Besson | On the sher information matrix for multivariate elliptically contoured distributions[END_REF] and [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF][START_REF] Parker | Methods and bounds for waveform parameter estimation with a misspecied model[END_REF] would be to derive SB-type formulas for parametric estimation problems involving CES distributed data under model misspecication.

In this section we aim exactly at lling this gap and obtaining general misspecied SB formulas for CES distributed data.

Remark : Throughout this section, we consider only the case of real parameter vectors. This is not a limitation, since we can always maps a complex vector in a real one simply by stacking its real and the imaginary parts. Clearly, the proposed derivation of the SB-type formulas could also be developed directly in the complex eld by means of the Wirtinger calculus as in [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF].

Notation: Throughout this section, let f (t) be a real scalar function, than f (t) df (t)/dt. Let A(θ) be a matrix (or possibly vector or even scalar) function of the vector θ, then A 0 A(θ 0 ) while A 0 i ∂A(θ)

∂θ i | θ=θ 0 and A 0 ij ∂ 2 A(θ)
∂θ i ∂θ j | θ=θ 0 , where the vector θ 0 will be always explicitly dened.

Problem setup

Let x = {x l } L l=1 , with x l ∈ C N , be a set of L independent complex random vectors (or snapshots)

representing the available observations. We assume that each snapshot is sampled from a CES distribution [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF][START_REF] Richmond | Adaptive array signal processing and performance analysis in non-Gaussian environments[END_REF], i.e., x l ∼ CES N (µ l , Σ, g), then its pdf can be expressed as:

p X (x l ) p X (x l ; µ l , Σ) = c N,g |Σ| -1 g((x l -µ l ) H Σ -1 (x l -µ l )) (4.1)
where c N,g is a normalizing constant, g(t) : R + → R is the density generator, µ l = E p {x l } is the mean value and Σ is a positive denite Hermitian matrix called scatter matrix. In the following, we always assume that the scatter matrix Σ is of full rank, i.e., rank (Σ) = N . From the Stochastic Representation Theorem [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF], a CES distributed random vector can be expressed as:

x l = d µ l + RTu l , (4.2)
where:

• u l ∼ U (CS N ) is a N -dimensional complex random vector uniformly distributed on the unit hyper-sphere with N -1 topological dimension. As reported in [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF] (Lemma 1), E p {u l } = 0 and E p u l u H l = (1/N )I where I is the identity matrix of a suitable dimension.

• R √ Q is a real and non-negative random variable called modular variate, while Q is called second order modular variate. Moreover, under the assumption that rank(Σ) = N , we have that:

Q l (x l -µ l ) H Σ -1 (x l -µ l ) d = Q, ∀l ∈ N. (4.3)
As shown in [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF], the pdf of Q has a one-to-one relation with density generator:

p Q (t) = δ -1 N,g t N -1 g(t), (4.4)
where δ N,g ∞ 0 t N -1 g(t)dt < ∞. As a consequence of (4.3) and (4.4), the expectation of functions of the quadratic form Q l , say h(Q l ), can be explicitly derived as:

E Q {h(Q l )} ∞ 0 h(t)p Q (t)dt = δ -1 N,g ∞ 0 h(t)t N -1 g(t)dt, ∀l ∈ N.

MISSPECIFIED CRAMÉR-RAO BOUND AND ITS SLEPIAN-BANGS-TYPE FORMULAS FOR CES DISTRIBUTIONS 91

It is clear from (4.5) that such expectation does not depend on the index l, since the pdf in (4.4) of the quadratic form Q l is invariant w.r.t. l. For this reason, to avoid confusion, in the rest of the section we always indicate

E Q {h(Q l )} simply as E Q {h(Q)}.
• T is a complex N × N matrix with rank (T) = N , such that Σ = TT H .

• If E Q {Q} < ∞ and rank(Σ) = N , then the covariance matrix M = E p {(x l -µ l )(x l -µ l ) H } of x l can be decomposed as M = σ 2 Σ, where (see Theorem 4 in [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]):

σ 2 E Q {Q} N , (4.6)
and σ 2 can be interpreted as the statistical power of the CES distributed vector x l .

From the Stochastic Representation Theorem, it is clear that the representation of a CES distributed vector x l is not uniquely determined by (4.2). In fact,

x l d = µ l + RTu l d = µ l + (c -1 R)(cT)u l , ∀c > 0.
From a dierent, yet equivalent, standpoint, this identiability problem can be understood as an implicit consequence of the functional expression of a CES distribution.

It is immediate to verify from (4.1) that CES N (µ l , Σ, g(t)) ≡ CES N (µ l , c 2 Σ, g(t/c 2 )), ∀c > 0.

As amply discussed in [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF], this identiability issue can be solved by posing a constraint on the modular variate R (and consequently, through (4.4), on the density generator g(t)), or on the scatter matrix Σ. For convenience, we choose to put a constraint on R 2 = Q. Specically in the rest of this section, we always assume that E Q {Q} = N and consequently, from (4.6), M = Σ, i.e., the scatter matrix equates the covariance matrix.

Due to the independence assumption, the joint pdf of the set x is given by the product of the marginal pdfs of each snapshot x l given in (4.1):

p X (x) p X (x 1 , . . . , x L ; µ 1 , . . . µ L , Σ) = (c N,g ) L |Σ| -L L l=1 g((x l -µ l ) H Σ -1 (x l -µ l )) (4.7)
In th following we consider the following mismatched situation:

• The acquired dataset x = {x l } L l=1 is characterized by the true but unknown joint pdf given in (4.7). In particular, each data snapshot x l follows a CES distribution with mean value µ l , scatter matrix Σ and density generator g(t), i.e., x l ∼ CES N (µ l , Σ, g), l = 1, . . . , L.

• In order to derive an inference algorithm, we assume that each snapshot of the dataset

x is sampled from a CES distribution with a density generator w(t), possibly dierent from g(t) for all t ∈ R + , and a mean value γ l (θ) and a scatter matrix Π(θ) parametrized by a deterministic parameter vector θ ∈ Θ ⊂ R d to be estimated. In particular, we allow the assumed marginal model CES N (γ l (θ), Π(θ), w) to dier from the true one, CES N (µ l , Σ, g) for every θ ∈ Θ.

This is a recurring scenario in array processing applications, where the mean value and/or the scatter matrix of the acquired snapshot vectors are assumed to be parametrized by a deterministic parameter vector whose components represent the Doppler frequency, the DoAs of potential sources and so on.

The assumed marginal pdf of each snapshot x l can then be expressed as:

f X (x l ; θ) f X (x l ; γ l (θ), Π(θ)) = c N,w |Π(θ)| -1 w((x l -γ l (θ)) H Π(θ) -1 (x l -γ l (θ))) (4.
x = {x l } L l=1 reads f X (x; θ) f X (x 1 , . . . , x L ; γ l (θ), Π(θ)) = (c N,w ) L |Π(θ)| -L L l=1 w((x l -γ l (θ)) H Π(θ) -1 (x l -γ l (θ))) (4.9)
This scenario clearly represents an estimation problem in non-Gaussian data and in the presence of model misspecication. Let θf θf (x) be a, possibly mismatched, estimator of the parameter vector θ ∈ Θ derived under the assumed model f X (x; θ) in (4.9) while the data are characterized by the true model p X (x) in (4.7). Then, as discussed in [192194], under suitable regularity conditions, a lower bound on the error covariance matrix of any misspecied (MS)unbiased (see [START_REF] Vuong | Cramér-rao bounds for misspecied models[END_REF][START_REF] Fortunati | The Misspecied Cramér-Rao bound and its application to scatter matrix estimation in complex elliptically symmetric distributions[END_REF]) mismatched estimator θf exists and it is given by the MCRB dened as:

C p θf , θ 0 E p ( θf -θ 0 )( θf -θ 0 ) T ≥ 1 L A -1 (θ 0 )B(θ 0 )A(θ 0 ) -1 (4.10)
where θ 0 is the so-called pseudo-true parameter vector dened in [START_REF] Huber | The behavior of maximum likelihood estimates under nonstandard conditions[END_REF], [START_REF] White | Maximum likelihood estimation of misspecied models[END_REF] and [START_REF] Vuong | Cramér-rao bounds for misspecied models[END_REF] as the parameter vector that minimize the Kullback-Leibler distance (KLD) between the true and the assumed models:

θ 0 arg min θ∈Θ {D (p X f θ )} = arg min θ∈Θ {-E p {ln f X (x l ; θ)}} , (4.11)
where D(p X f θ ) E p {ln(p X (x l )/f X (x l ; θ))}. The matrices A(θ 0 ) and B(θ 0 ) are dened as:

[A(θ 0 )] i,j E p ∇ T θ ∇ θ ln f X (x l ; θ 0 ) i,j = E p ∂ 2 ln f X (x l ; θ) ∂θ i ∂θ j θ=θ 0 (4.12)
and

[B(θ 0 )] i,j E p ∇ θ ln f X (x l ; θ 0 )∇ T θ ln f X (x l ; θ 0 ) i,j = E p ∂ ln f X (x l ; θ) ∂θ i θ=θ 0 • ∂ ln f X (x l ; θ) ∂θ j θ=θ 0 . (4.13)
Our goal is to provide a general closed-form expressions of A(θ 0 ) and B(θ 0 ) in the aforementioned context. In other words, we aim at deriving two SB-type formulas for the evaluation of the matrices A(θ 0 ) and B(θ 0 ) that represent a generalization of the classical FIM for estimation problem under model misspecication. It is important to note that in the following derivations we always assume the existence and the uniqueness of the pseudo-true parameter vector θ 0 dened in (4.11). As it is clear from (4.10) and as it is amply discussed in [192194], nding the θ 0 that minimizes the KLD between the true and the assumed model is a prerequisite for the evaluation of the MCRB since all the derivatives involved in the two matrices A(θ 0 ) and B(θ 0 ) have to be evaluated at θ 0 . Finally, it is worth noticing that, under suitable regularity conditions on the true and the assumed pdfs, the denition of the pseudo-true parameter vector θ 0 in (4.11) can be expressed in an equivalent form as: This Section focuses on the derivation of the SB formulas for A(θ 0 ) and B(θ 0 ) dened in (4.12) and (4.13), provided that there exists a unique pseudo-true parameter vector θ 0 satisfying (4.11).

∂D (p X f θ ) ∂θ i θ=θ 0 = -E p ∂ln f X (x; θ) ∂θ i θ=θ 0 = 0, i = 1, . . . , d.
We start by evaluating explicit expressions for the following quantities:

V ij (θ 0 ) = ∂ln f X (x; θ) ∂θ i θ=θ 0 • ∂ln f X (x; θ) ∂θ j θ=θ 0 , i, j = 1, . . . , d, (4.15)
and

H ij (θ 0 ) = ∂ 2 ln f X (x; θ) ∂θ i ∂θ j θ=θ 0
, i, j = 1, . . . , d.

(4.16)

The entries of the matrices A(θ 0 ) and B(θ 0 ) can then be obtained by taking the expectation operator, w.r.t. the true distribution p

X (x), of H ij (θ 0 ) and V ij (θ 0 ), i.e., [A(θ 0 )] ij = E p {H ij (θ 0 )} and [B(θ 0 )] ij = E p {V ij (θ 0 )}.
Unfortunately, this expectation can be evaluated in closed form only in two particular cases, as we will detail below. In all the other cases, numerical integration techniques, e.g., the Monte Carlo integration, could be exploited.

Evaluation of (4.15), (4.16) and of their expectation w.r.t. the true data distribution According to the general mismatched estimation problem discussed in the previous section, we consider here the general case in which, for each available snapshot x l , a parametric CES model f X (x l ; θ) = CES N (γ l (θ), Π(θ), w) is assumed, while actually each observation vector is distributed according to a dierent, possibly non-parametric CES data model, i.e.,

x l ∼ p X (x l ) = CES N (µ l , Σ, g).

From (4.9), it is immediate to verify that:

∂ ln f X (x; θ) ∂θ i θ=θ 0 = -L ∂ ln |Π(θ)| ∂θ i θ=θ 0 + L l=1 φ(G l (θ 0 )) ∂G l (θ) ∂θ i θ=θ 0 = -Ltr(P 0 i ) + L l=1 φ(G 0 l ) ∂G 0 l ∂θ i (4.17)
where

∂ ln |Π| ∂θ i θ=θ 0 = tr Π -1 0 Π 0 i = tr(P 0 i ), (4.18)
with φ(t) = w (t)/w(t) and

P 0 i Π -1/2 0 Π 0 i Π -1/2 0 , (4.19) G 0 l G l (θ 0 ) (x l -γ 0 l ) H Π -1 0 (x l -γ 0 l ), (4.20)
where, for notational simplicity, γ 0 l γ l (θ 0 ) and Π 0 Π(θ 0 ). Then, an explicit expression of V ij (θ 0 ) in (4.15) is given by:

V ij (θ 0 ) =L 2 tr(P 0 i )tr(P 0 j ) -Ltr(P 0 i ) L l=1 φ(G 0 l ) ∂G 0 l ∂θ j -Ltr(P 0 j ) L l=1 φ(G 0 l ) ∂G 0 l ∂θ i + + L l=1 L m=1 φ(G 0 l )φ(G 0 m ) ∂G 0 l ∂θ i ∂G 0 m ∂θ j , (4.21)
where

∂G 0 l ∂θ i = -2Re (x l -γ 0 l ) H Π -1 0 ∂γ 0 l ∂θ i -(x l -γ 0 l ) H S 0 i (x l -γ 0 l ), (4.22
S 0 i = Π -1 0 Π 0 i Π -1 0 . (4.23)
The term H ij (θ 0 ) in (4.16) can be obtained, through direct calculation, from (4.17) as:

H ij (θ 0 ) = Ltr(P 0 i P 0 j -P 0 ij ) + L l=1 φ (G 0 l ) ∂G 0 l ∂θ j ∂G 0 l ∂θ i + L l=1 φ(G 0 l ) ∂ 2 G 0 l ∂θ i ∂θ j (4.24)
where

P 0 ij Π -1/2 0 Π 0 ij Π -1/2 0 (4.25)
and φ (t) = w (t)/w(t) -(w (t)) 2 /w 2 (t). After having obtained the terms V ij (θ 0 ) and H ij (θ 0 ), we have to evaluate their expectations w.r.t. the true distribution p X (x). Since all the derivatives in (4.21) and (4.24) have to be evaluated in the pseudo-true parameter vector θ 0 , we can exploit the equality established in (4.14). Specically, from (4.17), we have that:

∂D (p X f θ ) ∂θ i θ=θ 0 = Ltr(P 0 i ) - L l=1 E p φ(G 0 l ) ∂G 0 l ∂θ i = 0, i = 1, . . . , d, (4.26)
and consequently,

L l=1 E p φ(G 0 l ) ∂G 0 l ∂θ i = Ltr(P 0 i ), i = 1, . . . , d. (4.27)
Now, taking the expectation operator w.r.t. p X (x) of the term V ij (θ 0 ) in (4.21) and by exploiting the equality in (4.27), the matrix B(θ 0 ) can be expressed as:

[B(θ 0 )] ij = -L 2 tr(P 0 i )tr(P 0 j ) + L l=1 L m=1 E p φ(G 0 l )φ(G 0 m ) ∂G 0 l ∂θ i ∂G 0 m ∂θ j . (4.28)
Similarly, the matrix A(θ 0 ) can be obtained by taking the expectation operator w.r.t. p X (x) of the term H ij (θ 0 ) in (4.24) as:

[A(θ 0 )] ij = Ltr(P 0 i P 0 j -P 0 ij ) + L l=1 E p φ (G 0 l ) ∂G 0 l ∂θ j ∂G 0 l ∂θ i + L l=1 E p φ(G 0 l ) ∂ 2 G 0 l ∂θ i ∂θ j . (4.29)
As we can see from (4.21) and (4.24), the expressions of V ij (θ 0 ) and H ij (θ 0 ) are highly involved from an analytical standpoint, and consequently it is impossible to derive in closed form their expectations in the general case. As discussed in Appendix 6.4, closed form expressions can be obtained when the random terms in (4.21) and (4.24) satisfy certain independence conditions.

Fortunately, there are two scenarios of great practical interest in which such conditions are met and consequently a closed form expression for A(θ 0 ) and B(θ 0 ) can be obtained. These special cases are Scenario 1: The true marginal pdf is an unspecied CES model, i.e., p X (x l ) = CES N (µ l , Σ, g), while the assumed pdf is a parametric complex Gaussian model, i.e., f X (x l ; θ) = CN (γ l (θ), Π(θ)).

Scenario 2: The Scenario 2 is characterized by two assumptions: A1 : The true and the assumed pdfs share the same parametric mean value γ l (θ) and the same parametric scatter matrix Π(θ) while the misspecication is caused by a wrong assumption on the density generator w(t). More formally, posit that the true marginal pdf of each snapshot is given by a parametric CES distribution such that p X (x l ) p X (x l ; θ) = CES N (µ l , Σ, g), where µ l = γ l ( θ) and Σ = Π( θ) for a given θ ∈ Θ. The assumed model is instead another parametric CES distribution that share the same parameterization of the true one but may have a dierent density generator, i.e., f X (x l ; θ) = CES N (γ l (θ), Π(θ), w), θ ∈ Θ, and g(t) = w(t), ∀t ∈ R + .

A2 : The true parameter vector θ and the pseudo-true parameter vector θ 0 are equals. In particular, θ is the solution of the optimization problem in (4.11).

Note that these two assumptions are veried for the scatter matrix estimation problem discussed in [START_REF] Fortunati | The Misspecied Cramér-Rao bound and its application to scatter matrix estimation in complex elliptically symmetric distributions[END_REF].

The Scenario 1 describes a common practice in array processing applications. In fact, when the true data model is unknown, a prevalent choice is to assume a simple Gaussain model that guarantees an easy derivation and a consequent real-time implementation of the estimation algorithm.

The Scenario 2 is a bit dierent, since it imply the a-priori knowledge of the functional form of the true parametric mean value γ l (θ) and of the parametric scatter matrix Π(θ). There are, however, a variety of practical applications in which this a-priori information is indeed available to the user. As an example, one can think of array signal processing applications in which the apriori knowledge of the array geometry leads to a correct specication of the parametrized mean value and covariance matrix of the collected snapshots, while the uncertainty on the statistical disturbance model could cause a wrong choice of the density generator. Of course, there are cases in which also the knowledge of the array geometry could be incorrect or partial, and consequently the assumed mean value and scatter matrix dier from the true ones (see [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF] for more details).

Scenario 1 In this subsection we provide the SB formulas, i.e., the closed form expressions of the matrices A(θ 0 ) in (4.12) and B(θ 0 ) in (4.13), for the Scenario 1. We start by noticing that the assumed complex Gaussian model belongs to the CES class, i.e., f X (x l ; θ) = CES N (γ l (θ), Π(θ), w), where w(t) = exp(-t). Consequently, it is immediate to verify that φ(t) = -1 and φ (t) = 0. The matrix B(θ 0 ) can be evaluated from (4.28) by using the fact that φ(G l ) = -1:

[B(θ 0 )] i,j = -L 2 tr(P 0 i )tr(P 0 j ) + L l=1 L m=1 E p ∂G 0 l ∂θ i ∂G 0 m ∂θ j . (4.30)
Following the procedure discussed in Appendix 6.4.2, the matrix B(θ 0 ) can be expressed as:

[B(θ 0 )] i,j = 2 L l=1 Re Π 0 i Π -1 0 r 0 l + ∂γ 0 l ∂θ i H Π -1 0 ΣΠ -1 0 Π 0 j Π -1 0 r 0 l + ∂γ 0 l ∂θ j + + L E Q {Q 2 } N (N + 1) -1 tr(S 0 i Σ)tr(S 0 j Σ) + L E Q {Q 2 } N (N + 1) tr(S 0 i ΣS 0 j Σ), (4.31)
where

r 0 l µ l -γ 0 l . (4.32)
Note that all the derivatives have to be evaluated in the pseudo-true parameter vector θ 0 dened in (4.11).

CHAPTER 4. PERFORMANCE ANALYSIS WITH APPLICATIONS TO ARRAY PROCESSING

Let us evaluate the matrix A(θ 0 ). From (4.24) and φ(G l ) = -1, and φ (G l ) = 0, through direct calculation (see Appendix 6.4.2) we obtain:

α ij l (θ 0 ) E p ∂ 2 G 0 l ∂θ i ∂θ j = 2Re ∂γ 0 l ∂θ i H Π -1 0 ∂γ 0 l ∂θ j + + 2Re r 0 l H S 0 j ∂γ 0 l ∂θ i + S 0 i ∂γ 0 l ∂θ j -Π -1 0 ∂ 2 γ 0 l ∂θ i ∂θ j + + r 0 l H Π -1/2 0 (P 0 i P 0 j + P 0 j P 0 i -P 0 ij )Π -1/2 0 r 0 l + + tr (P 0 i P 0 j + P 0 j P 0 i -P 0 ij )Π -1/2 0 ΣΠ -1/2 0 . (4.33)
Moreover, as in (4.31), we used constraint on Q l . Finally, by inserting (4.33) in (4.24), we obtain:

[A(θ 0 )] i,j = Ltr(P 0 i P 0 j -P 0 ij ) - L l=1 α ij l (θ 0 ). (4.34)
The expressions (4.31) and (4.34) represent the SB formulas for the mismatched Scenario 1.

Scenario 2 This Subsection focuses on the Scenario 2, i.e., the case in which the true and the assumed pdfs are CES distributions that share the same parametrized mean value and scatter matrix but are characterized by dierent density generators. From the proof provided in Appendix 6.4.3, the matrices B(θ 0 ) and A(θ 0 ) can be expressed respectively as:

[B(θ 0 )] i,j = [B( θ)] i,j = 2 N E Q {Qφ 2 (Q)} L l=1 Re ∂γ 0 l ∂θ i H Π -1 0 ∂γ 0 l ∂θ j + L E Q {Q 2 φ 2 (Q)} N (N + 1) -1 tr(P i )tr(P j ) + L E Q {Q 2 φ 2 (Q)} N (N + 1) tr(P i P j ), (4.35)
[A(θ 0 )] i,j = [A( θ)] i,j = 2 N E Q {Qφ (Q)} + N E Q {φ(Q)} L l=1 Re ∂γ 0 l ∂θ i H Π -1 0 ∂γ 0 l ∂θ j + L E Q {Q 2 φ (Q)} N (N + 1) -1 tr(P 0 i P 0 j ) + L E Q {Q 2 φ (Q)} N (N + 1) tr(P 0 i )tr(P 0 j ). (4.36)
The expressions (4.35) and (4.36) represent the SB formulas for the mismatched Scenario 2. It is worth to recall that the previous two formulas can be applied only if, in the particular estimation problem an hand, the pseudo-true parameter vector equates the true parameter vector, i.e., when θ 0 = θ. The reason for this restriction is claried in Appendix 6.4.3.

Relationship to previous results

It is possible to show that the SB formulas derived for the Scenarios 1 and 2 encompass all the previously derived SB formulas as special cases. In particular, both the SB formula for CES distributions under perfect model specication, proposed in [START_REF] Besson | On the sher information matrix for multivariate elliptically contoured distributions[END_REF], and the SB formulas for the scatter/covariance matrix estimation of CES distributed vectors under misspecication of the density generator, proposed in [START_REF] Fortunati | The Misspecied Cramér-Rao bound and its application to scatter matrix estimation in complex elliptically symmetric distributions[END_REF][START_REF] Greco | Maximum likelihood covariance matrix estimation for complex elliptically symmetric distributions under mismatched conditions[END_REF], can be obtained as special cases of the SB 4.1. MISSPECIFIED CRAMÉR-RAO BOUND AND ITS SLEPIAN-BANGS-TYPE FORMULAS FOR CES DISTRIBUTIONS 97

formulas shown here for the Scenario 2. In addition, the SB formulas for misspecied Gaussian models, proposed in [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF], are special cases of the SB formulas discussed here for the Scenario 1. We note, in passing, that the SB-type formulas derived in [START_REF] Parker | Methods and bounds for waveform parameter estimation with a misspecied model[END_REF] can also be obtained as a particular case of the ones proposed in (4.31) and (4.34). This can be easily done by partitioning the unknown parameter vector as θ = [η T , ν T , ω T ] T and by taking into account the particular parameterization of the assumed mean value and of the assumed scatter matrix supposed in [START_REF] Parker | Methods and bounds for waveform parameter estimation with a misspecied model[END_REF],

i.e., γ(θ) = γ(η, ω) and Π(θ) = Π(ν, ω). For sake of clarity, these discussions and proofs are reported in Appendices 6.4.4-6.4.1. that all the other parameters were assumed to be known (e.g., in [START_REF] Bacharach | Weiss-weinstein bound on multiple change-points estimation[END_REF][START_REF] Bacharach | Prior inuence onweiss-weinstein bounds for multiple change-point estimation[END_REF], the means and variances of the Gaussian distributions associated with the dierent segments were assumed to be known). Consequently, in this section, we ll this lack. More precisely, we propose a hybrid bound for multiple change-point estimation when the discrete change locations and the continuous parameters of the distributions associated with the dierent segments are both unknown. To achieve this, we propose the Hybrid Cramér-RaoWeiss-Weinstein bound (HCRWWB), in which the CRB part is associated with the continuous parameters, and the WWB part is associated with the discrete parameters.

PERFORMANCE ANALYSIS

For sake of clarity, in the following we focus on the change-points application, nevertheless, the same methodology can be adapted to any parametric estimation problem in the presence of unknown continuous and discrete parameters.

Problem statement

Let us now introduce the observation model used throughout this section. We consider a time series x = [x 1 , . . . , x N ] T with independent random variables x n belonging to the observation space denoted by Ω ⊂ R N . This time series is submitted to multiple abrupt changes, that arise at unknown time instants t 1 , . . . , t Q , also referred to as change-points or changes. The total number of changes Q is assumed to be known. Thus, the observation model can be written as

           x n ∼ f (x n ; η 0 ) for n = 1, . . . , t 1 x n ∼ f (x n ; η 1 ) for n = t 1 + 1, . . . , t 2 . . . x n ∼ f (x n ; η Q ) for n = t Q + 1, . . . , N (4.37)
in which f (x n ; η q ) denotes the distribution of the random variables x n in the (q + 1)-th segment, namely the segment delineated by the two consecutive change-points t q and t q+1 , with q ∈ {0, . . . , Q} (that is the set of successive integers between 0 and Q), t 0 0 and t Q+1 N . These distributions f (. ; η q ) are parameterized by parameter vectors η q = [η q 1 , . . . , η q L ] T ∈ R L (for instance, in the case of Gaussian distributions, L = 2 and the parameter vector η q includes the mean and and the variance of this Gaussian distribution for the (q + 1)-th segment). We assume that all the distributions f (. ; η q ), for q ∈ {0, . . . , Q}, belong to the same family.

The parameter estimation problem in such scenario consists in estimating i) the change-point locations t q , q = 1, . . . , Q, and ii) the parameter vectors η q , q = 0, . . . , Q. Thus, the unknown parameter vector to estimate is θ

= [η T , t T ] T ∈ Θ ⊂ R L(Q+1)+Q , with η = [η T 0 , . . . , η T Q ] T ∈ Θ η ⊂ R L(Q+1) and t = [t 1 , . . . , t Q ] T ∈ Θ t ⊂ Z Q .
The sets Θ, Θ η and Θ t denote the parameter spaces for θ, η and t, respectively, that is, Θ = Θ η × Θ t .

The purpose of this section is to assess estimation performance for the parameter vector θ by providing lower bounds on the MSE for a family of estimators θ(x) of θ. Specically, we adopt a hybrid point of view, in the sense that we consider the parameter vectors η q stacked in η as being unknown and deterministic, with true values η q (accordingly, the true value of the full parameter vector is denoted by η ), and the parameter vector t as being random. Consequently, the estimator θ(x) will be hybrid as well, for example it can be the ML-MAP estimator (Maximum Likelihood-Maximum A Posteriori) [START_REF] Yeredor | The joint MAP-ML criterion and its relation to ML and to extended leastsquares[END_REF]. The context of hybrid estimation is appropriate in cases where no a priori information is available on some of the unknown parameters. In addition, interestingly, the hybrid point of view makes it possible to get rid of some integrals, and thus, to obtain genuine closed-form expressions for the bound.

Since parameter vector t is random, in agreement with the Bayesian framework, it is assigned a prior distribution denoted by π(t). A convenient and commonly used prior is the uniform random walk. In other words, t q is assumed to be given by t q = t q-1 + q , q = 1, . . . , Q, where t 0 0 and q are i.i.d. uniformly distributed random variables on the set of integers {1, . . . , τ }. The value of τ can be freely chosen as long as the last change t Q occurs at least before the end of the observations, i.e., the maximum possible value for τ is τ max = (N -1)/Q , with . denoting the oor function. Consequently, the joint (discrete) prior distribution for parameter vector t is given by

π(t) = 1 τ Q Q q=1 I {t q-1 +1,...,t q-1 +τ } (t q ). (4.38)
The support of this prior law is denoted by T t = [t 1 , . . . , t Q ] T ∈ Z Q ∀q ∈ {1, . . . , Q}, t q ∈ {t q-1 + 1, . . . , t q-1 + τ }, t 0 = 0, t Q < N , and it corresponds to the set of the possible segmentations of the observation window {1, . . . , N } in exactly Q segments, each with maximum length τ .

From model (4.37) and the aforementioned assumptions, the likelihood of the observations can be written as

f (x | t ; η) = Q q=0 t q+1 n=tq+1 f (x n ; η q ) (4.39)
with t 0 0 and t Q+1 N . Note that from (4.38) and (4.39), it is also possible to write the joint distribution between the observations x and the parameter vector t, for a given parameter vector η, as f (x, t ; η) f (x | t ; η) π(t).

Finally, note that the random variable x n can be either absolutely continuous or discrete, depending on the application. In the following, we assume that it is continuous. However, the extension to the discrete case is straightforward.

Background on the covariance inequality

Let us consider an estimation problem, with an unknown parameter vector θ ∈ Θ which can be either deterministic, or random, or hybrid the latter case being considered in here. Let θ(x) an estimator of θ, i.e., a measurable function Ω → Θ. Let v(x, θ) a real measurable function, such that (i) the covariance matrix E{v(x, θ)v T (x, θ)} is positive denite, and (ii) the matrices E{ θ(x)v T (x, θ)} and E{θv T (x, θ)} have nite elements. The following matrix inequality then holds and is commonly referred to as the covariance inequality [START_REF] El Korso | Deterministic performance bounds on the mean square error for near eld source localization[END_REF]:

E θ(x) -θ θ(x) -θ T CV -1 C T (4.40) in which C E{ θ(x)-θ v T (x, θ)}, V E{v(x, θ)v T (x, θ
)} and the matrix inequality A B denotes the so-called Lowner partial ordering, i.e., the dierence A -B is a nonnegative matrix. In the sequel, we denote by S R the set of symmetric matrices, by S R + the set of nonnegative matrices and by S R ++ the set of positive denite matrices, all with size R × R.

At this point, it is worth noticing that, without any further assumptions on the vector function v(x, θ), the matrix C in (4.40) generally depends on θ(x), hence the right-hand side of (4.40) is not an interesting lower bound on the MSE of θ(x). However, for some well-chosen vector functions v(x, θ), and for some adequate set of estimators θ, it is possible to obtain a matrix C that does not depend on θ(x). In that case, the right-hand side of (4.40) is a lower bound on the MSE and applies to any estimator in the aforementioned set. We now explain how we dene the vector functions v(x, θ) to obtain our proposed lower bound.

The hybrid Cramér-RaoWeiss-Weinstein bound

In order to make the formulation of the bound more generic, we use notations that slightly dier from those in Section 4.2.1. Let us consider an R-dimensional hybrid unknown parameter vector θ = [θ T d , θ T r ] T belonging to the parameter space Θ ⊂ R R . Note that the term hybrid here means that θ d is deterministic, belonging to a subset Π d of R P , and whose true value is θ d , while θ r random, belonging to R Q (such that R = P +Q). For a given value of the deterministic parameter vector θ d , we dene a prior distribution π(θ r ; θ d ) for the random parameter vector θ r , whose support is denoted by Π r ⊂ R Q . Note that the prior might explicitly depend on θ d in some cases, see [START_REF] Basseville | Detection of Abrupt Changes, Theory and Application[END_REF]. Let us denote by f (x | θ r ; θ d ) the likelihood of the observations, so that the function

f (x, θ) f (x, θ r ; θ d ) = f (x | θ r ; θ d ) π(θ r ; θ d )
denotes the joint pdf of x and θ r parameterized by θ d . The following relations can be established between these notations and those from Section 4.2.1:

θ d ≡ η, θ r ≡ t, P ≡ L(Q + 1), Q ≡ Q, R ≡ L(Q + 1) + Q and π(θ r ; θ d ) = π(θ r ) ≡ π(t). Let us dene Θ {θ ∈ Θ | f (x, θ) > 0, for almost all x ∈ Ω}. The statistical expectation of a vectorial functional g θ d (x, θ r ) parameterized by θ d , w.r.t. the joint pdf f (x, θ) = f (x, θ r ; θ d ) is denoted by E x,θr;θ d {g θ d (x, θ r )}. We assume that, for any non-empty set (E, F ) ⊂ Ω × Π r , the integral E×F f (x, θ r ; θ d ) dx dθ r exists. Finally, let us denote by θ θ(x) = [ θT d , θT r ] T a joint (hybrid) estimator of [θ T d , θ T r ] T , i.e.
, θd θd (x) is an estimator of the deterministic parameter vector θ d and θr θr (x) is an estimator of a realization of the random vector θ r .

In order to obtain a lower bound on the estimation error of the hybrid parameter vector θ, the idea is to combine two dierent lower bounds w.r.t. θ d and θ r respectively. This results 102 CHAPTER 4. PERFORMANCE ANALYSIS WITH APPLICATIONS TO ARRAY PROCESSING in a hybrid lower bound for the estimation of the parameter vector θ. Such kind of lower bound has already been proposed in the literature [START_REF] Ren | Hybrid BarankinWeiss-Weinstein bounds[END_REF]. In this section, we propose to combine the (deterministic) Cramér-Rao bound with the (Bayesian) Weiss-Weinstein bound [START_REF] Weiss | A lower bound on the mean square error in random parameter estimation[END_REF]. A fully Bayesian version of this bound was rst proposed in a recursive form in [START_REF] Bell | Combined Cramér-Rao/Weiss-Weinstein bound for tracking target bearing[END_REF], and was we recently adapted to the o-line change-point problem for Poisson data in [START_REF] Bacharach | A Bayesian Lower Bound for Parameter Estimation of Poisson Data Including Multiple Changes[END_REF].

In order to derive the proposed hybrid bound, the vector function v(x, θ) : Ω × Θ → R M (with M ≥ R) is constructed in two parts since the estimation is hybrid: the P rst components of v(x, θ) are related to the deterministic parameters, while the J last components of v(x, θ) (P + J = M ) are related to the random parameters. Theoretically, J can be chosen such that J ≥ Q (which implies M ≥ R), but in the change-point framework, for similar reasons as explained in [208, Sect. II.B.] for the Barankin bound, taking J = Q (thus M = R) makes it possible to obtain closed-form expressions for the bound, so it can be computed in an ecient manner. Thus, for p = 1, . . . , P , we set

[v(x, θ)] p =      ∂ ln f (x, θ r ; θ d ) ∂θ d,p θ d =θ d , if θ ∈ Θ 0, if θ / ∈ Θ , (4.41)
in which θ d,p denotes the p-th component of the deterministic parameter vector θ d , and the derivatives are evaluated at the true value θ d of the parameter vector θ d . For q = P + 1, . . . , R,

we set [START_REF] Weinstein | A general class of lower bounds in parameter estimation[END_REF][START_REF] Ren | Hybrid BarankinWeiss-Weinstein bounds[END_REF] [

v(x, θ)] q =      f sq (x, θ r + h q ; θ d ) f sq (x, θ r ; θ d ) - f 1-sq (x, θ r -h q ; θ d ) f 1-sq (x, θ r ; θ d ) , if θ ∈ Θ 0, if θ / ∈ Θ , (4.42)
in which s q ∈ ]0 , 1[, and the vector h q ∈ R Q can be freely chosen in the set H θr h ∈ R Q θ r + h ∈ Π r . For convenience, we also dene the vector h q as the augmented vector h q with size R, completed with zeros, i.e., h q [0 T , h T q ] T . Without any further assumption, the matrix C in (4.40) still depends on θ(x), which is unwanted. We then make the two following additional assumptions:

1. We consider the class of estimators that are unbiased w.r.t. θ d , i.e.,

E x,θr;θ d { θ(x) -θ| θ d } = [0 T , d T ] T (4.43)
in which d is an arbitrary vector with size Q, independent of θ.

2. As for the classical CRB, we assume that, for any

θ d ∈ Π d , E x|θr;θ d ∂ ln f (x | θ r ; θ d ) ∂θ d = 0. (4.44)
Thus, assuming that both conditions (4.43) and (4.44) are satised, we obtain the following block-diagonal matrix C:

C = I 0 0 C 22 (4.45)
where I denotes the P × P identity matrix and the columns of C 22 are given, for 1 ≤ q ≤ Q and θ ∈ Θ , by

c q = E x,θr;θ d θ(x) -θ θ d f sq (x, θ r + h q ; θ d ) f sq (x, θ r ; θ d ) - f 1-sq (x, θ r -h q ; θ d ) f 1-sq (x, θ r ; θ d ) = [0 h q µ(s q , h q )] , (4.46)
in which we have dened, for h r ∈ H θr and s ∈ ]0 , 1[,

µ(s, h r ) Ω×R Q f s (x, θ r + h r ; θ d ) f 1-s (x, θ r ; θ d ) dx dθ r = E x,θr;θ d f s (x, θ r + h r ; θ d ) f s (x, θ r ; θ d ) . (4.47)
We can now dene the matrix V by splitting it into four blocks, as follows ii) (block V 22 ) for (q, q ) ∈ {P + 1, . . . , R} 2 and θ ∈ Θ , we have

V = V 11 V 12 V T 12 V 22 (4.
[V 22
] q,q = E x,θr;θ d f sq (x, θr + hq; θ d ) f sq (x, θr; θ d ) -f 1-sq (x, θrhq; θ d ) f 1-sq (x, θr; θ d ) × f s q (x, θr + h q ; θ d ) f s q (x, θr; θ d ) -f 1-s q (x, θrh q ; θ d ) f 1-s q (x, θr; θ d ) = ξ(sq, s q , hq, h q ) -ξ(sq, 1 -s q , hq, -h q ) -ξ(1 -sq, s q , -hq, h q ) + ξ(1 -sq, 1 -s q , -hq, -h q ) (4.50)

in which we have dened

ξ(α, β, h a , h b ) E x,θr;θ d f α (x, θ r + h a ; θ d ) f β (x, θ r + h b ; θ d ) f α+β (x, θ r ; θ d ) . (4.51)
This block corresponds to the usual Weiss-Weinstein lower bound for the estimation error associated with an estimator of θ r [START_REF] Weinstein | A general class of lower bounds in parameter estimation[END_REF]. Note that µ(s, h) = ξ(s, 0, h, 0). Thus, only the calculation of ξ(α, β, h a , h b ) is required for the determination of the elements in this block.

iii) (block V 12 ) for (p, q) ∈ {1, . . . , P } × {P + 1, . . . , R} and θ ∈ Θ This maximization operation can be very time consuming. For this reason we describe in the next section a method to make this computation feasible.

[V 12 ] p,q = E x,θr;θ d ∂ ln f (x, θ r ; θ d ) ∂θ r,p θ d f sq (x, θ r + h q ; θ d ) f sq (x, θ r ; θ d ) - f 1-sq (x, θ r -h q ; θ d ) f 1-sq (x, θ r ; θ d ) .

Practical computation of the Hybrid Cramér-RaoWeiss-Weinstein bound

Written as it is in (4.53), both obtaining closed-form expressions for W (H, s) and computing the supremum are infeasible tasks, even for simple problems. In order to overcome this diculty, two solutions are commonly adopted: (i) the set H Q θr is restricted to diagonal matrices H, i.e., the components of h q are all zero, except its q-th element [START_REF] Weinstein | A general class of lower bounds in parameter estimation[END_REF][START_REF] La Rosa | Barankin-type lower bound on multiple change-point estimation[END_REF]; and (ii) it is often noticed after extensive numerical experiments that the value s j = 0.5, for all j ∈ {1, . . . , Q}, leads to the tightest bound [START_REF] Weinstein | A general class of lower bounds in parameter estimation[END_REF][START_REF] Xu | Performances bounds on matched-eld methods for source localization and estimation of ocean environmental parameters[END_REF], which reduces the task dimensionality by a factor Q. It is also worth noticing that, in the framework of model (4.37), the set H Q θr is countable and nite, since its elements lie in the discrete set Θ t . This is also helpful to keep a reasonable complexity.

By setting s q = 0.5, ∀q (which is expected to yield the tightest lower bound), the supremum operation is computed over the set of matrices W W (H, s) ∈ S R ++ ; H ∈ H Q θr that is a (discrete) subset of S R , w.r.t. the Lowner partial ordering, denoted by . As a consequence, the uniqueness of the supremum might not be guaranteed, in particular if the set of upper bounds associated with W does not have a unique least element. If the supremum exists, it may belong to W, in which case it is the greatest element. On the contrary, when a greatest element does not exist, the set W may have several maximal elements, i.e., elements that are not lower than any other in the set. In such a case, a unique supremum does not exist and we have to nd a minimal element in the set of upper bounds of W. Formally, such a minimal upper bound B veries B W (upper bound). Moreover, if there exists a smaller element B such that B B W, then necessarily B = B (minimal). As already mentioned, B may not be unique without any additional constraint. If it is unique, then it is the unique least upper bound of W, that is its supremum.

Despite the previous comments, following [START_REF] La Rosa | Barankin-type lower bound on multiple change-point estimation[END_REF]Sect. III.D.], there is one way of obtaining a suitable minimal upper bound of W, such that it is dened in a unique manner, namely by dening it as the matrix associated with the Lowner-John hyper-ellipsoid. More explicitly, for each matrix A ∈ S R ++ , one can associate a centered hyper-ellipsoid, dened as the set E(A) y ∈ R R y T A -1 y ≤ 1 . Referring to [START_REF] La Rosa | Barankin-type lower bound on multiple change-point estimation[END_REF]Lemma 3], one can show that, for any positive denite matrices A and A , we have the equivalence: A

A i E(A) ⊆ E(A ). Thus, given a nite family of matrices {A i } i∈I , one can nd the minimum volume hyper-ellipsoid, denoted E(A LJ ), that covers the set of hyper-ellipsoids associated with each matrix A i , that is C E {E(A i ) ; i ∈ I}. The minimum volume ellipsoid E(A LJ ) is called the Lowner-John ellipsoid of C E , and one can show that the matrix A LJ is a minimal upper bound of the set {A i ; i ∈ I} (see [START_REF] La Rosa | Barankin-type lower bound on multiple change-point estimation[END_REF]Theorem 4]).

Finally, we compute a minimal upper bound of the set W by nding the matrix W LJ associated with the Lowner-John ellipsoid of W. This corresponds to a convex optimization problem [220, p. 411] that can be solved eciently using a semidenite programming procedure, such as those provided by the CVX package [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF].

PERFORMANCE ANALYSIS IN THE PRESENCE OF UNKNOWN CONTINUOUS AND DISCRETE PARAMETERS: APPLICATION TO CHANGE-POINTS ESTIMATION 105

We have presented the general expression of the bound, that is applicable to any hybrid estimation problem, and seen how to compute it eciently. We now give its expression for the problem exposed in Section 4.2.1.

Expressions of the Hybrid Cramér-RaoWeiss-Weinstein bound for the change-point problem

The derivation of the expression of the HRCWWB for the problem explained above is technical and rather complex. For sake of clarity, the main steps leading to these expressions are reported in Appendix. 6.5. Furthermore, in Appendix. 6.5.2, we explicit these general expressions to the case of Gaussian and Poisson distributions.

Numerical results

This section presents some simulation results that enable us to assess the tightness of the proposed bound. It is compared in terms of global mean square error (GMSE) with the so-called ML-MAP estimator, for the distributions discussed in section 6.5.2, namely the Gaussian and Poisson distribution. All the cases discussed in Appendix 6.5 were simulated with N = 100 observations and Q = 3 changes, and the GMSE of the ML-MAP estimator was obtained by computing the empirical MSE through 1000 Monte-Carlo simulations. At each Monte-Carlo run, the Q = 3

changes were generated according to the prior distribution (4.38).

ML-MAP estimator

The ML-MAP estimator can be used when some elements of the parameter vector are deterministic and the others are random variables [START_REF] Yeredor | The joint MAP-ML criterion and its relation to ML and to extended leastsquares[END_REF]. In our case, it is dened as follows:

(η MLMAP , tMLMAP ) arg max leading to ηMLMAP = η(x ; tMLMAP ).

Changes in the mean of a Gaussian distribution

In this classical change-point estimation problem, the parameter vector η contains the means µ q , q = 0, . . . , Q of each segment, and possibly the corresponding variances σ 2 q , q = 0, . . . , Q

(even if they remain constant), depending on whether they are assumed to be known or not. Figs.

4.1 and 4.2 display the GMSE and the associated HCRWWB for the means µ 1 , . . . , µ 4 of each segment, and for the change-point locations t 1 , t 2 , t 3 versus SNR, respectively. These quantities correspond to the diagonal elements of the matrices in (4.40). Variations of the variance estimates and their HCRWWBs are very similar to Fig. 4.1. Due to the lack of space, such graphs are not reported here. Note that the global root mean square error (GRMSE) of the estimated change-point locations was plotted instead of the GMSE, for a more relevant assessment of the gap with the bound. The SNR in the X-axis corresponds to snr m q,q dened in (6.82), quantifying the amount of change. More precisely, the Q changes generated for this experiment have all an equal amount of change, i.e., snr m 0,1 = snr m 1,2 = snr m 2,3 = snr m for a given snr m such that µ q = µ q-1 + (-1) q-1 snr m for q = 1, . . . , Q, according to (6.82). Fig. 4.1 clearly shows a threshold eect for the ML-MAP of the mean estimates, whose GMSEs move away from the HCRWWB for SNR lower than 15 dB. The threshold is lower regarding the last segment mean, around 10 dB. The HCRWWB renders this behavior very slightly, as can be seen from the tiny bulge in the shape of the bound, for SNRs between -10 dB and +5 dB. For SNRs higher than 15 dB, the GMSEs and HCRWWBs become much closer one to the other. The small gap remaining in between, at high SNR, comes from the fact that the CRB part in the HCRWWB actually corresponds, as already mentioned, to a so-called modied CRB (mCRB), in the sense of [START_REF] Gini | The modied Cramér-Rao bound in vector parameter estimation[END_REF]. Let us recall that the mCRB cannot be expected to be as tight as the classical CRB or the true Bayesian CRB, since the Fisher information is averaged over all the possible values of the nuisance parameter in this case t. With regards to estimation performance of change-point locations displayed in Fig. 4.2, both shapes of the GRMSEs and HCRWWBs highlight the existence of a non-information zone at low SNR, noticeable from the curve atness. In this area, the HCRWWB shows that the early changes are better estimated than the later ones. This is an eect of the chosen prior, which confers t q an increasing support as q tends to Q: according to (4.38), the rst change t 1 can fall between 1 and τ , while t 2 can fall between 2 (if t 1 = 1 = 1 and 2 = t 2 -t 1 = 1) and 2τ (if t 1 = 1 = τ and 2 = t 2 -t 1 = τ ). Outside this non-information zone, the dierence between GMSE and HCRWWB decreases: it is in the order of 9 or 10 samples for snr m = 0 dB while it is lower than 2 samples for snr m > 10 dB, and lower than 0.1 samples for snr m = 15 dB. Curves cannot be shown for higher SNRs because both empirical GRMSE and HCRWWB tend drastically to zero, due to the discrete nature of the change-point locations.

It is worth mentioning that the gap between the change-point location estimates and the bound is due to the discrete nature of these parameters. Indeed, discrete parameter estimation 4.2. PERFORMANCE ANALYSIS IN THE PRESENCE OF UNKNOWN CONTINUOUS AND DISCRETE PARAMETERS: APPLICATION TO CHANGE-POINTS ESTIMATION 107 is not the most usual estimation framework, and the classical convergence theorems (regarding the MLE for instance) no longer apply (see [START_REF] Choirat | Estimation in discrete parameter models[END_REF] for general considerations on discrete parameter estimation). This case is treated similarly to the previous one, with the only dierence that instead of changes in the mean, we study changes in the variance of the observations. Here, the parameter vector η includes at least the variances on each segment σ 2 q , q = 0, . . . , Q, and possibly the means µ q , q = 0, . . . , Q if they are unknown. The proposed simulations were obtained for unknown means, all equal to 1. Here, the SNR corresponds to snr v q,q dened in (6.82): it is the amount of change in terms of variance such that σ 2 q = snr v × σ 2 q-1 , for q = 1, . . . , Q, and for a given snr v . Estimated MSEs and the corresponding bounds for the variances and change locations are displayed in Figs.

Changes in the mean rate of a Poisson distribution

In the case of changes in the mean rate of a Poisson distribution, the parameter vector η includes the mean rates on each segment λ q , q = 0, . . . , Q. The results regarding mean rate and changepoint estimates are displayed in Figs. 4.5 and 4.6. Note that the SNR in these gures corresponds to the following denition (also used in [START_REF] Ferrari | Barankin lower bound for change points in independent sequences[END_REF][START_REF] La Rosa | Barankin-type lower bound on multiple change-point estimation[END_REF][START_REF] Bacharach | A Bayesian Lower Bound for Parameter Estimation of Poisson Data Including Multiple Changes[END_REF]) snr P q,q (λ q -λ q ) 2 /λ 2 q . Hence, for a given snr P , we have set λ q = λ q-1 (1 + √ snr P ), q = 1, . . . , Q. The main topic of this research is the estimation and performance analysis under non-standard conditions, related mainly to the non-Gaussian noise distribution, without assuming restrictive and classical regularity conditions and/or under misspecied models. As it can be seen from the above discussions, these studies are introduced mainly for array processing applications, however, it can be straightforwardly adapted to any other application as long as the observations are described by a parametric data model. Specically, the presented contributions can be categorized as follows • In the context of radar application, we have introduced ecient structured scatter matrix estimation algorithms for centered CES distributions. First, based on the EXIP principle, we have derived a two-step estimation procedure for the structured scatter matrix of CES distributions referred to as SESAME (StructurEd ScAtter Matrix Estimator). Then, we have conducted a theoretical asymptotic analysis to demonstrate the consistency, the unbiasedness, the asymptotic Gaussianity and eciency of the SESAME method. In addition, we have analytically expressed the coecients appearing in the SESAME's criterion for common CES distributions. Second, we have proposed a recursive application of the SESAME algorithm, which procedure leads to I-SESAME, possessing the same asymptotic performance but leading empirically to a faster convergence to the CRB than the initial SESAME.

In the same vein, we designed novel estimation procedure in the context of parameterized mean and variance under CG distribution clutter model without secondary data (i.e., without assuming the existing of target-free signals) and devised new subspace estimators using the MMSD framework, which minimizes the average natural distance between the true range space of interest and its estimate, under CG distribution.

• In the context of radio-interferometer application, we proposed a novel iterative scheme for parallel calibration of next generation radio astronomical arrays. The proposed algorithms take into account the direction dependent eects, which aect the apparent directions of the calibration sources, as well the parameter values variation across wavelength. The proposed PMCA (Parallel Multi-wavelength Calibration Algorithm) iteratively estimates the complex direction independent antenna gains and their noise powers, whereas, it jointly estimates the directions of the calibrators and their associated direction dependent gain. The parallelism in the PMCA ensures that the algorithm is scalable to large data sets, which can be paralleled over a number of nodes. The algorithm is set up in such a way that the use of data remains local, i.e., only the node on which a specic part of the data resides needs access to that data.

In the same vein, in Section 3.2, we carried out a robust analysis and algorithm design in the Jones matrices framework under CG distributed data.

• In the context of misspecied models, we formulated new Slepian-Bangs expressions for CES distributed data under misspecied model and thus we lled a theoretical and practical gap in the signal processing literature. Furthermore, we have shown that the proposed SB formulas encompass all the previously derived expressions as special cases. Moreover, these new SB formulas involved relatively slight modications with respects to the classical counterpart obtained for Gaussian and correctly specied data models with only an expectation of some scalar functions to derive or to evaluate numerically. The practical importance of the proposed expressions is in the fact that allow us to easily evaluate misspecied Cramér-Rao bounds for a lot of applications that are characterized by a non-Gaussian and heavy-tailed data behaviour along with a model misspecication.

• Finally, we derived closed-form expressions of the HCRWWB (Hybrid Cramér-Rao-Weiss-Weinstein bound) for signals including multiple change-points, as it is suited for an unknown parameter vector containing both continuous and discrete parameters. Numerical simulations assess the tightness of the bound in the cases of Gaussian (changes in the mean and the variance) and Poisson observations (changes in the mean rate). The proposed bound can be used as a reference to compare performance of various estimators for signals that include multiple change-points, in terms of MSE. In addition, a semidenite programming formulation to compute the tightest proposed HCRWWB was proposed. This latter consists of nding the unique minimum volume covering the set constituted by hyper-ellipsoid elements which are generated using the derived candidate HCRWWB matrices w.r.t. the so-called Loewner partial ordering.

Let us now take a look at some of the aforementioned work's perspectives: Perspectives w.r.t. joint parameterized mean and covariance/scatter matrix estimation: "On the way to abandon secondary data" Adaptive detection of targets embedded in a complex environment (involving strong clutter, jammers, etc.) is a major issue in statistical array processing. In this context, most of the commonly used detections processes are based on a 2-step approach, namely, performing

• The use of statistical methods to learn the unknown interference statistical parameters (covariance matrix, rank, etc.).

• The use of a statistical test (binary hypotheses) to perform the detection on a tested sample. This methodology has driven a signicant amount of research over the last decades, as numerous problems are to be dealt with: estimation and detection method design, robustness to various noise distributions and/or corruption by outliers, limited sample support, introduction of prior knowledge, and performance analysis. Two major hypotheses are however made in the previously described detection scheme:

• The secondary samples are assumed target-free.

• They are assumed homogeneous and stationary (following the same underlying distribution).

These hypotheses, inherently restrictive, are most of the time not realistic. In the following, two alternatives are presented. The rst alternative relies on a robust statistical-based method. And the second one aims to explore new methodologies, namely geometrical formulation, in order to alleviate the performance loss of statistical-based methods in mismatched hypothesis/situations. Bilodeau [223, p. 214]. The abandon of secondary data is challenging, compensation of this lack will be investigated throw prior knowledge on the collected data (e.g. prior p.d.f., structured mean and scatter matrix, etc...). Specically, this future work will be devoted to develop new estimations procedures to perform a simultaneous estimation of targets and clutter parameters in CES noise modeling. Consider that the target response p (Θ), of unknown parameters vector Θ (e.g. power, DoA...), is corrupted by an additive clutter plus noise of total covariance matrix Σ (Ω), with unknown parameters vector Ω. The parameterization of Σ is crucial since it allows for considering prior knowledge on the covariance structure (e.g., low rank, Toeplitz, persymetric...), hence developing estimation schemes that need few samples to be ecient. Prior information can come from physical considerations on the system, such as its geometry or physical models.

The framework of robust estimation considers the following loss function, that can be viewed as generalized log-likelihood

L (Θ, Ω) = M K K k=1 ln g (z k -p (Θ)) H Σ -1 (Ω) (z k -p (Θ)) + ln |Σ (Ω) |.
(5.1)

Robust constrained estimators are dened as the solutions of min

Θ,Ω L (Θ, Ω) s.t. Θ ∈ S θ Ω ∈ S ω (5.2)
where Θ ∈ S θ and Ω ∈ S ω are appropriate constrained sets corresponding to the considered signal and covariance modeling. These estimators are known to posses desirable practical properties, such as robustness to outliers, robustness to various signal distributions and asymptotic eciency.

In some cases, prior information is available on the parameters we want to estimate. For example, one can have acess to previous estimates of the target location, or knowledge on the covariance matrix model.Taking into account such prior information into the estimation process increases the performance of the corresponding adaptive processes and allows for using less samples. A natural way to do so is to use shrinkage (or regularization) methods, i.e., build estimators that are solutions of min Θ,Ω αL (Θ, Ω) + βL 0 (Θ, Ω)

s.t. Θ ∈ S θ Ω ∈ S ω (5.3)
where L 0 is a penalty function that shrinks the estimate towards prior values Θ 0 and Ω 0 . These methods have been recently studied in the context of robust estimation for the unstructured case (non parametered covariance and mean estimation) [START_REF] Ollila | Regularized M -estimators of scatter matrix[END_REF].

However, solving the aforementioned optimization problem is a challenging task since both the objective function and the constraint set are generally nonconvex. Hence, one has to develop specic algorithms depending on the considered covariance structure and signal characterization. Specically, the design of such algorithms will be tackled using the so-called Majorization-Minimization [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF] framework that has been shown to be an ecient tool in order to design fast and nearly optimal schemes in the context of structured covariance matrix estimation with signal free samples (p (Θ) = 0) [START_REF] Soloveychik | Group symmetric robust covariance estimation[END_REF][START_REF] Sun | Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions[END_REF]. As regards to the problem of estimating jointly the signal and covariance parameters, we note that [START_REF] Wang | Maximum likelihood estimation of compound Gaussian clutter and taget parameters[END_REF] tackled this problem by devising Expectation-Maximization type algorithms. Nevertheless, these algorithms do not take into account the covariance matrix structure. Moreover, they are restricted to a special, linear signal model (GMANOVA), thus cannot be directly applied to general MIMO-Radar models, nor to the DoD/DoA estimation problems, which are highly nonlinear and of interest to the practitioner [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF].

1.2-Robust Subspace Clustering for RADAR and Teledection

This perspective aims to exploit and extend recent advances in machine learning for the purpose of radar detection in complex interference environments. More specically, a focus on the use of robust Principal Component Analysis and robust union-of-subspaces recovery algorithms will be considered. This approach implies a radical change of paradigm and formulates a detection problem as a sparse regression one, rather than a statistical test. Indeed, the matrix of concatenated samples can be modeled as a sum of low-rank matrices (union of subspaces containing interferences) plus a sparse matrix times a dictionary of steering-vectors (representing the targets contribution). Recent advances allow to infer the recovery of this two components from a data matrix, as this problem is currently intensively investigated in the computer vision and machine learning community. Hence, this recovery allows us to build detection maps, as non-zero coecients of the sparse matrix indicates the presence of targets. Our preliminary work suggests that this track is promising (see Fig. 5.1). If this change in paradigm seems appropriate and promising, it raises new issues, notably when it comes to the integration of a priori knowledge and robust costs functions. We will thus pursue and extend these works for the considered radar detection applications. An eort will also be made to apply the developed algorithms on real data sets.

Perspectives w.r.t. the calibration of the next generation radio-interferometer

As explained above a robust calibration in the context of radio-interferometer is vital in order to obtain an accurate high resolution images. Using the same notation as in section 3.1, we recall that the sample covariance matrix reads

R(p f ) = G(g f )A(d f )Γ(γ f )Σ(σ f ) G(g f )A(d f )Γ(γ f )Σ(σ f ) H + R U (r f ) + Σ n (σ n f ) (5.4)
with f ∈ {f 1 , . . . , f J } and the unknown parameter vector is p

f = [d T f , g T f , (σ n f ) T , γ T f ]
T . Depending on the Landsdal's regime (cf., section 3), matrices G and Γ, are diagonal (regime 3), or full but structured (regime 4) [START_REF] Van Der Veen | Signal processing tools for radio astronomy[END_REF].

• From regime 3 to regime 4: The discussed method in Section 3.1 is adapted to the socalled regime 3. A direct extension to the regime 4, based on Section 3.1, will cause a heavy computational cost. I propose to handle this issue using partly calibrated arrays framework. Specically, each station will be calibrated separately using the same techniques as in Section 3.1, then, a consensus will be adopted to calibrate the overall arrays w.r.t. the gain based parameters. Regarding, the estimation of the DoA, which is the most heavy operation, I will rely on a mixed norm minimization which exploits block-sparsity and lowrank structure. This have been recently used successfully in [START_REF] Steens | Block-and rank-sparse recovery for direction nding in partly calibrated arrays[END_REF], in which the authors proposed an ecient implementation by nding a suitable compact equivalent problem reformulation. Finally, a projection technique [START_REF] Forster | Generalized rectication of cross spectral matrices for arrays of arbitrary geometry[END_REF] can be used in order to constrain the known geometry of the overall array.

• Robusteness: An attempt to take into account the heterogeneity of the observations in the context of the radio-interferometers was proposed in [START_REF] Kazemi | Robust radio interferometric calibration using the tdistribution[END_REF] by considering a Student-t distribution and in our paper using the CG distribution in [START_REF] Ollier | Robust calibration of radio interferometers in non-gaussian environment[END_REF]. Nevertheless, such work remains very marginal and conned only to the Jones matrix formulation. Consequently, it is natural to extend the approach proposed in Section 3.1 considering a CG distribution and more generally a CES distribution. Nevertheless, considering a CES distribution leads, very often, to no closed-form expression and thus weighs down the complexity. One alternative is to consider the following cost function

min p f ∈{f 1 ,...,f J } κ f (p f ) (5.5) s.t. C(p) = 0 (5.6) with κ f (p f ) = L R f (p f ) -Rf Ω f , Ω f = σ n f σ nT f -1 2 with Σ n f = diag(σ n f
) in which the operator L (.) will lead to a tradeo between robustness and low complexity. 116 CHAPTER 6. APPENDICES 6.1.2 Detail of calculus of Table 1 Complex Generalized Gaussian distribution The density generator function is given by:

g(t) = exp - t s b with s > 0, b > 0 and C m,g = sΓ(m)b -m/s π m Γ(m/s)
We obtain u ML (t) = s b t s-1 and ψ ML (t) = tu(t) = s b t s . To ensure the convergence of the ML scatter matrix estimator, the parameter s has to be strictly lower than 1. Thus, by recognizing a generalized Gamma distribution, we have

A ML = E ψ 2 ML (Q) = s 3 Γ(m/s)b 2+m/s R + x m-1+2s e - x s b dx = s 3 Γ(m/s)b 2+m/s Γ(m/s + 2)b 2+m/s s = m(m + s)
Complex W -distribution The density generator function can be written by:

g(t) = t s-1 e - t s b with s > 0, b > 0 and Cm,g = sΓ(m)b -(s+m-1)/s π m Γ s + m -1 s We obtain u ML (t) = s b t s-1 -(s -1)t -1 and ψ ML (t) = tu(t) = s b t s -(s -1).
To satisfy the Maronna's condition dened in [START_REF] Maronna | Robust M-estimators of multivariate location and scatter[END_REF], the parameter is necessarily s < 1. Again, the generalized Gamma distribution allows us to write

A ML = s Γ s + m -1 s b s + m -1 s R + s b x s -s + 1 2 x s+m-2 e - x s b dx = (m + 2s -1)(m + s -1) -2(s -1)(m + s -1) + (s -1) 2 = (m + s -1)s + m 2
Complex K-distribution The density generator function is:

g(t) = √ t ν-m K ν-m 2 √ νt with ν > 0 and C m,g = 2 ν (ν+m)/2 π m Γ(ν)
where K λ (•) denotes the modied Bessel function of the second kind. Straightforward, we get

u ML (t) = √ ν t K ν-m-1 2 √ νt K ν-m 2 √ νt and ψ ML (t) = √ νt K ν-m-1 2 √ νt K ν-m 2 √ νt
. Therefore, we obtain

AML = E ψ 2 ML (Q) = 2ν (ν+m)/2+1 Γ(ν)Γ(m) R + √ x m+ν K 2 ν-m-1 (2 √ νx) Kν-m (2 √ νx) dx = 1 2 m+ν Γ(ν)Γ(m) R + x m+ν+1 K 2 ν-m-1 (x) Kν-m (x)

dx

The integral expression can not be simplied into an explicit form.

Complex t-distribution

The reader is referred to [START_REF] Fortunati | Matched, mismatched, and robust scatter matrix estimation and hypothesis testing in complex t-distributed data[END_REF]. 4) and ( 7)]

√ N h ( r FP ) d → N (0, R ∞ )
where the asymptotic covariance R ∞ is given by R

∞ = 2 D h ΣD H h + D h ΩD T h (6.5)
with D h denotes the R-derivative ∂h ∂s of h at point s [START_REF] Kreutz-Delgado | The complex gradient operator and the CR-calculus[END_REF]. We note that

∂µ µ=µ e = 2 CRB -1 + CRB -1 = 4CRB -1
since CRB is real-valued. Finally, we obtain -√ N g N (µ e ) d → N 0, 4CRB -1 .

Complex Generalized Bingham Langevin distribution sampling

In this Appendix, we show how to sample a unitary random matrix X ∈ U

p CGBL (X) = p CGBL (X|C, K) = P p=1 p CVGBL (x p |c p , A p ) ∝ exp    P p=1 Re{c H p x p } + x H p A p x p    = P p=1 exp κRe{c H p x p } + x H p A p x p (6.7)
where, x p and c p stands for the p-th column of respectively X and C. Later, we emphasize briey the generation of the columns of X which are random unit norm vectors drawn from the (vector) CGBL distribution. All the upcoming framework shown below are the results of [START_REF] Ho | Simulation of the matrix bingham-von mises-sher distribution with applications to multivariate and relational data[END_REF] and [START_REF] Mardia | Directional Statistics[END_REF]. In the following framework, we dene the vector Bingham Langevin distribution [START_REF] Ho | Simulation of the matrix bingham-von mises-sher distribution with applications to multivariate and relational data[END_REF].

Next, we propose an algorithm for generating the columns of the unitary random matrix X using the vBL distribution. Finally, we deduce the generation of the unitary matrix X.

The vector Bingham Langevin distribution

The BL distribution [START_REF] Ho | Simulation of the matrix bingham-von mises-sher distribution with applications to multivariate and relational data[END_REF] is a probability distribution on the set of unitary real vectors which combines linear and quadratic terms. For instance, a given real unitary vector v follows the vBL distribution, i.e., v ∼ vBL(c, A). Its pdf reads as follows p vBL ∝ exp{c T v + v T Av} where A is a symmetric matrix and c is a real vector. As described in [START_REF] Ho | Simulation of the matrix bingham-von mises-sher distribution with applications to multivariate and relational data[END_REF], a Markov chain Monte Carlo method (MCMC) for sampling from p vBL is derived. All the details about the generation of the vBL(c, A) are available in [START_REF] Ho | Simulation of the matrix bingham-von mises-sher distribution with applications to multivariate and relational data[END_REF]. The vBL can be generated from the acceptance-rejection scheme [START_REF] Kent | A new method to simulate the bingham and related distributions in directional data analysis with applications[END_REF].

The vector complex vector generalized Bingham Langevin distribution

Let's start rst with dening the relation between the CVGBL distribution and the vBL distribution. Based on [START_REF] Mardia | Directional Statistics[END_REF], for a given complex unitary random vector v such that v ∼ CVGBL(c, A)

and its pdf reads as: for the p-th antenna and J i,q (θ i,q ) = q i 1 q i 2 q i 3 q i 4 for the q-th antenna, i.e., θ i,p = [p i 1 , p i 2 , p i 3 , p i 4 ] T and θ i,q = [q i 1 , q i 2 , q i 3 , q i 4 ] T . Using these latter notation, we obtain

p CVGBL (v) ∝ exp Re{c H v} + v H Av
u ipq (θ ip ) = Σ iq θ ip . (6.8)
where

Σ i,q =     α i,q β i,q 0 0 0 0 α i,q β i,q γ i,q ρ i,q 0 0 0 0 γ i,q ρ i,q     (6.9) in which α i,q = q * i 1 c i 1 + q * i 2 c i 3 , β i,q = q * i 1 c i 2 + q * i 2 c i 4 , γ i,q = q * i 3 c i 1 + q * i 4 c i 3 and ρ i,q = q * i 3 c i 2 + q * i 4 c i 4 .
We also obtain

u iqp (θ ip ) = Υ iq θ * ip . (6.10)
where

Υ i,q =     λ i,q µ i,q 0 0 ν i,q ξ i,q 0 0 0 0 λ i,q µ i,q 0 0 ν i,q ξ i,q     (6.11) in which λ i,q = q i 1 c i 1 + q i 2 c i 2 , µ i,q = q i 1 c i 3 + q i 2 c i 4 , ν i,q = q i 3 c i 1 + q i 4 c i 2 and ξ i,q = q i 3 c i 3 + q i 4 c i 4 .
Finally, the cost function

φ i (θ ip ) = M q=1 q>p w ipq -u ipq (θ ip ) H (β i τ pq Ω pq ) -1 w ipq -u ipq (θ ip ) + M q=1 q<p w iqp -u iqp (θ ip ) H (β i τ qp Ω qp ) -1 w iqp -u iqp (θ ip ) + Constant.
(6.12) can be re-written as

φ i (θ i,p ) = w i,p -u i,p (θ i,p ) H A i,p w i,p -u i,p (θ i,p ) + wi,p -ũi,p (θ i,p ) H Ãi,p wi,p -ũi,p (θ i,p ) + Constant (6.13) where w i,p = [w T i,p(p+1) , . . . , w T i,pM ] T , u i,p (θ i,p ) = [u T i,p(p+1) (θ i,p ), . . . , u T i,pM (θ i,p )] T and A i,p = bdiag{β i τ p(p+1) Ω, . . . , β i τ pM Ω} -1 . Furthermore, we have wi,p = [w * T i,1p , . . . , w * T i,(p-1)p ] T , ũi,p (θ i,p ) = [u * T i,1p (θ i,p ), . . . , u * T i,(p-1)p (θ i,p )] T and Ãi,p = bdiag{β i τ 1p Ω * , . . . , β i τ (p-1)p Ω * } -1 .
We make use of (6.8) in what follows where

u i,p (θ i,p ) =    u i,p(p+1) (θ i,p ) . . . u i,pM (θ i,p )    =    Σ i,p+1 θ i,p . . . Σ i,M θ i,p    = Σ i θ i,p (
Σ i = [Σ T i,p+1 , • • • , Σ T i,M ] T . Likewise, we use (6.10) in ũi,p (θ i,p ) =    u * i,1p (θ i,p ) . . . u * i,(p-1)p (θ i,p )    =    Υ * i,1 θ i,p . . . Υ * i,p-1 θ i,p    = Υ i θ i,p (6.15) in which Υ i = [Υ * T i,1 , • • • , Υ * T i,p-1 ]
T . Inserting (6.14) and (6.15) into (6.13) and taking the derivative w.r.t. θ i,p leads to the expressions in (3.53), using the fact that A i,p and Ãi,p are Hermitian.

6.4 SB-formula for misspecied models: Proofs and relationship to previous results 6.4.1 Some considerations on the expectation of (4.15) and (4.16)

The aim of this appendix is to show under which conditions the expectations E p {V ij (θ 0 )} and E p {H ij (θ 0 )} could be evaluated in closed form. Keeping in mind the Stochastic Representation of a CES distributed random vector x l given in (4.2), let us dene the vector z l = Tu l . Then, by recalling that r 0 l µ l -γ 0 l (see (4.32)), the term ∂G 0 l ∂θ i in (4.22) can be rewritten as:

∂G 0 l ∂θ i = - ∂γ 0 l ∂θ i H Π -1 0 (Q 1/2 l z l + r 0 l ) -(Q 1/2 l z l + r 0 l ) H Π -1 0 ∂γ 0 l ∂θ i -(Q 1/2 l z l + r 0 l ) H S 0 i (Q 1/2 l z l + r 0 l ), (6.16)
where S 0 i is dened in (4.23). Let us dene the following vector and scalar quantities:

h il = S 0 i r 0 l + Π -1 0 ∂γ 0 l ∂θ i (6.17)
and

A il = (r 0 l ) H S 0 i r 0 l + (r 0 l ) H Π -1 0 ∂γ 0 l ∂θ i + ∂γ 0 l ∂θ i H Π -1 r 0 l . (6.18)
Then (6.16) becomes

∂G 0 l ∂θ i = -Q l z H l S 0 i z l -Q 1/2 l z H l h il -Q 1/2 l h H il z l -A il . (6.19)
Note that, according to the Stochastic Representation Theorem [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF], we have that Q

1/2 l = d R and Q l = d Q. Consider now the term E p φ(G 0 l ) ∂G 0 l ∂θ i
. By exploiting the expression in (6.19), we have that:

E p φ(G 0 l ) ∂G 0 l ∂θ i = -E p {φ(G 0 l )Q l z H l S 0 i z l } -E p {φ(G 0 l )Q 1/2 l z H l h il } -E p {φ(G 0 l )Q 1/2 l h H il z l } -A il E p {φ(G 0 l )}. (6.20)
It must be noted now that, since G 0 l has been dened in (4.20) as,

G 0 l (x l -γ 0 l ) H Π -1 0 (x l -γ 0 l ) = (Q 1/2 l z l + r 0 l ) H Π -1 0 (Q 1/2 l z l + r 0 l ) (
E p φ(G 0 l )φ(G 0 m ) ∂G 0 l ∂θ i ∂G 0 m ∂θ j for E p {V ij (θ 0 )} in (4.21), E p φ (G 0 l ) ∂G 0 l ∂θ j ∂G 0 l ∂θ i and E p φ(G 0 l ) ∂ 2 G 0 l ∂θ i ∂θ j for E p {H ij (θ 0 )} in (4. 24
). There are, however, two important cases in which φ(G 0 l ) is statistically independent of Q 1/2 l (and consequently of Q l ) and z l : these are the Scenario 1 and the Scenario 2 which are discussed in the following Appendices.

Proof of the SB formulas for the Scenario 1

As discussed in Section 4.1.2, in the Scenario 1, the true marginal pdf is an unspecied CES pdf p X (x l ) = CES N (µ l , Σ, g), while the assumed one is a parametric complex Gaussian pdf f X (x l ; θ) = CN (γ l (θ), Π(θ)). As a consequence, φ(G l (θ)) ≡ -1, ∀θ ∈ Θ and then it trivially satisfy the independence condition discussed in 6.4.1. In order to obtain a closed form expression of the matrix B(θ 0 ), the second term in (4.30) can be rewritten as:

L l=1 L m=1 E p ∂G 0 l ∂θ i ∂G 0 m ∂θ j = L l=1 E p ∂G 0 l ∂θ i ∂G 0 l ∂θ j + L l=1 L m=1, m =l E p ∂G 0 l ∂θ i ∂G 0 m ∂θ j . (6.22)
By using (6.19) and the Stochastic Representation Theorem, when l = m we have that:

E p ∂G 0 l ∂θ i ∂G 0 m ∂θ j =E 2 Q {Q}E p z H l S 0 i z l E p z H m S 0 j z m + E Q {Q} A jm E p z H l S 0 i z l + A il E p z H m S 0 j z m + A il A jm , (6.23)
where A il is dened in (6.18) and where we used the facts that z l and z m are independent and E p {z l } = E p {z m } = 0. When l = m, the term in (6.23) has to be recast as:

E p ∂G 0 l ∂θ i ∂G 0 l ∂θ j =E Q {Q 2 }E p z H l S 0 i z l z H l S 0 j z l + E Q {Q} A jl E p z H l S 0 i z l + A il E p z H l S 0 j z l + E p z H l h il h H jl z l + E p h H il z l z l h jl + A il A jl , (6.24)
where h il is dened in (6.17) and where we used the fact that the third-order moments of u l (and then of z l ) vanish (see Lemma 1 in [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]). Moreover, as a consequence of the circularity property of z l [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF], the following equalities hold:

E p {z l z T l } = E p {z * l z H l } = 0, (6.25)
E p {z H l h il z H l h jl } = E p {h T il z * l z H l h jl } = 0, (6.26) E p {h H il z l h H jl z l } = E p {h H il z l z T l h * jl } = 0. (6.27)
From the properties of the trace operator, we have that:

E p z H l S 0 i z l = N -1 tr(S 0 i Σ), (6.28)
E p z H l h il h H jl z l = N -1 tr(h il h H jl Σ), (6.29)
6.4. SB-FORMULA FOR MISSPECIFIED MODELS: PROOFS AND RELATIONSHIP TO PREVIOUS RESULTS 123

E p h H il z l z l h jl = N -1 tr(h jl h H il Σ), (6.30)
while, by exploiting exactly the same procedure used in [START_REF] Greco | Cramér-Rao lower bounds on covariance matrix estimation for complex elliptically symmetric distributions[END_REF] and [START_REF] Greco | Maximum likelihood covariance matrix estimation for complex elliptically symmetric distributions under mismatched conditions[END_REF], it can be shown that:

E p z H l S 0 i z l z H l S 0 j z l = tr(S 0 i Σ)tr(S 0 j Σ) + tr(S 0 i ΣS 0 j Σ) N (N + 1)
.

(6.31)

Another useful relation can be obtained by exploiting the equality in (4.27). In fact, by combining (6.20) with (4.27), we have that:

E Q {Q}E p {z H l S 0 i z l } = -tr(P 0 i ) -A il , i = 1, . . . , d, (6.32)
and then, from (6.28) and from the identiablity condition for CES distribution assumed in this section, i.e. E Q {Q} = N , we obtain:

A il = -tr(P 0 i ) -tr(S 0 i Σ), i = 1, . . . , d. (6.33)
Using (6.32) and the condition E Q {Q} = N (and consequently E 2 Q {Q} = N 2 ), the term in (6.23) can be easily evaluated as:

E p ∂G 0 l ∂θ i
∂G 0 m ∂θ j = tr(P 0 i )tr(P 0 j ).

(6.34)

Similarly, using (6.33), (6.29), (6.30), (6.31) and the condition E Q {Q} = N , the term in (6.24) can be expressed as:

E p ∂G 0 l ∂θ i ∂G 0 l ∂θ j =tr(P 0 i )tr(P 0 j ) + tr(h il h H jl Σ) + tr(h jl h H il Σ) + E Q {Q 2 } N (N + 1) -1 tr(S 0 i Σ)tr(S 0 j Σ) + tr(S 0 i ΣS 0 j Σ) (6.35)
Finally, by combining all the previous results, we have that:

L l=1 L m=1 E p ∂G 0 l ∂θ i ∂G 0 m ∂θ j =L 2 tr(P 0 i )tr(P 0 j ) + L l=1 tr(h il h H jl Σ) + tr(h jl h H il Σ) + L E Q {Q 2 } N (N + 1) -1 tr(S 0 i Σ)tr(S 0 j Σ) + L E Q {Q 2 } N (N + 1) tr(S 0 i ΣS 0 j Σ), (6.36)
from which (4.31) follows immediately. Regarding the calculation of the matrix A(θ 0 ) in (4.34), no particular simplication can be made. Specically, the term E p ∂ 2 G 0 l ∂θ i ∂θ j in (4.33) has to be obtained through direct calculation from (6.16). Since the derivation is long, tedious and does not add any insightful considerations about the problem at hand, we decided to not report it here.

Proof of the SB formulas for the Scenario 2

In Scenario 2 we suppose that the true marginal pdf is given by a parametric CES distribution such that p X (x l ) p X (x l ; θ) = CES N (µ l , Σ, g), where µ l = γ l ( θ) and Σ = Π( θ) for a given θ ∈ Θ. The assumed pdf is itself a parametric CES distribution that share the same parametrization of the true one but may have a dierent density generator, i.e. f X (x l ; θ) = CES N (γ l (θ), Π(θ), w), θ ∈ Θ, and possibly g(t) = w(t), ∀t ∈ R + . In order to guarantee the 124 CHAPTER 6. APPENDICES correct identiability of the true and the assumed CES distributions, as before, we may impose a constraint of both the modular variate, that is: E p {(x l -µ l ) H Σ -1 (x l -µ l )} E p {(x l -γ l ( θ)) H Π( θ) -1 (x l -γ l ( θ))} = E Q {Q l } = N, (6.37) E f θ {(x l -γ l (θ)) H Π(θ) -1 (x l -γ l (θ))} = N, ∀θ ∈ Θ. (6.38) Let us suppose now that, for a given true pdf p X (x l ) p X (x l ; θ) and for a given assumed pdf f X (x l ; θ), the pseudo-true parameter vector θ 0 equates the true parameter vector θ, and consequently, from (4. Under this assumption, we have that: where γ 0 l γ l (θ 0 ) and Π 0 Π(θ 0 ). It is worth to highlight that the equality chains in (6.40) and (6.41) hold true if and only if θ 0 = θ. We are quite condent that, in the context of Scenario 2 (i.e. when the misspecication is only due to a wrong assumption on the density generator), the equality θ 0 = θ always holds true.

G 0 l (x l -γ 0 l ) H Π -1 0 (x l -γ 0 l ) = (x l -γ l ( θ)) H Π( θ) -1 (x l -γ l ( θ)) = (x l -µ l ) H Σ -1 (x l -µ l ) Q l = d Q = R 2 ,
As a consequence of (6.40), the expectation of functions of G 0 l , say h(G 0 l ), can be explicitly derived as:

E Q {h(G 0 l )} ∞ 0 h(t)p Q (t)dt = δ -1
N,g ∞ 0 h(t)t N -1 g(t)dt = E Q {h(Q)}, ∀l ∈ N.

(6.42)

Since such expectation does not depend on the index l and on the unknown parameter vector θ, in the rest of this Appendix we always indicate E Q {h(G 0 l )} simply as E Q {h(Q)}.

As in above for the Scenario 1, in order to obtain a closed form expression of the matrix B(θ 0 ), the second term in (4.28) can be rewritten as: Let us now recall from [START_REF] Besson | On the sher information matrix for multivariate elliptically contoured distributions[END_REF] the following three equalities:

E p {u H l P i u l u H l P j u l } = tr(P i )tr(P j ) + tr(P i P j ) N (N + 1) for some vector a and hermitian matrix D.

By substituting (6.44) in (6.43) and by exploiting the equalities (6.45), (6.46) and (6. + L E Q {Q 2 φ 2 (Q)} N (N + 1) (tr(P i )tr(P j ) + tr(P i P j )) + L(L -1)

E 2 Q {Qφ(Q)} N 2
tr(P i )tr(P j ).

(6.48)

Let us now impose the equality in (4.27). In particular, we have that:

L l=1 E p φ(G 0 l ) ∂G 0 l ∂θ i = -E Q {Qφ(Q)} L l=1 E p ∂γ 0 l ∂θ i H Π -1/2 0 u l + u H l Π -1/2 0 ∂γ 0 l ∂θ i -E Q {Qφ(Q)} L l=1 E p u H l P 0 i u l = -LN -1 E Q {Qφ(Q)}tr(P i ), (6.49)
where we used the fact that the fact that E p {u l } = 0 while the last equality follows from (6.46).

Then, by exploiting the equality (4.27), we obtain the following relation:

E Q {Qφ(Q)} = -N, (6.50)
where φ(t) = w (t)/w(t) and w(t) is the density generator of the assumed CES distribution f X (x l ; θ). Finally, by substituting (6.50) in (6.48), and then replacing the obtained term in (4.28), we get the colsed form expression of the matrix B(θ 0 ) given in (4.35). The evaluation of the matrix A(θ 0 ) in (4.36) follows directly from the direct calculation of the terms E p φ (G 0 l )

∂G 0 l ∂θ j ∂G 0 l ∂θ i and E p φ(G 0 l ) ∂ 2 G 0 l ∂θ i ∂θ j
. In particular, all the derivatives have to be evaluated from (4.22), while, to evaluate the expectation operator w.r.t. the true distribution p X (x), one has to use the Stochastic Representation in (6.41) by keeping in mind that θ = θ 0 . Since this calculation is tedious and not informative, we decided to not report it here.

The SB formula under correctly specied CES models

The SB formula for correctly specied CES model has been derived in [START_REF] Besson | On the sher information matrix for multivariate elliptically contoured distributions[END_REF]. Using the formalism introduced in above, two parametric CES models are said to be correctly specied if there exist a vector θ ∈ Θ, such that the assumed CES distribution in (4.9) equates the true CES distribution in (4.7). More formally, the CES model f X (x; θ) is said to be correctly specied if there exists θ ∈ Θ, such that f X (x; θ) = p X (x) and, in particular, γ l ( θ) = µ l , Π( θ) = Σ and g(t) ≡ w(t).

As proved in [START_REF] Vuong | Cramér-rao bounds for misspecied models[END_REF], under correctly specied model, we have that θ = θ 0 and B( θ) = -A( θ), where B( θ) is the classical FIM evaluated at the true parameter vector θ.

In the following, we will show that the SB formula derived in [START_REF] Besson | On the sher information matrix for multivariate elliptically contoured distributions[END_REF] can be considered as a special case of the one obtained in Section 4.1.2. In fact, according to the matched model assumption, we can dene the true model as p X (x) = f X (x; θ) = CES N (γ l ( θ), Π( θ), w), while the assumed parametric model is f X (x; θ) = CES N (γ l (θ), Π(θ), w) with θ ∈ Θ. With this in mind, we can exploit the result in eq. ( 11) of [START_REF] Besson | On the sher information matrix for multivariate elliptically contoured distributions[END_REF], that is E Q {Q l φ(Q l )} = -N , where φ(t) g (t)/g(t). Finally, as shown in [START_REF] Vuong | Cramér-rao bounds for misspecied models[END_REF], we have that B(θ 0 ) = B( θ) = -A( θ) = FIM( θ).

Then by exploiting (4.35), we obtain: is exactly the same as the FIM given in eq. ( 20) of [START_REF] Besson | On the sher information matrix for multivariate elliptically contoured distributions[END_REF].

[B( θ)] i,j = [FIM( θ)] i,j = 2 N E Q {Q φ2 (Q)} L l=1 Re ∂ γl ∂θ i H Π-1 ∂ γl ∂θ j + L E Q {Q 2 φ2 (Q)} N (N +
6.4.5 The SB formulas for scatter matrix estimation under misspecication of the density generator

In [START_REF] Fortunati | The Misspecied Cramér-Rao bound and its application to scatter matrix estimation in complex elliptically symmetric distributions[END_REF][START_REF] Greco | Maximum likelihood covariance matrix estimation for complex elliptically symmetric distributions under mismatched conditions[END_REF], SB formulas for the scatter matrix estimation in CES distributed vectors under misspecication of the density generator have been proposed. It is easy to verify that this scenario is a special case of the more general Scenario 2 discussed in Section 4.1.2. In particular, according to [START_REF] Fortunati | The Misspecied Cramér-Rao bound and its application to scatter matrix estimation in complex elliptically symmetric distributions[END_REF][START_REF] Greco | Maximum likelihood covariance matrix estimation for complex elliptically symmetric distributions under mismatched conditions[END_REF], the dataset x is considered to be composed of L independent, zero mean snapshots distributed as x l ∼ CES N (µ, Σ, g) where µ ≡ 0 and Σ = Π( θ), for a given θ ∈ Θ that is the true parameter vector. For each snapshot, we assume a marginal pdf f X (x l ; θ) = CES N (0, Π(θ), w),

i.e. we misspecied the density generator. Since we aim at nding SB formulas for the estimation of the scatter matrix itself, we have that θ vecs(Π(θ)) ≡ vecs(Π), where vecs is the operator that maps a symmetric N × N matrix Π in a N (N + 1)/2-dimensional vector whose entries are the elements of the upper (or lower) sub-matrix of Π. A similar notation holds for the true parameter vector, and in particular θ = vecs( Π). Note that, in the following, we assume that the entries of θ (and then of the scatter matrix) are real numbers. Finally, the results in [START_REF] Fortunati | The Misspecied Cramér-Rao bound and its application to scatter matrix estimation in complex elliptically symmetric distributions[END_REF][START_REF] Greco | Maximum likelihood covariance matrix estimation for complex elliptically symmetric distributions under mismatched conditions[END_REF] can be readily derived by (4.35) and (4.36) by posing L = 1,

∂γ 0 l ∂θ i ≡ 0 and E Q {Qφ(Q)} = -N : [B(θ 0 )] i,j = E Q {Q 2 φ 2 (Q)} N (N + 1)
tr(P 0 i )tr(P 0 j ) + tr(P 0 i P 0 j ) -tr(P 0 i )tr(P 0 j ), (6.52) [A(θ 0 )] i,j = E Q {Q 2 φ (Q)} N (N + 1) tr(P 0 i P 0 j ) + tr(P 0 i )tr(P 0 j ) -tr(P 0 i P 0 j ), (6.53) where, in this case, the matrix P 0 i (or P 0 j ) can be expressed as P 0 i = Π -1 0 A i = Π-1 A i where Π 0 is the matrix that minimizes the KLD between the true and the assumed distributions, and then, due to the Assumption A2, it is equal to the true covariance matrix, i.e. Π 0 = Π = Π( θ) = Σ. A i is a matrix dened as A i ∂Π(θ) ∂θ i ≡ ∂Π ∂θ i .

6.4. SB-FORMULA FOR MISSPECIFIED MODELS: PROOFS AND RELATIONSHIP TO PREVIOUS RESULTS 127 6.4.6 The SB formulas for misspecied Gaussian models

In this subsection, we briey show how to obtain the SB formulas provided in [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF] by using the general results discussed here for the Scenario 1. The SB formulas for misspecied Gaussian models can, in fact, be derived as a special case of the ones given in (4.31) and (4.34) by positing as true distribution the parametric complex Gaussian model, i.e. p X (x l ) = CES N (µ l , Σ, exp(-t)) = CN (µ l , Σ), while the assumed marginal distribution is still given by f X (x l ; θ) = CES N (γ l (θ), Π(θ), exp(-t)) = CN (γ l (θ), Π(θ)). Note that the true and the assumed density generators are equal and it is simply given by g(t) = w(t) = exp(-t). This fact can be used to evaluate the term E Q {Q 2 } as:

E Q {Q 2 } = ∞ 0 exp(-t)t N +1 δ -1 N,g = -[exp(-t)t N +1 δ -1 N,g ] ∞ 0 + (N + 1) ∞ 0
t N exp(-t)δ -1 N,g dt.

(6.54)

Since, [exp(-t)t N +1 δ -1 N,g ] ∞ 0 = 0, we have that

E Q {Q 2 } = (N + 1) ∞ 0 t N exp(-t)δ -1 N,g dt = -[exp(-t)t N δ -1 N,g ] ∞ 0 + N (N + 1)
∞ 0 t N -1 exp(-t)δ -1 N,g dt = N (N + 1), (6.55) in which we used the fact that ∞ 0 t N -1 exp(-t)δ -1 N,g dt = 1. Consequently, the entries of the matrix B(θ 0 ) in (4.31) can be readily expressed as:

[B(θ 0 )] i,j = 2 L l=1 Re Π 0 i Π -1 0 r 0 l + ∂γ 0 l ∂θ i H Π -1 0 ΣΠ -1 0 Π 0 j Π -1 0 r 0 l + ∂γ 0 l ∂θ j + + Ltr(S 0 i ΣS 0 j Σ), (6.56)
where r 0 l µ l -γ l (θ 0 ).

The generalized Slepian formulas

In this particular case, we have that all the derivatives of Π w.r.t. θ are zero. In particular, S i = S j = P i = P j = P ij = 0. Then, by posing L = 1 as in [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF], the entries of the matrices B(θ 0 ) and A(θ 0 ) can be obtained from (6.56) and (4.34) respectively, as:

[B(θ 0 )] i,j = 2Re ∂γ 0 ∂θ i H Π -1 0 ΣΠ -1 0 ∂γ 0 ∂θ j (6.57)

and

[A(θ 0 )] i,j = -2Re ∂γ 0 ∂θ i H Π -1 0 ∂γ 0 ∂θ j -2Re (r 0 ) H Π -1 0 ∂ 2 γ 0 ∂θ i ∂θ j , (6.58)
where r 0 µγ(θ 0 ). It is immediate to verify that these two expression are exactly the same of eq. ( 38) and [START_REF] Mardia | Directional Statistics[END_REF] derived in [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF].

The generalized Bangs formulas

Here we suppose that the mean value γ is independent of the parameter vector θ, consequently we have that ∂γ ∂θ i = 0 and ∂ 2 γ ∂θ i ∂θ j = 0. Then, if L = 1, the entries of the matrix B(θ 0 ) in (6.56) can be readily expressed as:

[B(θ 0 )] i,j = 2Re Π 0 i Π -1 0 r 0 H Π -1 0 ΣΠ -1 0 Π 0 j Π -1 0 r 0 l + tr(S 0 i ΣS 0 j Σ), (6.59) 128 CHAPTER 6. APPENDICES from which one deduces straightforwardly eq. ( 44) of [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF] with r 0 µγ(θ 0 ). Finally, let us now derive the matrix A(θ 0 ) for the particular scenario at hand. From (4.34), we get:

[A(θ 0 )] i,j =tr(P 0 i P 0 j -P 0 ij ) -r 0 l H Π -1/2 0 (P 0 i P 0 j + P 0 j P 0 i -P 0 ij )Π -1/2 0 r 0 l + + tr (P 0 i P 0 j + P 0 j P 0 i -P 0 ij )Π which is the same as that given in eq. ( 46) of [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecication[END_REF]. . We rst do it in the general case, that is for any distribution of the observations f (.; η q ) (assumed to be known). We proceed in the following way: rst, we derive the block V 11 as given by (4.49). Then, we compute the block V After some tedious calculus, we show that the matrix block V 11 is block-diagonal, i.e., can be written

V 11 = diag τ + 1 2 F (η 0 ), . . . , τ + 1 2 F (η Q-1 ), N - Q(τ + 1) 2 F (η Q ) (6.62)
in which F (η q ) is the (L × L) Fisher information matrix for the parameter vector η q ∈ R L , and for one observation in the q-th segment, i.e., for q ∈ {0, . . . , Q} and ( , ) ∈ {1, . . . , L} 2 , F (η q ) , = -E x;η q ∂ 2 ln f (x ; η) ∂η q, ∂η q, = -Ω ∂ 2 ln f (x ; η q ) ∂η q, ∂η q, f (x ; η q ) dx. Invoking the the block diagonal structure of the block matrix V 11 leads to a decoupling between parameter vectors η q associated with each segment. As mentioned earlier, this block corresponds to the Weiss-Weinstein bound on the MSE of the parameter vector t for given parameters η 0 , . . . , η Q . In other words, we follow exactly the same methodology as in [START_REF] Bacharach | Weiss-weinstein bound on multiple change-points estimation[END_REF]. The main result of that paper states that the block V 22 is tridiagonal, i.e., for any (q, q ) ∈ {1, . . . , Q} 2 such that |q -q | > 1, we have [V 22 ] q,q = 0.

(6.64)

The diagonal terms of V 22 correspond to the numerator of ( 26) in [START_REF] Bacharach | Weiss-weinstein bound on multiple change-points estimation[END_REF], i.e., are dened as [V 22 ] q,q = u D (τ , h q ) • ρ |hq| q (ε hq (2s q )) + ρ |hq| q (ε hq (2s q -1)) -2 u D (τ , 2h q ) • ρ 2|hq| q (ε hq (s q )) (6.65) in which, i) for τ ∈ N, q ∈ {1, . . . , Q} and h q ∈ Z,

u D (τ , h q )      (τ -|hq|) 2 τ 2 ifq < Qand|h q | < τ τ -|h Q | τ ifq = Qand|h Q | < τ 0 and|h q | τ , (6.66)
ii) for q ∈ {1, . . . , Q} and s ∈ ]0 , 1[, ρ q (s) Ω f s (x ; η q-1 ) f 1-s (x ; η q ) dx (6.67)

and iii) for h ∈ Z and s ∈ ]0 , 1[, ε h (s) s ifh > 0 1 -s ifh < 0. From the terms in the block V 22 , we can compute the block C 22 since we have seen that µ(s, h) = ξ(s, 0, h, 0). In addition, due to the structure chosen for vectors h q , namely only its q-th component is nonzero, the block matrix C 22 is diagonal. The resulting expression of µ(s, h) corresponds to equation [START_REF] Ottersten | Covariance matching estimation techniques for array signal processing applications[END_REF] in [START_REF] Bacharach | Weiss-weinstein bound on multiple change-points estimation[END_REF]. After plugging this expression into (4.46), we obtain [C 22 ] q,q = h q u D (τ , h q ) ρ |hq| q (ε hq (s q )). (6.69) Finally, the remaining nonzero terms in the block V 22 are the superdiagonal ones (which equal, by symmetry, the subdiagonal ones). Referring to [209, equations (31) and [START_REF] Haardt | Subspace methods and exploitation of special array structures[END_REF]], the superdiagonal terms in the block V 22 are given, for q ∈ {1, . . . , Q -1}, by [V 22 ] q,q+1 = sign(h q h q+1 ) u S (τ , h q , h q+1 ) Υ q (τ , s q , s q+1 , h q , h q+1 ) ρ |hq| q (ε hq (s q )) ρ |h q+1 | q+1 (ε h q+1 (s q+1 )) (6.70) in which we have used the following denitions i) for q ∈ {1, . . . , Q -1} and (h q , h q+1 ) ∈ Z 2 , u S (τ , h q , h q+1 )

     (τ -|hq|)(τ -|h q+1 |) τ 3 ifq < Qand max(|h q |, |h q+1 |) < τ τ -|h Q | τ 2 ifq = Qand|h Q | < τ 0 if max(|h q |, |h q+1 |) τ , (6.71)
ii) for q ∈ {1, . . . , Q -1} and (s, s ) ∈ ]0 , 1[ 2 , r q (s, s ) Ω f s (x ; η q-1 ) f s (x ; η q ) f 1-s-s (x ; η q+1 ) dx, (6.72) 130 CHAPTER 6. APPENDICES iii) for q ∈ {1, . . . , Q -1} and (s, s ) ∈ ]0 , 1[ 2 , R q (s, s ) ρ q (s) ρ q+1 (s ) r q (s, s -s) , (6.73) and iv) for q ∈ {1, . . . , Q -1}, (h q , h q+1 ) ∈ Z 2 and (s q , s q+1 ) ∈ ]0 , 1 ε hq (sq), ε h q+1 (s q+1 )

1 -Rq ε hq (sq), ε h q+1 (s q+1 ) (6.74) in which, for x ∈ R, (x) + max(x, 0).

Block V 12

In this subsection, we are interested in the elements of the block V 12 . For p ∈ {1, . . . , L(Q + 1)} and q ∈ {1, . . . , Q}, and setting p = Lq + with q ∈ {0, . . . , Q} and ∈ {1, . . . , L}, (4.52) can be rewritten as [V 12 ] p,q = E x,t;η ∂ ln f (x, t; η) ∂η q, η f sq (x, t + h q ; η ) f sq (x, t ; η ) -f 1-sq (x, th q ; η ) f 1-sq (x, t ; η ) . where, for q ∈ {1, . . . , Q}, the L × 1 vectors v q and w q have components that can be written, for ∈ {1, . . . , L}, on the one hand, as v q, = [V 12 ] L(q-1)+ ,q = -h q u D (τ , h q ) ρ |hq|-1 q (ε hq (s q )) ϕ η q-1, ,q (ε hq (s q )) (6.77) and, on the other hand, w q, = [V 12 ] Lq+ ,q = h q u D (τ , h q ) ρ |hq|-1 q (ε hq (s q )) ϕ η q, ,q (ε hq (s q )) (6.78) where, for j ∈ {1, . . . , Q},  ∈ {j -1, j}, k ∈ {1, . . . , L}, and s ∈ ]0 , 1[, ϕ η ,k ,j (s) is dened by

ϕ η ,k ,j (s) Ω ∂ ln f (x ; η )
∂η ,k η f s (x ; η j-1 ) f 1-s (x ; η j ) dx. To conclude, equations (6.62), (6.64), (6.65), (6.69), (6.70), (6.76), (6.77) and (6.78) provide all the expressions necessary to determine the elements of the matrix W (H, s) CV -1 C T in (4.53). It is worth noticing that, due to the structure of the matrices V 11 , V 12 and V 22 , the inversion of V should not be particularly dicult from a computational point of view.

In the next section, we give more explicit expressions of these elements for widely encountered distributions in signal processing applications, namely, Gaussian [START_REF] Basseville | Detection of Abrupt Changes, Theory and Application[END_REF]235238] and Poisson [239 243] distributions. For each of these cases, we give closed-form expressions for i) F (η q ), dened in (6.63) and which leads to (6.62); ii) ρ q (s), dened in (6.67), which directly leads to (6.65) and (6.69), and partly to (6.70), (6.77) and (6.78); iii) R q (s, s ), dened in (6.72), which leads to (6.70); and nally for iv) ϕ η q, ,q (s), dened in (6.79), which leads to (6.77) and (6.78).

Gaussian case

The Gaussian distribution has perhaps the widest range of applications (see for instance [START_REF] Basseville | Detection of Abrupt Changes, Theory and Application[END_REF][START_REF] Chen | Parametric Statistical Change Point Analysis[END_REF][START_REF] Fotopoulos | Exact asymptotic distribution of change-point MLE for change in the mean of Gaussian sequences[END_REF], and references therein for an overview of the potential applications). In such cases, the model (4.37) is relevant: we can consider that the signal is piecewise Gaussian, i.i.d., that is, for q ∈ {0, . . . , Q} and n ∈ {t q + 1, . . . , t q+1 } (with t 0 0 and t Q+1 N ), we have x n ∼ N (µ q , σ 2 q ).

In other words, both mean and variance are likely to change from one segment to another, which means that L = 2 and the parameter vector η q [µ q , σ 2 q ] T includes the mean µ q and the variance σ 2 q of the signal on the (q + 1)-th segment (q ∈ {0, . . . , Q}).

Straightforward computations lead to the following explicit expressions

• for F (η q ), q ∈ {0, . . . , Q},

F (η q ) = diag 1 σ 2 q , 1 2(σ 2 q ) 2 , (6.80)
• for ρ q (s), q ∈ {1, . . . , Q}, s ∈ ]0 , 1[, ρ q (s) = snr v q-1,q s s snr v q-1,q + 1 -s exp -s(1 -s) snr m q-1,q 2(s snr v q-1,q + 1 -s) ,

in which we have dened the following two signal-to-noise ratios (SNR), for (q, q ) ∈ {0, . . . , Q} 2 : snr m q,q (µ q -µ q ) 2 σ 2 q and snr v q,q σ 2 q σ 2 q , (6.82) which characterize the amount of change in the mean and variance, respectively.

• for R q (s, s ), q ∈ {1, . . . , Q -1} and (s, s ) ∈ ]0 , 1[ 2 , after tedious computations Rq(s, s ) = s snr v q-1,q+1 + (s -s) snr v q,q+1 + 1 -s (s snr v q-1,q + 1 -s)(s snr v q,q+1 + 1 -s ) × exp(-s(1 -s ) 2(s snr v q-1,q+1 + (s -s) snr v q,q+1 + 1 -s ) ( s snr v q-1,q+1 + 1 -s s snr v q-1,q + 1 -s snr m q-1,q + s snr v q-1,q+1 + 1 -s s snr v q,q+1 + 1 -s snr m q,q+1 -snr m q-1,q+1 )), (6.83)

• for ϕ µ q ,q (s), q ∈ {0, . . . , Q}, ∈ {1, . . . , L}, q ∈ {q -1, q}, s ∈ ]0 , 1[, we obtain: ϕ µ q ,q (s) = ρ q (s) σ 2 q sµ q-1 snr v q-1,q + (1 -s)µ q s snr v q-1,q + 1 -s -µ q (6.84) and ϕ σ 2 q ,q (s) = ρ q (s) 2σ 2 q sµ q-1 snr v q-1,q + (1 -s)µ q 2 σ 2 q (s snr v q-1,q + 1 -s) 2 + σ 2 q -2µ q(sµ q-1 snr v q-1,q + (1 -s)µ q ) σ 2 q (s snr v q-1,q + 1 -s) + µ 2 q σ 2 q -1 .

Poisson case

The case of Poisson observations is also of interest in a number of signal processing applications, for instance for the segmentation of astronomical time series, possibly multivariate [START_REF] Scargle | Studies in astronomical time series analysis: V. Bayesian blocks, a new method to analyze structure in photon counting data[END_REF][START_REF] Dobigeon | Joint segmentation of multivariate astronomical time series: Bayesian sampling with a hierarchical model[END_REF].

Let us assume that the observations are modeled according to (4.37), where the distribution on each segment is Poisson, i.e., for q ∈ {0, . . . , Q} and n ∈ {t q +1, . . . , t q+1 } (with t 0 0 and t Q+1 N ), we have x n ∼ P(λ q ), or equivalently, f (x n ; λ q ) = Pr(X = x n ; λ q ) = exp{-λ q }λ xn q /x n !. In this case, we have η q λ q in the (q + 1)-th segment, which is a scalar parameter(L = 1).

Similarly to the case of Gaussian observations and after some computations, we obtain the following explicit expressions:

• For F (λ q ), q ∈ {0, . . . , Q}, F (λ q ) = 1 λ q , (6.86)

• for ρ q (s), q ∈ {1, . . . , Q}, s ∈ ]0 , 1[, ρ q (s) = exp{-sλ q-1 -(1 -s)λ q + λ s q-1 λ 1-s q } (6.87)

• for R q (s, s ), q ∈ {1, . . . , Q -1} and (s, s ) ∈ ]0 , 1[ 2 , tedious computations lead to R q (s, s ) = exp -λ q 1 -λ q-1 λ q s λ q+1 λ q 1-s + λ q-1 λ q s λ q+1 λ q 1-s (6.88)

• for ϕ λ q ,q (s), q ∈ {0, . . . , Q}, ∈ {1, . . . , L}, q ∈ {q -1, q}, s ∈ ]0 , 1[, tedious computations yield ϕ λ q ,q (s) = ρ q (s) λ s q-1 λ 1-s q λ q -1 .

(6.89)

We nally obtain the HCRWWB for a Poisson distributed signal that includes Q changepoints by plugging these last expressions into the equations (6.62), (6.64), (6.65), (6.69), (6.70), (6.76), (6.77) and (6.78) from Section 6.5.1, and by applying the procedure described in Appendix 4.2.2.3, which leads to the tightest bound.

  Aliated member of the Sensor Array and Multichannel Technical Committee in the IEEE Signal Processing Society (IEEE SPS SAM TC), since Mars 2013. • Outstanding reviewer (top 10th percentile of reviewers) : Elsevier -Signal Processing 2016 & Elsevier -Digital Signal Processing 2016

  ≈ 7/year (IEEE Transactions on Signal Processing, Elsevier Signal Processing, Signal Processing Letters, IEEE Transactions on Information Theory, Digital Signal Processing). International Conferences: ≈ 20/year (IEEE conferences and EUSIPCO).

Chapter 1 Introduction

 1 In this chapter, I present briey the main contributions that I obtained during the past ve years, i.e., since I held the position of Assistant Professor at IUT de Ville d'Avray. This work concerns the study of some recent and ubiquitous problems related to array signal processing in non standard conditions: mismatched models, presence of outliers, non-Gaussian noise, low signal to noise ratio, low number of samples, etc. Before developing the aforementioned issues, I rst, present in the next subsection a summary of my Ph.D. thesis in order to highlight the dierence and evolution from my Ph.D. subject/work to my current research. 1.1 Summary of my Ph.D. thesis During my Ph.D. thesis I focused on the study of performance analysis in array signal process-

  These studies allowed me to obtain interesting results and a comfortable level of accomplishment with 4 IEEE Transactions on Signal Processing Journal, 2 Elsevier Signal Processing Journal and 1 EURASIP JSAP during my Ph.D. thesis [J1-J7]. Nevertheless, I felt the need to move forward and to explore new leads in the area of array signal processing. In this spirit, I decided to work on some recent and ubiquitous problems and emerging subjects related to array signal processing in non standard conditions, namely, mismatched models, presence of outliers, non-Gaussian noise, low signal to noise ratio, low number of samples, etc. 1.2 Overview of my current research activity Since September 1st, 2013, I am holding an assistant professor position at IUT de Ville d'Avray, and I am aliated to the Laboratoire Energétique Mécanique Electromagnétisme (LEME EA-4416). At the time of my recruitment, I was the sole researcher from CNU-61 in the LEME. While this can be seen as a disadvantage, I took it as an advantage by freely proposing new subjects, external collaborations, partners to my laboratory and by giving a signal processing touch to the LEME. This lead us to recruit, in our laboratory, permanent and aliated researchers as well Ph.D. students in the area of signal processing. Nowadays, I am the head of the Equipe Signal group of the LEME laboratory. Regarding, my research activity, while maintaining a parallel research activity for new problems related to performance analysis (mostly in non-Gaussian scenario), I guided my research interest towards robust and/or scalable parametric estimation procedure with applications in the eld of array signal processing.

20 CHAPTER 1 .

 201 INTRODUCTION in the Ph.D. thesis of V. Ollier (ANR MAGELLAN Project), the ARPE internship of M. Brossard (Grant from ENS Paris-Sacaly jointly with TU-Darmstadt as host university), the ANACONDA project (an I-code grant) and the ON-FIRE project (GDR-ISIS Jeunes Chercheurs grant for which I am the principal investigator) [L2,JS1,J18,J16,IC33,IC28,IC27]. These projects have been giving me the opportunity to collaborate with many researchers, to cite few; Prof. P. Larzabal, Prof. A. Ferrari, Prof. M. Pesavento, Prof. S. Wijnholds, Dr. R. Boyer, Dr. A. Breloy, Dr. R. Flamary, Dr. F. Iutzeler.

  In the Ph.D. work of L. Bacharach, and with collaboration with Dr. A. Renaux, Prof. E. Chaumette and Prof. J.-Y. Tourneret, we are focusing on the performance analysis of change-point estimation. It is worth mentioning that change-point estimation has received 22 CHAPTER 1. INTRODUCTION

1 , 1 evec R - 1 eH

 111 in which Y e = κ 1 W -1 e +κ 2 vec R -, W e = R T e ⊗ R eand ∂r(µ) ∂µ refers to the Jacobian matrix of r(µ).

Figure 2 . 1 :

 21 Figure 2.1: MSE of SESAME, COCA and RCOMET vs. number of snapshots

For a set

  of K independent and identically distributed (i.i.d.) zero-mean observations following a CG distribution, we have the representation S d = DT, where, the columns of D ∈ C M ×K are the i.i.d. Gaussian vectors of covariance matrix Σ, and T = diag( √ τ

Fig. 2 .

 2 Fig. 2.8 displays the AFE of the dierent estimators in a critically small sample support (K = P ) for varying SNR. Fig. 2.5 also displays the same results for a number of samples K = 4P = N and Fig. 2.6 illustrates the impact of the concentration parameter κ on the AFE.

Figure 2 . 4 :

 24 Figure 2.4: Fraction of energy of U in R(U) versus SNR, K = P = 5, ν = 0.1, κ = 100, CG distribution prior for S and CIB prior for the subspace U.

Figure 2 . 5 :

 25 Figure 2.5: Fraction of energy of U in R(U) versus SNR, K = N = 20, ν = 0.1, κ = 100, CG distribution prior for S and CIB prior for the subspace U.

Figure 2 . 6 :

 26 Figure 2.6: Fraction of energy of U in R(U) versus κ, K = N , ν = 0.1, SNR= 5dB, CG distribution prior for S and CIB prior for the subspace U.

Figure 2 . 7 :

 27 Figure 2.7: Output of the ANMF using, respectively, the SCM, the SFPE, the MMSD and the Prior Only as subspace estimation process in the ANMF detector.

  the k th target, respectively; the transmit and receive steering vectors are dened as

(2. 70 )

 70 This leads to the IMAPE that diers from the IMLE in the estimation of Σ, that is adjusted by 2.3. ROBUST PARAMETERIZED MEAN ESTIMATION WITHOUT SECONDARY DATA 53 the estimation of the shape and scale parameter, as

2. 3 .

 3 ROBUST PARAMETERIZED MEAN ESTIMATION WITHOUT SECONDARY DATA 55 the CRB.

Figure 2 . 8 :

 28 Figure 2.8: MSE vs. L under K-distributed clutter, SCR = 15dB.

Figure 2 Figure 2 . 10 :

 2210 Figure 2.9: MSE vs. SCR under Kdistributed clutter

Figure 2

 2 Figure 2.11: MSE vs. SCR under tdistributed clutter.

  synthesis [96100]. These aspects are intertwined and must be carried out to take advantage of the new radio interferometers. As an example, lack of calibration has dramatic eects in the image reconstruction by causing severe distortions. In this Chapter, we focus on calibration, which involves the estimation of all unknown perturbation eects, as it is a cornerstone of the imaging step [101103]. Array calibration aspects have been tackled for a few decades in the array processing community leading to a variety of calibration algorithms [104106]. Such algorithms can be classied into two dierent approaches depending on the presence [107110], or the absence [111116], of

1 .

 1 First, (i) suboptimality due to the consideration of only one wavelength at a time; (ii) the assumption of a centralized processor, i.e., a single compute agent simultaneous accessing all data; and (iii) the ineciency with respect to the DoA estimation performances in 60 CHAPTER 3. ROBUST AND SCALABLE PARAMETRIC CALIBRATION WITH APPLICATION TO RADIO ASTRONOMICAL ARRAYS the severe radio astronomical contexts (which is generally treated using the MUSIC-based algorithm).

[ 61 3. 1

 611 SJ1,J18,J16,ICS1,IC33,IC28,IC27,IC24].3.1. SPARSE AND PARALLEL MULTI-WAVELENGTH CALIBRATION ALGORITHM Sparse and parallel multi-wavelength calibration algorithmBefore presenting the proposed algorithm, let us summaries its main goal. We aim, in this section, at designing an iterative algorithm, namely the Parallel Multi-wavelength Calibration Algorithm (PMCA), that focuses on the calibration of a sensor array based radio interferometer, involving its individual antennas and propagation disturbances. In addition, we assume that the sensor array has an arbitrary geometry, identical elements and is simultaneously excited by inaccurately known calibration sources. The proposed PMCA overcomes the aforementioned limitations, by: i) reformulating the parametric model in the multi-wavelength scenario in order to exploit wavelength diversity; ii) relying on distributed and consensus algorithms; and iii) adapting the sparse reconstruction methods to the calibration of radio interferometers. From the parallel calibration perspective, the PMCA successively estimates the undirectional antenna gains along with the directional and noise parameters for multiple subbands, where we enforce the coherence over the wavelength of the estimates based on physical and astronomical phenomena[135138].

3. 1 . 4

 14 Proposed parallel multi-wavelength calibration algorithm 3.1.4.1 Overview of the proposed parallel multi-wavelength calibration algorithm

.14) 3 . 1 .

 31 SPARSE AND PARALLEL MULTI-WAVELENGTH CALIBRATION ALGORITHM 65 Consequently, we dene the cost function for the entire network as κ(p) = λ∈Λ κ λ (p λ ).

  1.1 of the PMCA. As shown in the box Algorithm.1.1, this optimization step is performed w.r.t. the direction independent (DI) gain parameters {g λ } λ∈Λ , while the remaining parameters {m λ , d λ,1 , . . . , d λ,Q , σ n λ } λ∈Λ of p are xed. During this step, each agent calibrates the sensor gains g λ using the data available locally and Algorithm.1.1a.

} 3

 3 λ∈Λ , the fusion center can update α with a closed-form expression, and then send its estimated value, α[t+1] , to all agents in the network. Specically, we obtain after some calculusα[t+1] = vec IU [t+1] T T = vec U [t+1] T T .Direction dependent parameter and noise power estimationIn this subsection, we describe Algorithm.1.2 of the PMCA dedicated to the estimation of the direction dependent (DD) parameters and noise powers, i.e., {m λ , d λ,1 , . . . , d λ,Q , σ n λ } λ∈Λ for xed {g λ } λ∈Λ , which is based mainly on a sparse representation framework with a parallel implementation.

λ 2 2 s

 2 .t. m 0, σ n λ 0, ∀λ ∈ Λ mq 0 = 1 for q = 1, . . . , Q.

(3. 30 ) 68 CHAPTER 3 .

 30683 ROBUST AND SCALABLE PARAMETRIC CALIBRATION WITH APPLICATION TO RADIO ASTRONOMICAL ARRAYS where

; 23 23

 23 Deduce locally Dλ λ∈Λ from D [k]

Figure 3 . 2 :

 32 Figure 3.2: LOFAR's Initial Test Station antenna locations [155].

3. 1 .K = 2 , ρ = 5 • 10 3 PK = 3 , ρ = 5 •Figure 3 . 3 :Figure 3 . 4 :

 1253353334 Figure 3.3: RMSE on {g λ } λ as function of wavelength and compared to the CBRs. The edge wavelengths have a higher error, particularly for K = 2 and K = 4, due to our erroneous choice of interpolating polynomial order compared to the real/true polynomial order K K = 3.

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: RMSE on the DI gains and antenna noise powers as function of number of samples N , and compared to their corresponding multi-constrained-CRB.

3. 1 . 73 Cyg

 173 SPARSE AND PARALLEL MULTI-WAVELENGTH CALIBRATION ALGORITHM

Figure 3 . 7 :1

 37 Figure 3.7: image of the uncalibrated covariance matrix at 49.4 MHz. The presence of diuse emission is noticeable.

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: Image results of the calibrated covariance matrix at 49.4 MHz from the proposed PMCA (left), the calibrated solution of the existing method used for LOFAR [127] (right). The diuse emission appears slightly brighter in the image calibrated using PMCA.

( 3 . 49 )

 349 Let us note w = [w T 1 , . . . , w T D ] T the complete data vector, whose covariance matrix, denoted as Ξ, has the following form Ξ = bdiag{β 1 Ψ, . . . , β D Ψ}.

Algorithm 2 : 2 while stop criterion unreached do 3 3 E 4 M

 2234 Relaxed concentrated ML based calibration algorithm input : D, M , B, C i , β i , x output : θ initialize: Ω ← Ω init , τ ← τ init , θ ← θ init 1 while stop criterion unreached do -step: ŵi obtained from (3.51), i ∈ {1, . . . , D} 4 -step: θi obtained as follows, i ∈ {1, . . . , D} 5 while stop criterion unreached do 6 θi,p obtained from (3.53), p ∈ {1, . . . , M } interest ε 3DC in a sequential manner. To do so, we use an iterative estimation procedure by optimizing a cost function w.r.t. one parameter while xing the others. 1) Estimation of g p : The diagonal elements of the gain matrix are given by ĝp = arg min gp κ(g p ) (3.54) where κ(g

Figure 3 .

 3 Figure 3.13, with and without the calibration source, before adding any perturbation. The inte-

Figure 3 . 10 :Figure 3 . 11 :Figure 3 . 12 :Figure 3 . 13 :Figure 3 . 14 :Figure 3 . 15 :

 310311312313314315 Figure 3.10: MSE of η [f 1 ] 1 vs. SNR, for D = 2 bright signal sources, M = 8 antennas and 4 weak

(4. 14 ) 4 . 1 . 93 4. 1 . 2

 14419312 MISSPECIFIED CRAMÉR-RAO BOUND AND ITS SLEPIAN-BANGS-TYPE FORMULAS FOR CES DISTRIBUTIONS Slepian-Bangs formulas under misspecied CES models

4. 1 .

 1 MISSPECIFIED CRAMÉR-RAO BOUND AND ITS SLEPIAN-BANGS-TYPE FORMULAS FOR CES DISTRIBUTIONS 95

99 4. 2

 992 IN THE PRESENCE OF UNKNOWN CONTINUOUS AND DISCRETE PARAMETERS: APPLICATION TO CHANGE-POINTS ESTIMATION Performance analysis in the presence of unknown continuous and discrete parameters: Application to change-points estimation Non-stationary signals are often encountered in signal processing applications. Abrupt changes are a common cause of such non-stationarity. The latter occurs when statistical properties, of random observations, change abruptly, i.e., quickly w.r.t. the sampling period. In this section, we focus on the o-line change-point estimation problem and aims at determining performance measure for this problem by deriving lower bounds on the MSE. In this context, the derivation of the CRB, which is the most widely used lower bound, is unfeasible. This is due to the discrete nature of the change-points, thus the required regularity conditions under which the CRB is a lower bound are not met, and the CRB of these parameters cannot be derived. To overcome this issue, there exist other lower bounds on the MSE which do not require the dierentiability of the log-likelihood. One can cite deterministic lower bounds including the family of Barankin bounds [200205], and Bayesian lower bounds such as Weiss-Weinstein bounds [206]. Some lower bounds have been derived for the change-point estimation problem. Specically, regarding deterministic lower bounds, Ferrari and Tourneret derived the Chapman-Robbins (that is a Barankin-type) bound for a single change-point estimation [207]. Almost one decade later, this work was extended by La Rosa et al. to the case of multiple change-points [208]. It turns out that these deterministic bounds i) only give a coarse insight on the change-point estimation behavior and ii) considered the change-points to be the unique unknown parameters, meaning

(4. 52 )

 52 This block corresponds to the cross-terms between the Cramér-Rao and the Weiss-Weinstein lower bounds.

  s * = 0.5 HCRWWB{µ2}, s * = 0.5 HCRWWB{µ3}, s * = 0.5 HCRWWB{µ4}, s * = 0.5

Figure 4 . 1 :

 41 Figure 4.1: Empirical GMSEs and HCR-WWBs for the mean estimates on each segment, for Q = 3 changes in the mean of Gaussian observations (N = 100).

  s * = 0.5 HCRWWB{t2}, s * = 0.5 HCRWWB{t3}, s * = 0.5

Figure 4 . 2 :

 42 Figure 4.2: Empirical GRMSEs and HCR-WWBs for the change-point locations estimates, for Q = 3 changes in the mean of Gaussian observations (N = 100).

4. 2 . 4 . 3 5 Figure 4 . 3 :

 243543 Figure 4.3: Empirical GMSEs and HCR-WWBs for the variance estimates on each segment, for Q = 3 changes in the variance of Gaussian observations (N = 100).

  s * = 0.5 HCRWWB{t2}, s * = 0.5 HCRWWB{t3}, s * = 0.5

Figure 4 . 4 :

 44 Figure 4.4: Empirical GRMSEs and HCR-WWBs for the change-point location estimates, for Q = 3 changes in the variance of Gaussian observations (N = 100).

2 MLMAP and σ 2 3 MLMAP, 1 MLMAP. For σ 2 4 MLMAP,

 2314 4.3 and 4.4. Note that results for the mean estimates are not shown here but are very similar to those obtained for variance estimates in Fig.4.3.The same remarks as those made for mean changes (see previous subsection) are valid in the present case, with the following slight dierences: the non-information zone, regarding the variance estimates σ 2 ranges from 0 dB to 25 dB, while it ranges from 0 dB to 15 dB regarding σ 2 one cannot distinguish such non-information zone. With regard to the change-point location estimates, the gap between the GRMSEs and the HCRWWBs becomes lower than 0.1 samples for snr v ≥ 25 dB.

108 CHAPTER 4 . 5 Figure 4 . 5 :

 1084545 Figure 4.5: Empirical GMSEs and HCR-WWBs for the mean rate estimates on each segment, for Q = 3 changes in the mean rate of Poisson observations (N = 100).

  s * = 0.5 HCRWWB{t2}, s * = 0.5 HCRWWB{t3}, s * = 0.5

Figure 4 . 6 :

 46 Figure 4.6: Empirical GRMSEs and HCR-WWBs for the change-point locations estimates, for Q = 3 changes in the mean rate of Poisson observations (N = 100).

110 CHAPTER 5 .

 1105 These two main steps are, respectively, based on Alternating Direction of Multiple Multipliers and Distributed Iterative Hard Thresholding procedures. This leads to a statistically ecient and robust scheme to unmodeled sources as shown by numerical simulations and 109 CONCLUSION & PERSPECTIVES by application to actual LOFAR data. PMCA oers a framework that is easily adaptable to variations on the scenarios presented in Section 3.1.1, while maximizing the use of a priori available knowledge of the physics of the instrument and the measurement process.

1. 1 -

 1 Parametric robust adaptive signal processing, with application to MIMO systems The main goal is to propose/study new estimations procedures and their performances in joint estimation of targets and clutter parameters in the context where signals are embedded in a correlated non-Gaussian noise. Very few works have been done yet in this eld as pointed out

Figure 5 . 1 :

 51 Figure 5.1: Detection using data from STAP DGA-MI: the result using (top) the generalized robust shrinkage estimator, (bottom) our robust subspace clustering algorithm [ICS3].

  6.14) 6.4. SB-FORMULA FOR MISSPECIFIED MODELS: PROOFS AND RELATIONSHIP TO PREVIOUS RESULTS 121

  [START_REF] Greenewald | Robust sar stap via kronecker decomposition[END_REF]:∂D (p X f θ ) ∂θ i θ= θ = -E p ∂ln f X (x; θ) ∂θ i θ= θ = 0, i = 1, . . . , d.

( 6 2 0

 62 .40) and consequently, the Stochastic Representation Theorem allows us to write the following equality:x l -γ 0 l = x l -µ l = d RTu l = d RΠ 1/u l .(6.41)

  obtained by rstly taking the derivative of the terms in (4.22) and then by substituting in the obtained expression the Stochastic Representation of the dierence vector x l -γ 0 l in (6.41), we have that: -Q l u H l P 0 i u l . (6.44) 6.4. SB-FORMULA FOR MISSPECIFIED MODELS: PROOFS AND RELATIONSHIP TO PREVIOUS RESULTS 125

E

  p {u H l P i u l } = E p {tr(P i u l u H l )} = N -1 tr(P i )

  [START_REF] Ginolhac | Approximate distribution of the low-rank adaptive normalized matched lter test statistic under the null hypothesis[END_REF], and remembering that E p u l u H l = (1/N )I (see Section 4.1.1), we get:

6. 5

 5 Closed form expressions of the hybrid Cramér-Rao-Weiss-Weinstein bound for the change-point problem In this Appendix, we give the main steps leading to the expressions of the HCRWWB for the change-point problem 6.5.1 General case Let us derive the expression of the bound for signals which include a given number Q of abrupt changes, i.e., those modeled by equation (4.37)

6. 5 .

 5 CLOSED FORM EXPRESSIONS OF THE HYBRID CRAMÉR-RAO-WEISS-WEINSTEIN BOUND FOR THE CHANGE-POINT PROBLEM 129 6.5.1.2 Blocks V 22 and C 22

[ 2 ,

 2 Υq(τ , sq, s q+1 , hq, h q+1 )2(τ -|hq| -|h q+1 |) + + (τ -|hq| -|h q+1 | + 1) + -(τ -max(|hq|, |h q+1 |)) + -1 -R 1-min(|hq |,|h q+1 |) q

(6. 75 )

 75 After some tedious calculus, the matrix block V 12 has the form

6. 5 . 131 6. 5 . 2

 513152 CLOSED FORM EXPRESSIONS OF THE HYBRID CRAMÉR-RAO-WEISS-WEINSTEIN BOUND FOR THE CHANGE-POINT PROBLEM Case of Gaussian and Poisson distributions

(6. 85 ) 132 CHAPTER

 85132 Using these expressions and plugging them into the appropriate equations from Appendix 6.5.1 lead to the HCRWWB for a Gaussian signal submitted to Q abrupt changes.

  Responsible of the students recruitment procedure for IUT GEII Ville d'Avray, approximatively 100 students/year (total number of applicant around 1300/year). Member of 7 Ph.D. defenses committee (Ph.D. defenses of A. Boudjellal Sept. 2015 at University of Orléans-France, X. Zhang Aug. 2016 at TU-Darmstadt-Germany, S. Chiarucci Dec. 2016 at University of Florence-Italy and A. Mennad Dec. 2017 at EMP-Algiers A. Bitar Jun. 2018 at SONDRA-France, V. Ollier Jul. 2018 at SATIE-France, L. Bacharcah

	Students recruitment
	Setting up of the timetable (room management, teachers, temporary sta, ...)
	Management of internship and apprenticeship research
	Organization of the defenses and preparation of the jury
	Budget management (around 75k euros/year)

of Paris Nanterre University since Sept 2014. EESC stands for Embedded Electronics and Communications Systems (15 lectures and approximatively 22 students/year). My missions are • • Member of the improvement council board of the Master of Science EESC (Embedded Electronics and Communications Systems) of Paris Nanterre University since Sept 2016. 0.11 Miscellaneous •

  of A into a vector. The operators ⊗, and • refer, respectively, to Kronecker, element-wise and Khatri-Rao (column-wise Kronecker) matrix product. E p {.} and E X/Y {.} denote respectively, the expectation under the pdf p and the conditional expectation. E {.} will be used when there is no ambiguity.

	In the following, vectors (respectively matrices) are denoted by boldface lowercase
	letters (respectively uppercase letters). The notation	d = indicates has the same distribution as".
	Convergence in distribution and in probability are, respectively, denoted by	d → and	P →. For a

matrix A, |A| and Tr (A) denote the determinant and the trace of A. A T (respectively A H ) stands for the transpose (respectively conjugate transpose) matrix. The vec-operator vec(A) stacks all columns

  Let N i.i.d. observations, y n ∼ CES m (0, R e , g) with N > m. We recall the xed-point equation satised by the unstructured ML estimator of R e :

	CHAPTER 2. ROBUST SCATTER MATRIX ESTIMATION AND SUBSPACE ESTIMATION WITH APPLICATION TO RADAR
	2.1.2.1 Step 1

  15) w.r.t. µ = [µ 1 , . . . , µ P ] T ∈ R P as already mentioned.

Theorem 1 The SESAME, µ, given by (2.15), is a consistent estimator of µ e . Likewise, R ( µ)

  2.1. EFFICIENT ESTIMATION OF COVARIANCE/SCATTER MATRICES WITH CONVEX STRUCTURE UNDER CES DISTRIBUTION 33The I-SESAME procedures Input N i.i.d. data, y n ∼ CES m (0, R e , g) with N > m Step 1 Compute R FP from y 1 , . . . , y N with(2.14)

  We generate 5000 sets of N independent m-dimensional t-distributed samples, y n ∼ Ct m,d (0, R e ) , n = 1, . . . , N with d = 5 degrees of freedom.

	28)
	2.1.5.2 Simulation for Toeplitz structure with t-distributed data
	We illustrate the SESAME performance for an Hermitian Toeplitz scatter matrix with t-distributed
	observations. For m = 4, the Toeplitz scatter matrix is generated from its rst row
	[1, -0.83 -0.20i, 0.78 + 0.37i, -0.66 -0.70i].

Table 2

 2

	recaps the averaging calculation time of the dierent algorithms. The proposed
	algorithms are globally time-ecient unlike Constrained Tyler and COCA estimators. The es-
	timation scheme with the joint-algorithm is slower than the one using the exact ML estimator,
	which makes sense since the degree of freedom of the t-distribution is also estimated. The COCA
	estimator suers from heavy computational cost, since the number of constraints grows linearly
	in N .

Table 2 :

 2 Algorithms running time.

	2.2. MINIMUM MEAN SQUARE DISTANCE ESTIMATION OF SUBSPACES IN PRESENCE OF NON GAUSSIAN SOURCES	37
	2.2 Minimum mean square distance estimation of subspaces in
	presence of non Gaussian sources
	In this section, we consider the subspace estimation which is an ubiquitous problem in signal pro-
	cessing, as it is often required to infer the low-dimensional subspace where information (usually
	the signals of interest) lies in. It is considered as a cornerstone of a plethora of applications and
	algorithms such as DoA estimation, signal detection and interference cancellation [32], to cite
	just a few. Furthermore, in modern array processing applications, this task reveals to be dicult
	since we frequently deal with non standard conditions, i.e., low sample support, low Signal to
	noise ratio (SNR), non-Gaussian observations, outliers in the sample set, etc.
	Most commonly, the signal subspace is estimated through the strongest eigenvectors of the
	singular value decomposition (SVD) of the SCM. This corresponds to the Maximum Likelihood
	Estimator (MLE) for the classical linear model with additive white Gaussian noise. The latter
	aords an accurate estimator for high SNR and/or for large number of samples. Nevertheless, it

  ∈C N ×P |U H U = I P the set of N × P unitary matrices, Y ∈ C N ×K the data matrix,

	CHAPTER 2. ROBUST SCATTER MATRIX ESTIMATION AND SUBSPACE ESTIMATION WITH APPLICATION TO RADAR
	2.2.1 Background theory
	Let us denote in the following, S + N the set of positive semi-denite matrices of dimension N × N , U N P = U
	to as Bingham-von-Mises-Fisher) distribution and propose a method to sample form
	this distribution. Finally, numerical simulations illustrates the performance of the proposed
	algorithms.

Table 1 :

 1 Φ, A) {Ap} = {φ p }A ∈ S + N where Φ = diag(φ) p CBL ∝ etr Re{C H U} + ΦU H AU Generalized Bingham CGB({Ap}) {Ap} ⊂ S + CL ∝ etr Re{C H U} Special cases of the CGBL distribution.

	2.2. MINIMUM MEAN SQUARE DISTANCE ESTIMATION OF SUBSPACES IN PRESENCE OF NON GAUSSIAN SOURCES	39
	Complex Distribution	Parameters	Probability density function
	Bingham Langevin CBL(C, N Bingham CB(Φ, A) A ∈ S + N Invariant Bingham CIB(κ, A) κ ∈ C and A ∈ S + N Langevin CL(C) C ∈ C M ×R	p CGB ∝ etr P p=1 u H p Apup p CB ∝ etr ΦU H AU p CSB ∝ etr κU H AU

p

  U|Y, τ t , λ t )) = UU H p(U|Y, τ t , λ t )dU

	2.2. MINIMUM MEAN SQUARE DISTANCE ESTIMATION OF SUBSPACES IN PRESENCE OF NON GAUSSIAN SOURCES	41
	with
	M(p((2.41)
	From (2.33) and (2.37), some manipulations allow to the posterior probability as:
	40)

Table 2 :

 2 A summary of the MMSD estimator in function of the prior of U.

  • , and of the second source are 45 • and 40 • . The coecients α 1 and α 2 are chosen to be 2 + 3j and 1 -0.5j, and the normalized Doppler frequencies f 1 f 2 are 0.3 and 0.8. There are L = 15 pulses per CPI, and each pulse contains T = 5 snapshots. For K-distributed clutter, we choose a = 2 and b = 10; and for t-distributed clutter, a = 1.1 and b = 2. The entries of the speckle covariance matrix Σ are generated by

	56	CHAPTER 2. ROBUST SCATTER MATRIX ESTIMATION AND SUBSPACE ESTIMATION WITH APPLICATION TO RADAR

1

The so-called ARL refers to the RL when angular parameters are considered as the only parameters of interest. and

  Figs. 2.8-2.11 investigates the performance of the proposed estimators, i.e., the IMMLE, IMLE, IMAPE, w.r.t. the classical CMLE, as well as with the derived CRB. In Figs. 2.8 and 2.9, the MSEs under K-distributed clutter are plotted, versus the pulse number L and the SCR, respectively; and in Figs. 2.10 and 2.11 the MSEs under t-distributed clutter, also versus the pulse number and the SCR, respectively. From Figs. 2.8-2.11, it becomes obvious that the conventional algorithm (CMLE) becomes poor when the clutter is a CG, and all the IMMLE, IMLE and the IMAPE algorithms lead to substantial superior performance. Figs. 2.8-2.11 also show that as few as two iterations for the IMLE and the IMAPE, and generally only one iteration for the IMMLE,

is enough for the respective estimators to have a satisfactory performance in terms of a resulting MSE appropriately close to the CRB, in asymptotic L and SCR cases. Furthermore, one can see from Figs. 2.8-2.11 that the performance of the IMMLE, which is a theoretically optimum estimator, is superior to that of the IMLE and the IMAPE, especially in severe scenarios. The reader is referred to Table

.

2.5 for a concise comparison between the IMMLE, IMAPE and IMLE.

Table 2 . 2

 22

	Repeat Step 1 and Step 2 until convergence
	Table 2.1: Summarization of the proposed algorithm

: Comparison between IMLE, IMAPE and IMMLE

  is to CHAPTER 3. ROBUST AND SCALABLE PARAMETRIC CALIBRATION WITH APPLICATION TO RADIO ASTRONOMICAL ARRAYS consider the problem as a consensus optimization problem for which the augmented Lagrangian reads

  3.1. SPARSE AND PARALLEL MULTI-WAVELENGTH CALIBRATION ALGORITHM 69 Algorithm 1.1: estimation of {g λ } λ∈Λ Input: Rλ λ∈Λ , p [i-1] , η g ;

  Deduce locally ĝλ λ∈Λ from α [t] using the constraints in (3.16); Output: ĝλ λ∈Λ ; Algorithm 1.1a: local estimation of {g λ } λ∈Λz

	2 repeat					λ	1]	λ∈Λ	;
	4 4	t = t + 1 ;
	6 6	Estimate locally g	[t] λ λ∈Λ	with Algorithm.1.1a ;
	8 8	Estimate α [t] at the fusion center with (3.22);
	10 10	Update locally the Lagrange multipliers y	[t] λ λ∈Λ	with (3.20);
	11 until λ∈Λ g [t-1] λ	-g	[t] λ 2	≤ λ∈Λ g	[t] λ 2	η g ;
	13 13

  CHAPTER 3. ROBUST AND SCALABLE PARAMETRIC CALIBRATION WITH APPLICATION TO RADIO ASTRONOMICAL ARRAYSestimator for the sequential updates of each block of unknown parameters, namely, τ , Ω and θ,

	are given by
	τ pq =	1 4	a H pq (θ)Ω -1 a pq (θ).	(3.43)
	Ω =	4 B pq	a pq (θ)a H pq (θ) a pq (θ) pq (θ) Ω a H -1	(3.44)
	followed by the normalization step
			Ω =	Ω tr Ω	.	(3.45)
	and nally, for a given Ω and τ , estimating θ is equivalent to the following minimization problem
	θ = arg min θ	pq	1 τ pq	a H

pq (θ)Ω -1 a pq (θ) .

  .58)2) Estimation of α i : We rst need to estimate ϕ i,p (we recall that ϕ i,p = ηi u p + ζ i v p ). || Ĵi,p -G p H i,p Z i,p (ϕ i,p )F i || 2F . Taking the derivative of κ(ϕ i,p ) w.r.t. ϕ i,p and setting the result to zero, we obtainTr j exp -j φi,p Ĵi,p F H i H H i,p G H p -j exp j φi,p G p H i,p F i

	3.2. ROBUST CALIBRATION OF RADIO INTERFEROMETERS IN NON-GAUSSIAN ENVIRONMENT 81
		ĴH i,p = 0	(3.60)
	which leads to
		Tr M i,p
	exp 2j φi,p =
		Tr M H i,p
	This is done as follows
	φi,p = arg min ϕ i,p	κ(ϕ i,p )	(3.59)

where κ(ϕ i,p ) =

Algorithm 3 :

 3 Case of structured Jones matrices input : D, M , B, C i , β i , x, Ĵi,p as output of Algorithm 3, i ∈ {1, . . . , D} and

	82	CHAPTER 3. ROBUST AND SCALABLE PARAMETRIC CALIBRATION WITH APPLICATION TO RADIO ASTRONOMICAL ARRAYS
	3.2.
			p ∈ {1, . . . , M }
			output : ε3DC
		initialize: ε3DC ← ε 3DC init 1 while stop criterion unreached do
		1 1	Obtain θi from (3.64), i ∈ {1, . . . , D}
		2 2	Obtain ĝp from (3.78), p ∈ {1, . . . , M }
		3 3	Obtain α i = [η i , ζ i ] T from (3.63), i ∈ {1, . . . , D}

  ,p as output of NSCA for i ∈ {1, . . . , D}, p ∈ {1, . . . , M } and f ∈ F

	MSCA 4: Multi-frequency structured calibration algorithm
	input : D, M , F , r [f ] p , H [f ] i,p ,	Ĵ[f]
	output : ˆ

scheme of the proposed multi-frequency structured calibration algorithm (MSCA) is described below. i

  88 CHAPTER 4. PERFORMANCE ANALYSIS WITH APPLICATIONS TO ARRAY PROCESSINGbest HCRWWB is proposed. This later is based on a semidenite programming formulation which consists in nding the unique minimum volume covering the set constituted by hyper-ellipsoid elements which are generated using the derived candidate HCRWWB matrices w.r.t. the so-called Loewner partial ordering.

	4.1. MISSPECIFIED CRAMÉR-RAO BOUND AND ITS SLEPIAN-BANGS-TYPE FORMULAS FOR CES DISTRIBUTIONS 89
	4.1 Misspecied Cramér-Rao bound and its Slepian-Bangs-type
	formulas for CES distributions
	This Chapter is based on the following articles [SJ5,J17,J14,ICS2,IC29,IC25,IC21,IC19].

  4.2. PERFORMANCE ANALYSIS IN THE PRESENCE OF UNKNOWN CONTINUOUS AND DISCRETE PARAMETERS: APPLICATION TO CHANGE-POINTS ESTIMATION 101 4.2.2 Proposed bound First, let us recall the general inequality leading to the proposed lower bound, namely the covariance inequality.

  4.2. PERFORMANCE ANALYSIS IN THE PRESENCE OF UNKNOWN CONTINUOUS AND DISCRETE PARAMETERS: APPLICATION TO CHANGE-POINTS ESTIMATION 103

  [START_REF] Meriaux | Ecient estimation of scatter matrix with convex structure under t-distribution[END_REF] for (p, p ) ∈ {1, . . . , P } 2 and θ ∈ Θ , the element [V 11 ] p,p in the matrix V is given by[V 11 ] p,p = E x,θr;θ d ∂ ln f (x, θ r ; θ d ) ∂θ d,pin which the last equality holds only if f (x, θ r ; θ d ) is twice dierentiable w.r.t. the vector θ. Note that this block corresponds to the so-called modied Cramér-Rao lower bound for the estimation error of an unbiased estimator of θ d[START_REF] Gini | The modied Cramér-Rao bound in vector parameter estimation[END_REF];

					48)
	in which
	i) (block V θ d	∂ ln f (x, θ r ; θ d ) ∂θ d,p	θ d	= -E x,θr;θ d	∂ 2 ln f (x, θ r ; θ d ) ∂θ d,p ∂θ d,p	θ d
					(4.49)

  104 CHAPTER 4. PERFORMANCE ANALYSIS WITH APPLICATIONS TO ARRAY PROCESSING Finally, for each value of H [h r1 , . . . , h rQ ] and s [s 1 , . . . , s Q ] T , we obtain a lower bound W (H, s) CV -1 C T on the mean square error. The proposed bound, that is the hybrid Cramér-RaoWeiss-Weinstein bound (HCRWWB), is the tightest of these lower bounds, i.e., it is obtained by maximizing W (H, s) w.r.t. H and s

	HCRWWB = sup	W (H, s).
	H∈H Q θr
	s∈]0,1[ Q

  Looking at(4.38), we can see that, as long as t belongs to its support T , the expression of the prior function π(t) does not explicitly depend on t. Consequently, the joint likelihood in (4.54) can be replaced with the classical likelihood f (x | t ; η). In addition, in the cases we study in the following sections (i.e., Gaussian and Poisson distributions), the maximization of the log likelihood w.r.t. η, for a given t, is not dicult and results in classical expressions of the empirical mean and/or variance for η(x ; t). Hence, it is possible to get rid of the dependence on η in (4.54), so that we obtain tMLMAP by using

		ln f (x, t ; η).	(4.54)
		η,t
	tMLMAP = arg max	ln f (x | t ; η(x ; t))	(4.55)
	t

  Furthermore, we have Σ = σ 1 W e +σ 2 vec (R e ) vec (R e ) H and the coecients dened by relations (2.10). To lighten the notations, we introduceV e = vec (R e ) vec (R e ) H U e = vec R -1 e vec R -1 (2κ 1 + mκ 2 ) U e Y e V e Y e = (κ 1 + mκ 2 ) 2 U e Y e ΣY e = κ 1 W -1 e + κ 2 U e = Y e

	Ω = ΣK and D h = where it is recalled that K is the commutation matrix. Moreover, we can easily check that ∂h ∂s = ∂r(µ) ∂µ H µ=µ e Y e KY T e = Y e K and K ∂r(µ) ∂µ * µ e = ∂r(µ) ∂µ µ e where A * refers to the conjugate matrix. Then, we obtain R∞ = 2 ∂r(µ) ∂µ H µ e YeΣYe ∂r(µ) ∂µ µ e + ∂r(µ) ∂µ H µ e YeΣYe ∂r(µ) ∂µ µ e . e H which satisfy the relations U e W e = vec R -1 e vec (R e ) H and V e W -1 e = vec (R e ) vec R -1 (6.6) e H . Moreover, we have vec (R e ) H vec R -1 e = m. After some calculus, we obtain: Y e W e Y e = κ 2 1 W -1 CHAPTER 6. APPENDICES Then, the equation (6.6) becomes e + κ 2 118 R∞ = 2 ∂r(µ) ∂µ H µ=µ e Ye ∂r(µ) ∂µ µ=µ e + ∂r(µ) µ=µ e ∂µ H Ye ∂r(µ)

  6.2. COMPLEX GENERALIZED BINGHAM LANGEVIN DISTRIBUTION SAMPLING 6.3 Notations and derivation of the block coordinate descent algorithm in the context of radio interferometers calibration We present here the steps to obtain (3.53) using the BCD in the context of radio interferometers calibration. Firstly, for sake of clarity, let us denote c i = [c i 1 , c i 2 , c i 3 , c i 4 ] T to refer to the four entries of the vectorization of source coherency matrix C i . Likewise, for the i-th source, we write J i,p (θ i,p ) = p i 1 p i 2 p i 3 p i 4

  and z l are mutually dependent and consequently the closed form evaluation of the expectation operator in (6.20) is not feasible. Similar considerations hold true for all the other terms involving the expectation operator w.r.t. the true pdf p X (x), that are

	122	CHAPTER 6. APPENDICES
	the random variable φ(G 0 l ) and the other random quantities Q 1/2 l
		6.21)

  It is immediate to verify that the matrix dened in(6.51)

		1)	-1 tr( Π-1	Πi )tr( Π-1	Πj )	(6.51)
	+ L	E Q {Q 2 φ2 (Q)} N (N + 1)	tr( Π-1	Πi	Π-1	Πj ),

where Π Π( θ) and Πi ∂Π(θ)

∂θ i | θ= θ.

  12 (equation (4.50)). Finally, we give the general expression of ξ(α, β, h a , h b ) from equation (4.51), which enables us to obtain the expressions of blocks C 22 and V 22 according to equations (4.46) and (4.50), respectively.6.5.1.1 Block V 11In this subsection, we give the expression of the elements of block V 11 for signals which include abrupt changes, i.e., those associated with model (4.37). Rewriting equation (4.49) with notations from Section 4.2.1, and setting the index changes p = Lq + and p = Lq + , with (q, q ) ∈ {0, . . . , Q} 2 and ( , ) ∈ {1, . . . , L} 2 , we have

[V 11 ] p,p = -E x,t;η ∂ 2 ln f (x, t ; η) ∂η q, ∂η q , η=η .

(6.61)

Let us recall that Z computational agents are disposed on a network, where the z-th agent, Az, accesses data for wavelengths λ ∈ Λz ⊂ Λ. In this perspective, we propose a parallel calibration scheme, in which the sensor gains corresponding to each wavelength λ are estimated locally and consensus is enforced among agents by imposing the aforementioned constraints C1-C5.

Educational and academic responsibilities

• Elected member of the CCD (Comité Consultatif des sections CNU 61/63) of Paris

CHAPTER 5. CONCLUSION & PERSPECTIVES

Again, one can notice that (5.5) is adapted to the regime 3 which can be extended to the regime 4 based on the above discussion.

• Finally, it is worth noting that in the next radio-interferometers generation, the sampling rate allows us to considerer multiple realisation of the SCM (Sample Covariance Matrix). A promising idea is to deal with the true distribution of the SCM which is in this case, Whisart distributed in a non-homogenous background. Taking such modeling, is theoretically interesting, nevertheless, it might lead to a non-scalable solution at hand. Consequently, eort have to be done in order to lighten the complexity cost.

Perspectives w.r.t. the study and analysis of algorithm based on misspecied models

Recent tools established on misspecied lower bound on the mean square error allow to predict more accurately the mean square error behavior than the classical lower bounds in presence of model errors. Specically, in [START_REF] Mennad | Slepian-bangs-type formulas and the related misspecied cramer-rao bounds for complex elliptically symmetric distributions[END_REF] we provided a new generalization of the classical SB (Slepian-Bang) formula to parametric estimation problems involving non-Gaussian, heavy-tailed, CES distributed data in the presence of model misspecication. While these bounds are helpful since model errors exist in practice due to system imperfections, still lot of perspectives are in order

• First, a natural extension of [START_REF] Smith | Statistical resolution limits and the complexied Cramér Rao bound[END_REF] to mismatched model is interesting. Specically, proposing and developing new complexied form of the MCRB (misspecied CRB) is of interest since, as its counterpart well-specied complexied CRB, it might provide a computationally ecient formulation from the perspective of derivations and analysis.

• Second, making the analogy w.r.t. to the well-specied bounds, we notice that there is a room of a plenty challenges. As an example, the constrained, the hybrid (i.e., the vector of parameters may contain deterministic and random parameters), the non-circular versions of the MCRB could be established. On the other hand, the extensions of our recent works [START_REF] Bacharach | Weiss-weinstein bound on multiple change-points estimation[END_REF][START_REF] Chaumette | A class of weiss-weinstein bounds and its relationship with the bobrovsky-mayer-wolf-zakai bounds[END_REF] is worth to be studied since it would allow us to propose a novel class of misspecied Weiss-Weinstein bounds (essentially free from regularity conditions on the probability density functions support and the natural of the parameter of interest which might include deterministic and random parameters).

• Finally, some practical examples can be considered in order to design and/or nd the targeted tradeo between the estimation accuracy and complexity of the model (and thus the complexity of its associated algorithm). As an example, this can provide a rigourous answer to some questions as the need of taking or not the structure of all/part parameters, at which frequency should we calibrate the system, at which level should we robustify/generalize the considered class of distribution to which belong the noise, etc. Chapter 6 Appendices 6.1 Some technical considerations w.r.t. the SESAME scheme 6.1.1 Practical implementation for holding the PSD constraint Through the mapping µ → R (µ), SESAME ensures that the estimate R S = R( µ) belongs to the desired convex subset S . In order to stress the reliance of the criterion J R ML , R (µ) dened in (2.15) on this mapping, we introduce the following notation:

Obviously, when it comes to scatter matrix estimation, the PSD constraint has to be complied.

Hence the considered subset is expressed S = S + M ∩ L , where S + M is the cone of PSD matrices, and L describes the particular structure of the CM (e.g., Toeplitz, banded, persymmetric, etc.).

However, in this case, the mapping µ → R (µ) does not usually possess an explicit expression. To overcome this issue, let us consider an auxiliary mapping µ → R L (µ) so that R L (µ) ∈ L , ∀µ.

The most common example is R L (µ) = P p=1 µ p B p , where {B p } form a basis of a linear convex set. We can then recast an equivalent SESAME as

which can be obtained using standard semi-denite program solvers. Nevertheless, this approach can be computationally demanding. Therefore, it is worth noting that if the constraint is omitted,

provides usually tractable expressions, and even closed form ones for linear structures. The major interest of this formulation is that (6.2) and (6.3) yield the same estimate if the PSD constraint is satised by R L ( μ) retrospectively. Moreover, for suciently large N , this relaxation provides a positive semi-denite matrix R L ( μ) with probability arbitrarily close to one, thanks to the consistency of SESAME (see Theorem 1.).

In Section 2.1.5, performance of these two proposed SESAME implementations are compared.

It is noted that imposing PSD constraint may be benecial for small sample support, but that performance of both implementations are generally equivalent.

We denote by v = g + ih where g and h are respectively the real and the imaginary parts of the complex random vector v. We assign w = {g 1 , h 1 , ..., g N , h N } is a 2N unit real random vector.

The generation of the vector v is equivalent to the generation of a real random vector w ∼ vBL(c, Ã) where c = {d 1 , f 1 , ..., d N , f N } is a 2N unit real random vector deduced from the unit vector c, in such a way, c = d + if , d and f are respectively the real and the imaginary parts of the complex random vector c and à is expressed in function of A, i.e., the block 2 × 2 matrices of the matrix à are dened as Ãij = δ ij cos(φ ij ) -sin(φ ij ) sin(φ ij ) cos(φ ij ) ∀i, j ∈ {1, ..., N }

where a ij = δ ij exp{iφ ij } are the elements of the matrix A. Note that à ∈ R 2N ×2N is a symmetric matrix. The following algorithm details the generation of the unit complex random vector v ∼ CVGBL(c, A) as described below:

Algorithm 5: The generation of the unit complex random vector v input : c, A output : v initialize: v (0) ← v init (unit norm vector) 1 while stop criterion unreached do 2 2 Compute the 2N real unit vector w from the vector v 3 3

Compute the 2N real unit vector c from the vector c 4 4

Compute the 2N × 2N real symmetric matrix à from the matrix A 5 5

Update the real unit random vector w = vBL(c, Ã) 6 6 Update the complex unit random vector v from w 7 end

The matrix complex generalized Bingham Langevin distribution

Inspiring from [START_REF] Ho | Simulation of the matrix bingham-von mises-sher distribution with applications to multivariate and relational data[END_REF], a Markov chain is generated which converges to CGBL(C, {A p }) in such a way the random unitary matrix X ∼ CGBL(C, {A p }). The procedure is detailed as follows:

Algorithm 6: The generation of the unitary matrix X input : C, {A p } output : X initialize: X (0) ← X init (a unitary matrix) 1 while stop criterion unreached do 2 for p ∈ {1, ..., P } in random order do 3 3 Compute the null space N of the matrix X [,-p] We denote by X [,-p] the matrix X after having removed the p-th column of X and the function.