Nizar Bouhlel
email: nizar.bouhlel@agrocampus-ouest.fr

David Rousseau
email: david.rousseau@univ-angers.fr

Makrides Andreas

A Generic Formula and Some Special Cases for the Kullback-Leibler Divergence between Central Multivariate Cauchy Distributions

Keywords: Multivariate Cauchy distribution (MCD), Kullback-Leibler divergence (KLD), multiple power series, Lauricella D-hypergeometric series

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

Introduction

Multivariate Cauchy distribution (MCD) belongs to the elliptical symmetric distributions [START_REF] Ollila | Complex Elliptically Symmetric Distributions: Survey, New Results and Applications[END_REF] and is a special case of the multivariate t-distribution [START_REF] Kotz | Multivariate T-Distributions and Their Applications[END_REF] and the multivariate stable distribution [START_REF] Press | Multivariate stable distributions[END_REF]. MCD has been recently used in several signal and image processing applications for which non-Gaussian models are needed. To name a few of them, in speckle denoizing, color image denoizing, watermarking, speech enhancement, among others. Sahu et al. in [START_REF] Sahu | Statistical modeling and Gaussianization procedure based de-speckling algorithm for retinal OCT images[END_REF] presented a denoizing method for speckle noise removal applied to a retinal optical coherence tomography (OCT) image. The method was based on the wavelet transform where the sub-bands coefficients were modeled using a Cauchy distribution. In [START_REF] Ranjani | Generalized SAR Despeckling Based on DTCWT Exploiting Interscale and Intrascale Dependences[END_REF], a dual tree complex wavelet transform (DTCWT)-based despeckling algorithm was proposed for synthetic aperture radar (SAR) images, where the DTCWT coefficients in each subband were modeled with a multivariate Cauchy distribution. In [START_REF] Sadreazami | Color image denoising using multivariate cauchy PDF in the contourlet domain[END_REF], a new color image denoizing method in the contourlet domain was suggested for reducing noise in images corrupted by Gaussian noise where the contourlet subband coefficients were described by the heavy-tailed MCD. Sadreazami et al. in [START_REF] Sadreazami | A Study of Multiplicative Watermark Detection in the Contourlet Domain Using Alpha-Stable Distributions[END_REF] put forward a novel multiplicative watermarking scheme in the contourlet domain where the watermark detector was based on the bivariate Cauchy distribution and designed to capture the across scale dependencies of the contourlet coefficients. Fontaine et al. in [START_REF] Fontaine | Cauchy Multichannel Speech Enhancement with a Deep Speech Prior[END_REF] proposed a semi-supervised multichannel speech enhancement system where both speech and noise follow the heavy-tailed multi-variate complex Cauchy distribution.

Kullback-Leibler divergence (KLD), also called relative entropy, is one of the most fundamental and important measures in information theory and statistics [START_REF] Cover | Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing)[END_REF][START_REF] Pardo | Statistical Inference Based on Divergence Measures[END_REF]. KLD was first introduced and studied by Kullback and Leibler [START_REF] Kullback | On Information and Sufficiency[END_REF] and Kullback [START_REF] Kullback | Information Theory and Statistics[END_REF] to measure the divergence between two probability mass functions in the case of discrete random variables and between two univariate or multivariate probability density functions in the case of continuous random variables. In the literature, numerous entropy and divergence measures have been suggested for measuring the similarity between probability distributions, such as Rényi [START_REF] Rényi | On Measures of Entropy and Information[END_REF] divergence, Sharma and Mittal [START_REF] Sharma | New non-additive measures of relative information[END_REF] divergence, Bhattacharyya [START_REF] Bhattacharyya | On a measure of divergence between two statistical populations defined by their probability distributions[END_REF][START_REF] Kailath | The Divergence and Bhattacharyya Distance Measures in Signal Selection[END_REF] divergence and Hellinger divergence measures [START_REF] Giet | A minimum Hellinger distance estimator for stochastic differential equations: An application to statistical inference for continuous time interest rate models[END_REF]. Other general divergence families have been also introduced and studied like the φ-divergence family of divergence measures defined simultaneously by Csiszár [START_REF] Csiszár | Eine informationstheoretische Ungleichung und ihre Anwendung auf den Beweis der Ergodizität von Markoffschen Ketten[END_REF] and Ali and Silvey [START_REF] Ali | A General Class of Coefficients of Divergence of One Distribution from Another[END_REF] where the KLD measure is a special case, the Bregman family divergence [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF], the R-divergences introduced by Burbea and Rao [START_REF] Burbea | On the convexity of some divergence measures based on entropy functions[END_REF][START_REF] Burbea | On the convexity of higher order Jensen differences based on entropy functions (Corresp.)[END_REF][START_REF] Burbea | Entropy differential metric, distance and divergence measures in probability spaces: A unified approach[END_REF], the statistical f -divergences [START_REF] Csiszar | Information-type measures of difference of probability distributions and indirect observation[END_REF][START_REF] Nielsen | On the chi square and higher-order chi distances for approximating f-divergences[END_REF] and recently the new family of a generalized divergence called the (h, φ)-divergence measures introduced and studied in Menéndez et al. [START_REF] Menéndez | Asymptotic behaviour and statistical applications of divergence measures in multinomial populations: A unified study[END_REF]. Readers are referred to [START_REF] Pardo | Statistical Inference Based on Divergence Measures[END_REF] for details about these divergence family measures.

KLD has a specific interpretation in coding theory [START_REF] Cover | Information theory and statistics[END_REF] and is therefore the most popular and widely used as well. Since information theoretic divergence and KLD in particular are ubiquitous in information sciences [START_REF] Mackay | Information Theory, Inference and Learning Algorithms[END_REF][START_REF] Ruiz | Information Theory in Computer Vision and Pattern Recognition[END_REF], it is therefore important to establish closed-form expressions of such divergence [START_REF] Nielsen | Statistical Divergences between Densities of Truncated Exponential Families with Nested Supports: Duo Bregman and Duo Jensen Divergences[END_REF]. An analytical expression of the KLD between two univariate Cauchy distributions was presented in [START_REF] Chyzak | A closed-form formula for the Kullback-Leibler divergence between Cauchy distributions[END_REF][START_REF] Nielsen | On f-divergences between Cauchy distributions[END_REF]. To date, the KLD of MCDs has no known explicit form, and it is in practice either estimated using expensive Monte Carlo stochastic integration or approximated. Monte Carlo sampling can efficiently estimate the KLD provided that a large number of independent and identically distributed samples is provided. Nevertheless, Monte Carlo integration is a too slow process to be useful in many applications. The main contribution of this paper is to derive a closed-form expression for the KLD between two central MCDs in a general case to benchmark future approaches while avoiding approximation using expensive Monte Carlo (MC) estimation techniques. The paper is organized as follows. Section 2 introduces the MCD and the KLD. Section 3 gives some definitions and propositions related to a multiple power series used to compute the closed-form expression of the KLD between two central MCDs. In Sections 4 and 5, expressions of some expectations related to the KLD are developed by exploiting the propositions presented in the previous section. Section 6 demonstrates some final results on the KLD computed for the central MCD. Section 7 presents some particular results such as the KLD for the univariate and the bivariate Cauchy distribution. Section 8 presents the implementation procedure of the KLD and a comparison with Monte Carlo sampling method. A summary and some conclusions are provided in the final section.

Multivariate Cauchy Distribution and Kullback-Leibler Divergence

Let X be a random vector of R p which follows the MCD, characterized by the following probability density function (pdf) given as follows [START_REF] Kotz | Multivariate T-Distributions and Their Applications[END_REF]

f X (x|µ, Σ, p) = Γ( 1+p 2 ) π p 2 Γ( 1 2 ) 1 |Σ| 1 2 1 [1 + (x -µ) T Σ -1 (x -µ)] 1+p 2 . ( 1
)
This is for any x ∈ R p , where p is the dimensionality of the sample space, µ is the location vector, Σ is a symmetric, positive definite (p × p) scale matrix and Γ(.) is the Gamma function. Let X 1 and X 2 be two random vectors that follow central MCDs with pdfs

f X 1 (x|Σ 1 , p) = f X 1 (x|0, Σ 1 , p) and f X 2 (x|Σ 2 , p) = f X 2 (x|0, Σ 2 , p)
given by [START_REF] Ollila | Complex Elliptically Symmetric Distributions: Survey, New Results and Applications[END_REF]. KLD provides an asymmetric measure of the similarity of the two pdfs. Indeed, the KLD between the two central MCDs is given by

KL(X 1 ||X 2 ) = R p ln f X 1 (x|Σ 1 , p) f X 2 (x|Σ 2 , p) f X 1 (x|Σ 1 , p)dx (2) = E X 1 {ln f X 1 (X)} -E X 1 {ln f X 2 (X)}. (3)
Since the KLD is the relative entropy defined as the difference between the cross-entropy and the entropy, we have the following relation:

KL(X 1 ||X 2 ) = H( f X 1 , f X 2 ) -H( f X 1 ) (4)
where

H( f X 1 , f X 2 ) = -E X 1 {ln f X 2 (X)} denotes the cross-entropy and H( f X 1 ) = -E X 1 {ln f X 1 (X)}
the entropy. Therefore, the determination of KLD requires the expression of the entropy and the cross-entropy. It should be noted that the smaller KL(X 1 ||X 2 ), the more similar are f X 1 (x|Σ 1 , p) and f X 2 (x|Σ 2 , p). The symmetric KL similarity measure between X 1 and X 2 is

d KL (X 1 , X 2 ) = KL(X 1 ||X 2 ) + KL(X 2 ||X 1
). In order to compute the KLD, we have to derive the analytical expressions of E

X 1 {ln f X 1 (X)} and E X 1 {ln f X 2 (X)} which depend, respec- tively, on E X 1 {ln[1 + X T Σ -1 1 X]} and E X 1 {ln[1 + X T Σ -1 2 X]}.
Consequently, the closed-form expression of the KLD between two zero-mean MCDs is given by

KL(X 1 ||X 2 ) = 1 2 log |Σ 2 | |Σ 1 | - 1 + p 2 E X 1 {ln[1 + X T Σ -1 1 X]} -E X 1 {ln[1 + X T Σ -1 2 X]} . ( 5
)
To provide the expression of these two expectations, some tools based on the multiple power series are required. The next section presents some definitions and propositions used for this goal.

Definitions and Propositions

This section presents some definitions and exposes some propositions related to the multiple power series used to derive the closed-form expression of the expectation

E X 1 {ln[1 + X T Σ -1 1 X]} and E X 1 {ln[1 + X T Σ -1
2 X]}, and as a consequence the KLD between two central MCDs.

Definition 1. The Humbert series of n variables, denoted

Φ (n)
2 , is defined for all x i ∈ C, i = 1, . . . , n, by the following multiple power series (Section 1.4 in [START_REF] Srivastava | Multiple Gaussian Hypergeometric Series; Ellis Horwood Series in Mathematics and Its Applications Statistics and Operational Research[END_REF])

Φ (n) 2 (b 1 , . . . , b n ; c; x 1 , . . . , x n ) = ∞ ∑ m 1 =0 .. ∞ ∑ m n =0 (b 1 ) m 1 . . . (b n ) m n (c) ∑ n i=1 m i n ∏ i=1 x m i i m i ! . ( 6
)
The Pochhammer symbol (q) i indicates the i-th rising factorial of q, i.e., for an integer i > 0 (q) i = q(q + 1) . . .

(q + i -1) = i-1 ∏ k=0 (q + k) = Γ(q + i) Γ(q) (7) 3.1. Integral Representation for Φ (n) 2
Proposition 1. The following integral representation is true for Real{c} > Real{∑ n i=1 b i } > 0 and Real{b i } > 0 where Real{.} denotes the real part of the complex coefficients

. . . ∆ 1 - n ∑ i=1 u i c- n ∑ i=1 b i -1 n ∏ i=1 u b i -1 i e x i u i du i = B b 1 , . . . , b n , c - n ∑ i=1 b i Φ (n) 2 (b 1 , . . . , b n ; c; x 1 , . . . , x n ) ( 8
)
where ∆ = {(u 1 , . . . , u n )|0 ≤ u i ≤ 1, i = 1, . . . , n; 0 ≤ u 1 + . . . + u n ≤ 1} and the multivariate beta function B is the extension of beta function to more than two arguments (called also Dirichlet function) defined as (Section 1.6.1 in [START_REF] Mathai | Special Functions for Applied Scientists[END_REF])

B(b 1 , . . . , b n , b n+1 ) = ∏ n+1 i=1 Γ(b i ) Γ(∑ n+1 i=1 b i ) . ( 9
)
Proof. The power series of exponential function is given by

e x i u i = ∞ ∑ m i =0 x m i i m i ! u m i i . ( 10
)
By substituting the expression of the exponential into the multiple integrals we have ..

∆ 1 - n ∑ i=1 u i c-∑ n i=1 b i -1 n ∏ i=1 u b i -1 i e x i u i du i = .. ∆ 1 - n ∑ i=1 u i c-∑ n i=1 b i -1 n ∏ i=1 ∞ ∑ m i =0 x m i i m i ! u m i +b i -1 i du i (11) = ∞ ∑ m 1 =0 .. ∞ ∑ m n =0 n ∏ i=1 x m i i m i ! × I D
where the multivariate integral I D , which is a generalization of a beta integral, is the type-1 Dirichlet integral (Section 1.6.1 in [START_REF] Mathai | Special Functions for Applied Scientists[END_REF]) given by

I D = . . . ∆ 1 - n ∑ i=1 u i c-∑ n i=1 b i -1 n ∏ i=1 u m i +b i -1 i du i = ∏ n i=1 Γ(b i + m i )Γ(c -∑ n i=1 b i ) Γ(c + ∑ n i=1 m i ) . ( 12
) Knowing that Γ(b i + m i ) = Γ(b i )(b i ) m i
, the expression of I D can be written otherwise

I D = ∏ n i=1 Γ(b i )Γ(c -∑ n i=1 b i ) Γ(c) ∏ n i=1 (b i ) m i (c) ∑ n i=1 m i . ( 13
)
Finally, plugging (13) back into [START_REF] Kullback | Information Theory and Statistics[END_REF] leads to the final result

Γ(c - n ∑ i=1 b i ) n ∏ i=1 Γ(b i ) Γ(c) +∞ ∑ m 1 ,..., m n =0 n ∏ i=1 (b i ) m i (c) ∑ n i=1 m i n ∏ i=1 x m i i m i ! = B b 1 , . . . , b n , c - n ∑ i=1 b i Φ (n) 2 (b 1 , . . . , b n ; c; x 1 , . . . , x n ) (14)
Given Proposition 1, we consider the particular cases n = {1, 2} one by one as follows:

Case n = 1 1 B(b 1 , c -b 1 ) 1 0 u b 1 -1 1 e x 1 u 1 (1 -u 1 ) c-b 1 -1 du 1 = ∞ ∑ m 1 =0 (b 1 ) m 1 (c) m 1 x m 1 1 m 1 ! = Φ (1) 2 (b 1 ; c; x 1 ) = 1 F 1 (b 1 , c; x 1 ) (15)
where 1 F 1 (.) is the confluent hypergeometric function of the first kind (Section 9.21 in [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]).

Case

n = 2 1 B(b 1 , b 2 , c -b 1 -b 2 ) u 1 ≥0,u 2 ≥0, u 1 +u 2 ≤1 u b 1 -1 1 u b 2 -1 2 e x 1 u 1 +x 2 u 2 (1 -u 1 -u 2 ) c-b 1 -b 2 -1 du 1 du 2 = ∞ ∑ m 1 =0 ∞ ∑ m 2 =0 (b 1 ) m 1 (b 2 ) m 2 (c) m 1 +m 2 x m 1 1 m 1 ! x m 2 2 m 2 ! = Φ (2) 2 (b 1 , b 2 ; c; x 1 , x 2 ) = Φ 2 (b 1 , b 2 , c; x 1 , x 2 ) ( 16
)
where the double series Φ 2 is one of the components of the Humbert series of two variables [START_REF] Humbert | The Confluent Hypergeometric Functions of Two Variables[END_REF] that generalize Kummer's confluent hypergeometric series 1 F 1 of one variable. The double series Φ 2 converges absolutely at any x 1 , x 2 ∈ C. N and given by

F (n) N (a; b 1 , . . . , b n ; c, c n ; x 1 , . . . , x n ) = x -a n +∞ ∑ m 1 ,..., m n =0 (a) ∑ n i=1 m i (a -c n + 1) ∑ n i=1 m i (a + b n -c n + 1) ∑ n i=1 m i ∏ n-1 i=1 (b i ) m i (c) ∑ n-1 i=1 m i n-1 ∏ i=1 x i x n m i 1 m i ! (1 -x -1 n ) m n m n ! . ( 17
)
The multiple power series [START_REF] Giet | A minimum Hellinger distance estimator for stochastic differential equations: An application to statistical inference for continuous time interest rate models[END_REF] is absolutely convergent on the region

|x i x -1 n | + |1 -x -1 n | < 1 in C n , ∀ i ∈ {1, . . . , n -1}.
The multiple power series F (n) N (.) can also be transformed into two other expressions as follows

F (n) N (a; b 1 , . . . , b n ; c, c n ; x 1 , . . . , x n ) = +∞ ∑ m 1 ,..., m n =0 (a -c n + 1) ∑ n-1 i=1 m i (b n ) m n (a) ∑ n i=1 m i (a + b n -c n + 1) ∑ n i=1 m i ∏ n-1 i=1 (b i ) m i (c) ∑ n-1 i=1 m i n-1 ∏ i=1 x m i i m i ! (1 -x n ) m n m n ! , (18)
= x 1-c n n +∞ ∑ m 1 ,..., m n =0 (a -c n + 1) ∑ n i=1 m i (b n -c n + 1) m n (a) ∑ n-1 i=1 m i (a + b n -c n + 1) ∑ n i=1 m i ∏ n-1 i=1 (b i ) m i (c) ∑ n-1 i=1 m i n-1 ∏ i=1 x m i i m i ! (1 -x n ) m n m n ! . ( 19
)
By Horn's rule for the determination of the convergence region (see [START_REF] Erdélyi | Higher Transcendental Functions[END_REF], Section 5.7.2), the multiple power series ( 18) and ( 19) are absolutely convergent on region

|x i | < 1, ∀ i ∈ {1, . . . , n -1}, |1 -x n | < 1 in C n .
Equation [START_REF] Csiszár | Eine informationstheoretische Ungleichung und ihre Anwendung auf den Beweis der Ergodizität von Markoffschen Ketten[END_REF] can then be deduced from ( 17) by using the following development where the

F (p)
N function can be written as

F (n) N (a; b 1 , . . . , b n ; c, c n ; x 1 , . . . , x n ) = x -a n +∞ ∑ m 1 ,..., m n-1 =0 (a) ∑ n-1 i=1 m i (a -c n + 1) ∑ n-1 i=1 m i (a + b n -c n + 1) ∑ n-1 i=1 m i ∏ n-1 i=1 (b i ) m i (c) ∑ n-1 i=1 m i × n-1 ∏ i=1 x i x n m i 1 m i ! ∞ ∑ m n =0 (α) m n (α -c n + 1) m n (α + b n -c n + 1) m n (1 -x -1 n ) m n m n ! (20)
and α = a + ∑ n-1 i=1 m i is used here to alleviate writing equations. Using the definition of Gauss' hypergeometric series 2 F 1 (.) [START_REF] Mathai | Special Functions for Applied Scientists[END_REF] and the Pfaff transformation [START_REF] Koepf | Hypergeometric Summation an Algorithmic Approach to Summation and Special Function Identities[END_REF], we can write

∞ ∑ m n =0 (α) m n (α -c n + 1) m n (α + b n -c n + 1) m n (1 -x -1 n ) m n m n ! = 2 F 1 α, α -c n + 1; α + b n -c n + 1; 1 -x -1 n (21) = x α n 2 F 1 α, b n ; α + b n -c n + 1; 1 -x n (22) = x α n ∞ ∑ m p =0 (α) m n (b n ) m n (α + b n -c n + 1) m n (1 -x n ) m n m n ! . ( 23
)
By substituting ( 23) into [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF], and using the following two relations:

(a) ∑ n-1 i=1 m i (α) m n = (a) ∑ n i=1 m i , (24)
(a + b n -c n + 1) ∑ n-1 i=1 m i (α + b n -c n + 1) m n = (a + b n -c n + 1) ∑ n i=1 m i (25)
we can get [START_REF] Csiszár | Eine informationstheoretische Ungleichung und ihre Anwendung auf den Beweis der Ergodizität von Markoffschen Ketten[END_REF].

The second transformation is given as follows

2 F 1 α, α -c n + 1; b n -c n + α + 1; 1 -x -1 n = x α-c n +1 n 2 F 1 b n -c n + 1, α -c n + 1; α + b n -c n + 1; 1 -x n (26) = x α-c n +1 n ∞ ∑ m n =0 (α -c n + 1) m n (b n -c n + 1) m n (α + b n -c n + 1) m n (1 -x n ) m n m n ! . ( 27
)
By substituting ( 27) into (20), we get [START_REF] Ali | A General Class of Coefficients of Divergence of One Distribution from Another[END_REF].

Lemma 1. The multiple power series F [START_REF] Lauricella | Sulle funzioni ipergeometriche a piu variabili[END_REF] when ac n + 1 = c and is given as follows

(n) N is equal to the Lauricella D-hypergeometric function F (n) D (see Appendix A)
F (n) N (a; b 1 , . . . , b n ; c, c n ; x 1 , . . . , x n ) = +∞ ∑ m 1 ,..., m n =0 (a) ∑ n i=1 m i ∏ n i=1 (b i ) m i (a + b n -c n + 1) ∑ n i=1 m i n-1 ∏ i=1 x m i i m i ! (1 -x n ) m n m n ! (28) = F (n) D (a, b 1 , . . . , b n ; a + b n -c n + 1; x 1 , . . . , x n-1 , 1 -x n ) (29)
Proof. By using Equation ( 18) of the multiple power series

F (n)
N and after having simplified (a

-c n + 1) ∑ n-1 i=1 m i to the numerator and (c) ∑ n-1 i=1 m i
to the denominator, we can get the result.

Integral Representation for F

(n+1) N Proposition 2. The following integral representation is true for Real{a} > 0, Real{a -c n+1 + 1} > 0, and Real{a -c n+1 + b n+1 + 1} > 0 Γ(a)Γ(a -c n+1 + 1) Γ(a -c n+1 + b n+1 + 1) F (n+1) N (a; b 1 , . . . , b n+1 ; c, c n+1 ; x 1 , . . . , x n+1 ) = ∞ 0 e -r r a-1 Φ (n) 2 (b 1 , . . . , b n ; c; rx 1 , . . . , rx n )U(b n+1 , c n+1 ; rx n+1 )dr ( 30
)
where U(•) is the confluent hypergeometric function of the second kind (Section 9.21 in [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]) defined for Real{b} > 0, Real{z} > 0 by the following integral representation

U(b, c; z) = 1 Γ(b) ∞ 0 e -zt t b-1 (1 + t) c-b-1 dt (31)
and

Φ (n) 2 (•) is defined by Equation (6).
Proof. The multiple power series

Φ (n)
2 and the confluent hypergeometric function U(•) are absolutely convergent on [0, +∞]. Using these functions in the above integral and changing the order of integration and summation, which is easily justified by absolute convergence, we get

∞ 0 e -r r a-1 Φ (n) 2 (b 1 , . . . , b n ; c; rx 1 , . . . , rx n )U(b n+1 , c n+1 ; rx n+1 )dr = ∞ ∑ m 1 =0 .. ∞ ∑ m n =0 (b 1 ) m 1 . . . (b n ) m n (c) ∑ n i=1 m i n ∏ i=1 x m i i m i ! I ( 32
)
where integral I is defined as follows

I = ∞ 0 e -r r a-1+ n ∑ i=1 m i U(b n+1 , c n+1 ; rx n+1 )dr. ( 33
)
Substituting the integral expression of U(•) in the previous equation and replacing α = a + ∑ n i=1 m i to alleviate writing equations, we have

I = 1 Γ(b n+1 ) ∞ 0 ∞ 0 e -(1+x n+1 t)r r α-1 t b n+1 -1 (1 + t) -(c n+1 -b n+1 -1) drdt. ( 34
)
Knowing that [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF] ∞

0 e -(1+x n+1 t)r r α-1 dr = Γ(α) (1 + x n+1 t) α ( 35
)
and

∞ 0 t b n+1 -1 (1 + t) c n+1 -b n+1 -1 (1 + x n+1 t) α dt = Γ(b n+1 )Γ(α -c n+1 + 1) Γ(α + b n+1 -c n+1 + 1) 2 F 1 (α, b n+1 ; α + b n+1 -c n+1 + 1; 1 -x n+1 ) (36)
the new expression of I is then given by

I = Γ(α)Γ(α -c n+1 + 1) Γ(α + b n+1 -c n+1 + 1) +∞ ∑ m n+1 =0 (α) m n+1 (b n+1 ) m n+1 (α + b n+1 -c n+1 + 1) m n+1 (1 -x n+1 ) m n+1 m n+1 ! . ( 37
)
Using the fact that

Γ(α) = Γ(a)(a) ∑ n i=1 m i and (a) ∑ n i=1 m i (α) m n+1 = (a) ∑ n+1 i=1 m i
, and developing the same method to Γ(α + b n+1c n+1 + 1), the final complete expression of the integral is then given by

Γ(a)Γ(a -c n+1 + 1) Γ(a + b n+1 -c n+1 + 1) ∞ ∑ m 1 =0 .. ∞ ∑ m n+1 =0 (b 1 ) m 1 . . . (b n ) m n (c) ∑ n i=1 m i (a -c n+1 + 1) ∑ n i=1 m i (b n+1 ) m n+1 (a) ∑ n+1 i=1 m i (a + b n+1 -c n+1 + 1) ∑ n+1 i=1 m i n ∏ i=1 x m i i m i ! × (1 -x n+1 ) m n+1 m n+1 ! = Γ(a)Γ(a -c n+1 + 1) Γ(a -c n+1 + b n+1 + 1) F (n+1) N (a; b 1 , . . . , b n+1 ; c, c n+1 ; x 1 , . . . , x n+1 ). ( 38
)

Expression of

E X 1 {ln[1 + X T Σ -1
1 X]} Proposition 3. Let X 1 be a random vector that follows a central MCD with pdf given by f

X 1 (x|Σ 1 , p). Expectation E X 1 {ln[1 + X T Σ -1 1 X]} is given as follows E X 1 ln[1 + X T Σ -1 1 X] = ψ 1 + p 2 -ψ 1 2 (39)
where ψ(.) is the digamma function defined as the logarithmic derivative of the Gamma function (Section 8.36 in [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]).

Proof. Expectation E X 1 {ln[1 + X T Σ -1 1 X]} is developed as follows E X 1 {ln[1 + X T Σ -1 1 X]} = A |Σ 1 | 1 2 R p ln[1 + x T Σ -1 1 x] [1 + x T Σ -1 1 x] 1+p 2 dx ( 40
)
where

A = Γ( 1+p 2 )π -1+p 2 .
Utilizing the following property log(x) f (x)dx = ∂ ∂a x a f (x) dx a=0 , as a consequence the expectation is given as follows

E X 1 {ln[1 + X T Σ -1 1 X]} = A |Σ 1 | 1 2 ∂ ∂a R p [1 + x T Σ -1 1 x] a-1+p 2 dx a=0 (41)
Consider the transformation y = Σ -1/2

1
x where y = [y 1 , y 2 , . . . , y p ] T . The Jacobian determinant is given by dy = |Σ 1 | -1/2 dx (Theorem 1.12 in [START_REF] Mathai | Jacobians of Matrix Transformations and Functions of Matrix Argument[END_REF]). The new expression of the expectation is given by

E X 1 {ln[1 + X T Σ -1 1 X]} = A ∂ ∂a R p [1 + y T y] a-1+p 2 dy a=0 . ( 42
)
Let u = y T y be a transformation where the Jacobian determinant is given by (Lemma 13.3.1 in [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF])

dy = π p 2 Γ( p 2 ) u p 2 -1 du. ( 43
)
The new expectation is as follows

E X 1 {ln[1 + X T Σ -1 1 X]} = Γ( 1+p 2 ) π 1/2 Γ( p 2 ) ∂ ∂a +∞ 0 u p 2 -1 (1 + u) a-1+p 2 du a=0 (44)
Using the definition of beta function, we can write that

+∞ 0 u p 2 -1 (1 + u) a-1+p 2 du = Γ( p 2 )Γ( 1 2 -a) Γ( 1+p 2 -a) . ( 45
)
The derivative of the last integral w.r.t a is given by

∂ ∂a +∞ 0 u p 2 -1 (1 + u) a-1+p 2 du a=0 = Γ( p 2 )Γ( 1 2 ) Γ( 1+p 2 ) ψ( 1 + p 2 ) -ψ( 1 2 ) (46)
Finally, the expression of

E X 1 {ln[1 + X T Σ -1 1 X]} is given by E X 1 {ln[1 + X T Σ -1 1 X]} = ψ 1 + p 2 -ψ 1 2 . ( 47
)

Expression of E

X 1 {ln[1 + X T Σ -1 2
X]} Proposition 4. Let X 1 and X 2 be two random vectors that follow central MCDs with pdfs given, respectively, by f

X 1 (x|Σ 1 , p) and f X 2 (x|Σ 2 , p). Expectation E X 1 {ln[1 + X T Σ -1 2 X]} is given as follows E X 1 {ln[1 + X T Σ -1 2 X]} = ψ 1 + p 2 -ψ 1 2 + ln λ p - ∂ ∂a F (p) D a, 1 2 , 1 2 , . . . , 1 2 , a + 1 2 p ; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p a=0 . ( 48
)
where λ 1 ,. . . , λ p are the eigenvalues of the real matrix Σ 1 Σ -1 2 , and F Proof. To prove Proposition 4, different steps are necessary. They are described in the following:

5.1. First Step: Eigenvalue Expression Expectation E X 1 {ln[1 + X T Σ -1
2 X]} is computed as follows

E X 1 {ln[1 + X T Σ -1 2 X]} = A |Σ 1 | 1 2 R p ln[1 + x T Σ -1 2 x] [1 + x T Σ -1 1 x] 1+p 2 dx ( 49
)
where

A = Γ( 1+p 2 )π -1+p 2 . Consider transformation y = Σ -1/2 1
x where y = [y 1 , y 2 , . . . , y p ] T . The Jacobian determinant is given by dy = |Σ 1 | -1/2 dx (Theorem 1.12 in [START_REF] Mathai | Jacobians of Matrix Transformations and Functions of Matrix Argument[END_REF]) and matrix

Σ = Σ 1 2 1 Σ -1 2 Σ 1 2
1 is a real symmetric matrix since Σ 1 and Σ 2 are real symmetric matrixes. Then, the expectation is evaluated as follows

E X 1 {ln[1 + X T Σ -1 2 X]} = A R p ln[1 + y T Σy] [1 + y T y] 1+p 2 dy. ( 50
)
Matrix Σ can be diagonalized by an orthogonal matrix P with P -1 = P T and Σ = PDP -1 where D is a diagonal matrix composed of the eigenvalues of Σ. Considering that y T Σy = tr(Σyy T ) = tr(PDP T yy T ) = tr(DP T yy T P), the expectation can be written as

E X 1 {ln[1 + X T Σ -1 2 X]} = A R p ln[1 + tr(DP T yy T P)] [1 + y T y] 1+p 2 dy. ( 51
)
Let z = P T y with z = [z 1 , z 2 , . . . , z p ] T be a transformation where the Jacobian determinant is given by dz = |P T |dy = dy. Using the fact that tr(DP T yy T P) = tr(Dzz T ) = z T Dz and y T y = z T P T Pz = z T z, then the previous expectation ( 51) is given as follows

E X 1 {ln[1 + X T Σ -1 2 X]} = A R p ln[1 + z T Dz] [1 + z T z] 1+p 2 dz (52) = A R .. R ln[1 + ∑ p i=1 λ i z 2 i ] [1 + ∑ p i=1 z 2 i ] 1+p 2 dz 1 . . . dz p (53)
where λ 1 ,. . . , λ p are the eigenvalues of Σ 1 Σ -1 2 .

Second Step: Polar Decomposition

Let the independent real variables z 1 , . . . , z p be transformed to the general polar coordinates r, θ 1 , . . . , θ p-1 as follows, where r > 0, -π/2 < θ j ≤ π/2, j = 1, . . . , p -2, -π < θ p-1 ≤ π [START_REF] Mathai | Jacobians of Matrix Transformations and Functions of Matrix Argument[END_REF],

z 1 = r sin θ 1 ( 54
)
z 2 = r cos θ 1 sin θ 2 (55)
z j = r cos θ 1 cos θ 2 . . . cos θ j-1 sin θ j , j = 2, 3, . . . , p -1 (56)

z p = r cos θ 1 cos θ 2 . . . cos θ p-1 . (57)
The Jacobian determinant according to theorem (1.24) in [START_REF] Mathai | Jacobians of Matrix Transformations and Functions of Matrix Argument[END_REF] is

dz 1 . . . dz p = r p-1 p-1 ∏ j=1 | cos θ j | p-j-1 drdθ j . ( 58
)
It is clear that with the last transformations, we get ∑ p i=1 z 2 i = r 2 and the multiple integral in (53) is then given as follows

E X 1 {ln[1 + X T Σ -1 2 X]} = A +∞ 0 r p-1 [1 + r 2 ] 1+p 2 π/2 -π/2 .. π -π p-1 ∏ j=1 | cos θ j | p-j-1 × ln 1 + r 2 (λ 1 sin 2 θ 1 + . . . + λ p cos 2 θ 1 . . . cos 2 θ p-1 ) dr p-1 ∏ j=1 dθ j . ( 59
)
By replacing the expression of sin 2 θ j by 1cos 2 θ j , for j = 1, . . . , p -1, we have the following expression

λ 1 sin 2 θ 1 + . . . + λ p cos 2 θ 1 . . . cos 2 θ p-1 = λ 1 + (λ 2 -λ 1 ) cos 2 θ 1 + . . . + (λ p -λ p-1 ) cos 2 θ 1 cos 2 θ 2 . . . cos 2 θ p-1 . (60)
Let x i = cos 2 θ i be a transformation to use where dx i = 2x 1/2 i (1x i ) 1/2 dθ i . Then the expectation given by the multiple integral over all θ j , j = 1, . . . , p -1 is as follows

2A +∞ 0 r p-1 [1 + r 2 ] 1+p 2 1 0 . . . 1 0 p-1 ∏ j=1 x p-j 2 -1 j (1 -x j ) -1 2 ln[1 + r 2 B p (x 1 , . . . , x p-1 )]drdx 1 . . . dx p-1 (61)
where

B p (x 1 , . . . , x p-1 ) = λ 1 + (λ 2 -λ 1 )x 1 + . . . + (λ p -λ p-1 )x 1 x 2 . . . x p-1 , p ≥ 1 and B 1 = λ 1 .
In the following, we use the notation B p instead of B p (x 1 , . . . , x p-1 ) to alleviate writing equations. Let t = r 2 be transformation to use. Then, one can write

= A +∞ 0 t p 2 -1 [1 + t] 1+p 2 1 0 . . . 1 0 p-1 ∏ j=1 x p-j 2 -1 j (1 -x j ) -1 2 ln[1 + tB p ]dtdx 1 . . . dx p-1 . (62)
In order to solve the integral in (62), we consider the following property given by log(x) f (x) dx = -∂ ∂a x -a f (x)dx a=0 and the following equation given as follows

1 + B p t -a = 1 Γ(a) +∞ 0 y a-1 e -(1+B p t)y dy. ( 63
)
Making use of the above equation, we obtain a new expression of (62) given as follows

E X 1 {ln[1 + X T Σ -1 2 X]} = - ∂ ∂a A Γ(a) +∞ 0 t p 2 -1 [1 + t]
∏ j=1 x p-j 2 -1 j (1 -x j ) -1 2 dx j dydt a=0 (64) = - ∂ ∂a A Γ(a) +∞ 0 t p 2 -1 [1 + t]
where H(t, y) is defined as

H(t, y) = 1 0 . . . 1 0 e -B p ty p-1 ∏ j=1 x p-j 2 -1 j (1 -x j ) -1 2 dx j . (66)

Third Step: Expression for H(t,y) by Humbert and Beta Functions

Let

x i = 1 -x i , i = 1, . . . , p -1 be transformations to use. Then (λ 2 -λ 1 )x 1 = (λ 2 -λ 1 )(1 -x 1 ) (67) (λ 3 -λ 2 )x 1 x 2 = (λ 3 -λ 2 )(1 -x 1 )[1 -x 2 ] (68) (λ 4 -λ 3 )x 1 x 2 x 3 = (λ 4 -λ 3 )(1 -x 1 )(1 -x 2 )[1 -x 3 ] (69) . . . = . . . (λ p -λ p-1 ) p-1 ∏ i=1 x i = (λ p -λ p-1 ) p-1 ∏ i=1 (1 -x i ). (70)
Adding equations from (67) to (70), we can state that the new expression of the function B p becomes

B p = λ p -(λ p -λ 1 )x 1 -(λ p -λ 2 )(1 -x 1 )x 2 -(λ p -λ 3 )(1 -x 1 )(1 -x 2 )x 3 -. . . -(λ p -λ p-1 )(1 -x 1 ) . . . (1 -x p-2 )x p-1 . (71)
Then, the multiple integral H(t, y) given by (66) can be written otherwise

H(t, y) = 1 0 . . . 1 0 e -B p ty p-1 ∏ j=1 (1 -x j ) p-j 2 -1 x -1 2 j dx 1 . . . dx p-1 . ( 72
)
Let the real variables x 1 , x 2 , . . . , x p-1 be transformed to the real variables u 1 , u 2 , . . . , u p-1 as follows

u 1 = x 1 ( 73
)
u 2 = (1 -x 1 )x 2 = (1 -u 1 )x 2 ( 74
)
u 3 = (1 -x 1 )(1 -x 2 )x 3 = (1 -u 1 -u 2 )x 3 (75) . . . u p-1 = p-2 ∏ i=1 (1 -x i )x p-1 = (1 - p-2 ∑ i=1 u i )x p-1 . ( 76
)
The Jacobian determinant is given by du

1 . . . du p-1 = p-1 ∏ j=1 1 - j-1 ∑ i=1 u i dx 1 . . . dx p-1 . (77)
Accordingly, the new expression of B p becomes

B p = λ p - p-1 ∑ i=1 (λ p -λ i )u i . ( 78
)
As a consequence, the new domain of the multiple integral (72

) is ∆ = {(u 1 , u 2 , . . . , u p-1 ) ∈ R p-1 ; 0 ≤ u 1 ≤ 1, 0 ≤ u 2 ≤ 1 -u 1 , 0 ≤ u 3 ≤ 1 -u 1 -u 2 , . . . , and 0 ≤ u p-1 ≤ 1 -u 1 - u 2 .
. .u p-2 }, and the expression of H(t, y) is given as follows

H(t, y) = . . . ∆ e -B p ty p-1 ∏ j=1 1 - j-1 ∑ i=1 u i -1 u j 1 -∑ j-1 i=1 u i -1 2 p-1 ∏ j=1 1 - u j 1 -∑ j-1 i=1 u i p-j 2 -1 du j (79) = . . . ∆ e -B p ty p-1 ∏ j=1 u -1 2 j 1 - j ∑ i=1 u i p-j 2 -1 1 - j-1 ∑ i=1 u i 1 2 - p-j 2 du 1 . . . du p-1 (80) = . . . ∆ e -B p ty 1 - p-1 ∑ i=1 u i p 2 - p-1 2 -1 p-1 ∏ j=1 u -1 2 j du j (81)
= e -λ p ty . . .

∆ 1 - p-1 ∑ i=1 u i -1 2 p-1 ∏ i=1 u -1 2 i e (λ p -λ i )u i ty du i . ( 82
)
Using Proposition 1, we subsequently find that (.) is the Humbert series of p -1 variables and B( 1 2 , . . . , 1 2 ) is the multivariate beta function. Applying the following successive two transformations r = ty (dr = tdy) and u = 1/t (du = -u 2 dt), the new expression of the expectation given by ( 65) is written as follows

H(t, y) = e -λ p
E X 1 {ln[1 + X T Σ -1 2 X]} = - ∂ ∂a A Γ(a) B 1 2 , . . . , 1 2 p +∞ 0 r a-1 e -λ p r × Φ (p-1) 2 1 2 , . . . , 1 2 p-1 ; p 2 ; (λ p -λ 1 )r, . . . , (λ p -λ p-1 )r +∞ 0 u a-1 2 (1 + u) -1+p 2 e -ru du dr a=0 . ( 84
)

Final Step

The last integral is related to the confluent hypergeometric function of the second kind U(.) as follows

+∞ 0 u a-1 2 (1 + u) -1+p 2 e -ru du = Γ(a + 1 2 )U(a + 1 2 , a + 1 - p 2 , r). ( 85
)
As a consequence, the new expression is

E X 1 {ln[1 + X T Σ -1 2 X]} = - ∂ ∂a A Γ(a + 1 2 ) Γ(a) B 1 2 , . . . , 1 2 × +∞ 0 r a-1 e -λ p r Φ (p-1) 2 1 2 , . . . , 1 2 ; p 2 ; 1; (λ p -λ 1 )r, . . . , (λ p -λ p-1 )r U(a + 1 2 , a + 1 - p 2 , r)dr a=0 . ( 86
)
Using the transformation r = λ p r and the Proposition 2, and taking into account the expression of A, the new expression becomes

E X 1 {ln[1 + X T Σ -1 2 X]} = - ∂ ∂a B(a + 1 2 , p 2 ) B( 1 2 , p 2 ) λ -a p × F (p) N a; 1 2 , . . . , 1 2 , a + 1 2 p ; p 2 , a - p 2 + 1; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , λ -1 p a=0 (87)
Knowing that

∂ ∂a B( p 2 , a + 1 2 ) B( p 2 , 1 2 ) a=0 = ψ 1 2 -ψ 1 + p 2 , and (88)
F (p) N a; 1 2 , . . . , 1 2 , a + 1 2 ; p 2 , a - p 2 + 1; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , λ -1 p a=0 = 1, ( 89
) the new expression of E X 1 {ln[1 + X T Σ -1 2 X]} becomes E X 1 {ln[1 + X T Σ -1 2 X]} = ψ 1 + p 2 -ψ 1 2 - ∂ ∂a λ -a p F (p) N a; 1 2 , . . . , 1 2 , a + 1 2 p ; p 2 , a - p 2 + 1; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , λ -1 p a=0 . ( 90
)
Applying the expression given by ( 18) of Definition 2 and relying on Lemma 1, the final result corresponds to the D-hypergeometric function of Lauricella F (p) D (.) given by

E X 1 {ln[1 + X T Σ -1 2 X]} = ψ 1 + p 2 -ψ 1 2 - ∂ ∂a λ -a p +∞ ∑ m 1 ,..., m p =0 (a) ∑ p i=1 m i (a + 1 2 ) m p ∏ p-1 i=1 ( 1 2 ) m i (a + 1+p 2 ) ∑ p i=1 m i p-1 ∏ i=1 1 - λ i λ p m i 1 m i ! (1 -λ -1 p ) m p m p ! a=0 (91) = ψ 1 + p 2 -ψ 1 2 - ∂ ∂a λ -a p F (p) D a, 1 2 , . . . , 1 2 , a + 1 2 p ; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p a=0 . ( 92
)
The final development of the previous expression is as follows

E X 1 {ln[1 + X T Σ -1 2 X]} = ψ 1 + p 2 -ψ 1 2 + ln λ p - ∂ ∂a F (p) D a, 1 2 , 1 2 , . . . , 1 2 , a + 1 2 p ; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p a=0 . ( 93
)
In this section, we presented the exact expression of

E X 1 {ln[1 + X T Σ -1 2 X]}.
In addition, the multiple power series F (p) D which appears to be a special case of F (p) N provides many properties and numerous transformations (see Appendix A) that make easier the convergence of the multiple power series. In the next section, we establish the KLD closed-form expression based on the expression of the latter expectation.

KLD between Two Central MCDs

Plugging [START_REF] Lauricella | Sulle funzioni ipergeometriche a piu variabili[END_REF] and (93) into (5) yields the closed-form expression of the KLD between two central MCDs with pdfs f X 1 (x|Σ 1 , p) and f X 2 (x|Σ 2 , p). This result is presented in the following theorem.

Theorem 1. Let X 1 and X 2 be two random vectors that follow central MCDs with pdfs given, respectively, by f X 1 (x|Σ 1 , p) and f X 2 (x|Σ 2 , p). The Kullback-Leibler divergence between central MCDs is

KL(X 1 ||X 2 ) = - 1 2 log p ∏ i=1 λ i + 1 + p 2 log λ p - ∂ ∂a F (p) D a, 1 2 , . . . , 1 2 , a + 1 2 p ; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p a=0 ( 94
)
where λ 1 ,. . . , λ p are the eigenvalues of the real matrix Σ 1 Σ -1 2 , and F (p) D (.) represents the Lauricella D-hypergeometric function defined for p variables.

Lauricella [START_REF] Lauricella | Sulle funzioni ipergeometriche a piu variabili[END_REF] gave several transformation formulas (see Appendix A), whose relations (A5)-(A7), and (A9) are applied to F (p) D (.) in (94). The results of transformation are as follows

F (p) D a, 1 2 , . . . , 1 2 , a + 1 2 ; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p = λ a+ p 2 p p-1 ∏ i=1 λ -1 2 i F (p) D 1 + p 2 , 1 2 , . . . , 1 2 , a + 1 2 ; a + 1 + p 2 ; 1 - λ p λ 1 , . . . , 1 - λ p λ p-1 , 1 -λ p (95) = λ 1 λ p -a F (p) D a, 1 2 , . . . , 1 2 , a + 1 2
; a + 1 + p 2 ; 1 - λ p λ 1 , . . . , 1 - λ 2 λ 1 , 1 - 1 λ 1 (96) = λ a p F (p) D a, 1 2 , . . . , 1 2
; a + 1 + p 2 ; 1 -λ 1 , 1 -λ 2 , . . . , 1 -λ p (97) = λ a p p ∏ i=1 λ -1 2 i F (p) D 1 + p 2 , 1 2 , . . . , 1 2 ; a + 1 + p 2 ; 1 - 1 λ 1 , 1 - 1 λ 2 , . . . , 1 - 1 λ p . ( 98
)
Considering the above equations, it is easy to provide different expressions of KL(X 1 ||X 2 )

shown in Table 1. The derivative of the Lauricella D-hypergeometric series with respect to a goes through the derivation of the following expression

∂ ∂a F (p) D a, 1 2 , 1 2 , . . . , 1 2 , a + 1 2
; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p a=0 (99) = +∞ ∑ m 1 ,..., m p =0 ∂ ∂a (a) ∑ p i=1 m i (a + 1 2 ) m p (a + 1+p 2 ) ∑ p i=1 m i a=0 p-1 ∏ i=1 1 2 m i 1 - λ i λ p m i 1 m i ! (1 -λ -1 p ) m p m p ! (100)
The derivative with respect to a of the Lauricella D-hypergeometric series and its transformations goes through the following expressions (see Appendix B for demonstration)

∂ ∂a (a) ∑ p i=1 m i (a + 1 2 ) m p (a + 1+p 2 ) ∑ p i=1 m i a=0 = ( 1 2 ) m p (1) ∑ p i=1 m i ( 1+p 2 ) ∑ p i=1 m i (∑ p i=1 m i ) , ( 101
) ∂ ∂a (a) ∑ p i=1 m i (a + 1+p 2 ) ∑ p i=1 m i a=0 = (1) ∑ p i=1 m i ( 1+p 2 ) ∑ p i=1 m i (∑ p i=1 m i ) , (102)
∂ ∂a

(a + 1 2 ) m p (a + 1+p 2 ) ∑ p i=1 m i a=0 = ( 1 2 ) m p ( 1+p 2 ) ∑ p i=1 m i m p -1 ∑ k=0 1 k + 1 2 - ∑ p i=1 m i -1 ∑ k=0 1 k + 1+p 2 , (103)
∂ ∂a

1 (a + 1+p 2 ) ∑ p i=1 m i a=0 = - 1
( 1+p 2 ) ∑ p i=1 m i ∑ p i=1 m i -1 ∑ k=0 1 k + 1+p 2 . ( 104
)
To derive the closed-form expression of d KL (X 1 , X 2 ) we have to evaluate the expression of KL(X 2 ||X 1 ). The latter can be easily deduced from KL(X 1 ||X 2 ) as follows

KL(X 2 ||X 1 ) = 1 2 log p ∏ i=1 λ i - 1 + p 2 log λ p + ∂ ∂a F (p) D a, 1 2 , . . . , 1 2 , a + 1 2 ; a + 1 + p 2 ; 1 - λ p λ 1 , . . . , 1 - λ p λ p-1 , 1 -λ p a=0 . (105)
Proceeding in the same way by using Lauricella transformations, different expressions of KL(X 2 ||X 1 ) are provided in Table 1. Finally, given the above results, it is straightforward to compute the symmetric KL similarity measure d KL (X 1 , X 2 ) between X 1 and X 2 . Technically, any combination of the KL(X 1 ||X 2 ) and KL(X 2 ||X 1 ) expressions is possible to compute d KL (X 1 , X 2 ). However, we choose the same convergence region for the two divergences for the calculation of the distance. Some expressions of d KL (X 1 , X 2 ) are given in Table 1.

Table 1. KLD and KL distance computed when X 1 and X 2 are two random vectors following central MCDs with pdfs f X 1 (x|Σ 1 , p) and f X 2 (x|Σ 2 , p). (115)

KL(X 1 ||X 2 ) = - 1 2 log p ∏ i=1 λ i + 1 + p 2 log λ p - ∂ ∂a F (p) D a, 1 2 , . . . , 1 2 , a + 1 2 p
; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p a=0 (106) = - 1 2 log p ∏ i=1 λ i - 1 + p 2 λ p 2 p p-1 ∏ i=1 λ -1 2 i ∂ ∂a F (p) D 1 + p 2 , 1 2 , . . . , 1 2 , a + 1 2 p
; a + 1 + p 2 ; 1 - λ p λ 1 , . . . , 1 - λ p λ p-1 , 1 -λ p a=0 (107) = - 1 2 log p ∏ i=1 λ i + 1 + p 2 log λ 1 - ∂ ∂a F (p) D a, 1 2 , . . . , 1 2 , a + 1 2 p
; a + 1 + p 2 ; 1 - λ p λ 1 , . . . , 1 - λ 2 λ 1 , 1 - 1 λ 1 a=0 (108) = - 1 2 log p ∏ i=1 λ i - 1 + p 2 ∂ ∂a F (p) D a, 1 2 , . . . , 1 2 p
; a + 1 + p 2 ; 1 -λ 1 , . . . , 1 -λ p a=0 (109) = - 1 2 log p ∏ i=1 λ i - 1 + p 2 p ∏ i=1 λ -1 2 i ∂ ∂a F (p) D 1 + p 2 , 1 2 , . . . , 1 2 p
; a + 1 + p 2 ; 1 - 1 λ 1 , . . . , 1 - 1 λ p a=0 (110) KL(X 2 ||X 1 ) = 1 2 p ∏ i=1 λ i - 1 + p 2 log λ p + ∂ ∂a F (p) D a, 1 2 , . . . , 1 2 , a + 1 2 p
; a + 1 + p 2 ; 1 - λ p λ 1 , . . . , 1 - λ p λ p-1 , 1 -λ p a=0 (111) = 1 2 log p ∏ i=1 λ i - 1 + p 2 λ - p 2 p p-1 ∏ i=1 λ 1 2 i ∂ ∂a F (p) D 1 + p 2 , 1 2 , . . . , 1 2 , a + 1 2 p
; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p a=0 (112) = 1 2 log p ∏ i=1 λ i - 1 + p 2 log λ 1 + ∂ ∂a F (p) D a, 1 2 ,
d KL (X 1 , X 2 ) = 1 + p 2 log λ p - ∂ ∂a F (p) D a, 1 2 , . . . , 1 2 , a + 1 2 p ; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p a=0 -λ - p 2 p p-1 ∏ i=1 λ 1 2 i × ∂ ∂a F (p) D 1 + p 2 , 1 2 , . . . , 1 2 , a + 1 2 p ; a + 1 + p 2 ; 1 - λ 1 λ p , . . . , 1 - λ p-1 λ p , 1 - 1 λ p a=0 (116) = - 1 + p 2 ∂ ∂a F (p) D a, 1 2 , . . . , 1 2 p
; a + 1 + p 2 ; 1 -λ 1 , . . . , 1 -λ p a=0 + p ∏ i=1 λ 1 2 i ∂ ∂a F (p) D 1 + p 2 , 1 2 , . . . , 1 2 p
; a + 1 + p 2 ; 1 -λ 1 , . . . , 1 -λ p a=0 (117) = - 1 + p 2 p ∏ i=1 λ -1 2 i ∂ ∂a F (p) D 1 + p 2 , 1 2
∂ ∂a 2 F 1 (a, 1
2 ; a + 1; 1 -λ) a=0 = ∞ ∑ n=1 1 2 n 1 n (1 -λ) n n! = -2 ln 1 + λ 1/2 2 . ( 120
)
Accordingly, the KLD is then expressed as

KL(X 1 ||X 2 ) = log (1 + λ 1 2 ) 2 4λ 1 2 (121) = log (1 + λ -1 2 ) 2 4λ -1 2 = KL(X 2 ||X 1 ). ( 122
)
We conclude that KLD between Cauchy densities is always symmetric. Interestingly, this is consistent with the result presented in [START_REF] Chyzak | A closed-form formula for the Kullback-Leibler divergence between Cauchy distributions[END_REF].

Case of p = 2

This case corresponds to the Bivariate Cauchy distribution. The KLD is then given by

KL(X 1 ||X 2 ) = - 1 2 log λ 1 λ 2 - 3 2 ∂ ∂a F 1 (a, 1 2 , 1 2
; a + 3

2 ; 1 -λ 1 , 1 -λ 2 ) a=0 (123)
where F 1 is the Appell's hypergeometric function (see Appendix A). The expression of the derivative of F 1 can be further developed

∂ ∂a F 1 (a, 1 2 , 1 2 ; a + 3 2 ; 1 -λ 1 , 1 -λ 2 ) a=0 = +∞ ∑ n,m=0 (
1) m+n ( 1 2 ) n ( 1 2 ) m ( 3 2 ) m+n 1 m + n (1 -λ 1 ) n n! (1 -λ 2 ) m m! . ( 124
)
In addition, when the eigenvalue λ i for i = 1, 2 takes some particular values, the expression of the KLD becomes very simple. In the following, we show some cases:

(λ 1 = 1, λ 2 = 1) or (λ 2 = 1, λ 1 = 1)
For this particular case, we have

∂ ∂a F 1 (a, 1 2 , 1 2 ; a + 3 2 ; 1 -λ i , 0) a=0 = ∂ ∂a 2 F 1 (a, 1
2 ; a + 3 2 ; 1 -λ i ) a=0 (125) = -ln λ i + 1 √ 1 -λ i ln 1 - √ 1 -λ i 1 + √ 1 -λ i + 2. ( 126
)
The demonstration of the derivation is shown in Appendix C.2. Then, KLD becomes equal to

KL(X 1 ||X 2 ) = ln λ i - 3 2 1 √ 1 -λ i ln 1 - √ 1 -λ i 1 + √ 1 -λ i -3. ( 127
) λ 1 = λ 2 = λ
For this particular case, we have

∂ ∂a F 1 (a, 1 2 , 1 2 ; a + 3 2 ; 1 -λ, 1 -λ) a=0 = ∂ ∂a 2 F 1 (a, 1; a + 3 2 ; 1 -λ) a=0 (128) = - 2 √ 1 -λ -1 ln( √ λ + √ λ -1) + 2. ( 129
)
For more details about the demonstration see Appendix C.3. The KLD becomes equal to

KL(X 1 ||X 2 ) = -ln λ + 3 √ 1 -λ -1 ln( √ λ + √ λ -1) -3. ( 130
)
It is easy to deduce that

KL(X 2 ||X 1 ) = ln λ + 3 √ 1 -λ ln( √ λ -1 + λ -1 -1) -3. ( 131
)
This result can be demonstrated using the same process as KL(X 1 ||X 2 ). It is worth to notice that KL(X 1 ||X 2 ) = KL(X 2 ||X 1 ) which leads us to conclude that the property of symmetry observed for the univariate case is no longer valid in the multivariate case. Nielsen et al. in [START_REF] Nielsen | On f-divergences between Cauchy distributions[END_REF] gave the same conclusion.

Implementation and Comparison with Monte Carlo Technique

In this section, we show how we practically compute the numerical values of the KLD, especially when we have several equivalent expressions which differ in the region of convergence. To reach this goal, the eigenvalues of Σ 1 Σ -1 2 are rearranged in a descending order λ p > λ p-1 > . . . > λ 1 > 0. This operation is justified by Equation ( 53) where it can be seen that the permutation of the eigenvalues does not affect the expectation result. Three cases can be identified from the expressions of KLD.

8.1. Case 1 > λ p > λ p-1 > . . . > λ 1 > 0
The expression of KL(X 1 ||X 2 ) is given by Equation (109) and KL(X 2 ||X 1 ) is given by (115).

Case λ

p > λ p-1 > . . . > λ 1 > 1 KL(X 1 ||X 2
) is given by the Equation (110) and KL(X 2 ||X 1 ) is given by (114).

Case λ

p > 1 and λ 1 < 1 This case guarantees that 0 ≤ 1 -λ j /λ p < 1, j = 1, . . . , p -1 and 0 ≤ 1 -1/λ p < 1.
The expression of the KL(X 1 ||X 2 ) is given by Equation (106) and KL(X 2 ||X 1 ) is given by (112) or (113). To perform an evaluation of the quality of the numerical approximation of the derivative of the Lauricella series, we consider a case where an exact and simple expression of ∂ ∂a {F (p) D (.)}| a=0 is possible. The following case where λ 1 = . . . = λ p = λ allows the Lauricella series to be equivalent to the Gauss hypergeometric function given as follows

F (p) D a, 1 2 , . . . , 1 2 p
; a + 1 + p 2 ; 1 -λ, . . . , 1 -λ = 2 F 1 a, p 2 ; a + 1 + p 2 ; 1 -λ . ( 132
)
This relation allows us to compare the computational accuracy of the approximation of the Lauricella series with respect to the Gauss function. In addition, to compute the numerical value the indices of the series will evolve from 0 to N instead of infinity. The latter is chosen to ensure a good approximation of the Lauricella series. Table 2 shows the computation of the derivative of F (p) D (.) and 2 F 1 (.), along with the absolute value of error | |, where p = 2, N = {20, 30, 40}. The exact expression of ∂ ∂a { 2 F 1 (.)}| a=0 when p = 2 is given by Equation (129). We can deduce the following conclusions. First, the error is reasonably low and decreases as the value of N increases. Second, the error increases for values of 1λ close to 1 as expected, which corresponds to the convergence region limit. In the following section, we compare the Monte Carlo sampling method to approximate the KLD value with the numerical value of the closed-form expression of the latter. The Monte Carlo method involves sampling a large number of samples and using them to calculate the sum rather than the integral. Here, for each sample size, the experiment is repeated 2000 times. The elements of Σ 1 and Σ 2 are given in Table 3. Figure 1 depicts the absolute value of bias, mean square error (MSE) and box plot of the difference between the symmetric KL approximated value and its theoretical one, given versus the sample sizes. As the sample size increases, the bias and the MSE decrease. Accordingly, the approximated value will be very close to the theoretical KLD when the number of samples is very large. The computation time of the proposed approximation and the classical Monte Carlo sampling method are recorded using Matlab on a 1.6 GHz processor with 16 GB of memory. For the proposed numerical approximation, the computation time is evaluated to 1.56 s with N = 20. The value of N can be increased to further improve the accuracy, but it will increase the computation time. For the Monte Carlo sampling method, the mean time values at sample sizes of {65,536; 131,072; 262,144} are {2.71; 5.46; 10.78} seconds, respectively. D . We have also proposed a simple scheme to compute easily the Lauricella series and to bypass the convergence constraints of this series. Codes and examples for numerical calculations are presented and explained in detail. Finally, a comparison is made to show how the Monte Carlo sampling method gives approximations close to the KLD theoretical value. As a final note, it is also possible to extend these results on the KLD to the case of the multivariate t-distribution since the MCD is a particular case of this multivariate distribution.

Lauricella has given several transformation formulas, from which we use the two following relationships. More details can be found in Exton's book [START_REF] Exton | Multiple Hypergeometric Functions and Applications[END_REF] on hypergeometric equations. x nx 1

x n -1 ,

x nx 2

x n -1 , . . . ,

x n -x n-1 x n -1 , x n x n -1 (A7) = (1 -x 1 ) c-a n ∏ i=1 (1 -x i ) -b i F (n) D c -a, c - n ∑ i=1 b i , b 2 , . . . , b n ; c; x 1 , x 2 -x 1 x 2 -1 , . . . , x n -x 1 x n -1 (A8) = (1 -x n ) c-a n ∏ i=1 (1 -x i ) -b i F (n) D c -a,

D

  (.) represents the Lauricella D-hypergeometric function defined for p variables.

  1 e -y H(t, y)dydt a=0

Figure 1 .

 1 Figure 1. Top row: Bias (left) and MSE (right) of the difference between the approximated and theoretical symmetric KL for MCD. Bottom row: Box plot of the error. The mean error is the bias. Outliers are larger than Q 3 + 1.5 × IQR or smaller than Q 1 -1.5 × IQR, where Q 1 , Q 3 , and IQR are the 25th, 75th percentiles, and the interquartile range, respectively.

FD( 1 -= ( 1 -x 1 )= ( 1 -

 1111 (a, b 1 , . . . , b n ; c; x 1 , . . . , x n ) = n ∏ i=1 x i ) -b i F (n) D ca, b 1 , . . . , b n ; x n ) -a F (n) D a, b 1 , . . . , b n-1 , c -

Table 2 .

 2 Computation of A = ∂ ∂a { 2 F 1 (.)}| a=0 and B = ∂ ∂a {F (.)}| a=0 when p = 2 and λ 1 = . . . = λ p = λ.

	(p)

D

  b 1 , . . . , b n-1 , c -

	n ∑ i=1	b i ; c;	x 1 -x n x 1 -1	, . . . ,	x n-1 -x n x n-1 -1	, x n .	(A9)
	Appendix B. Demonstration of Derivative

Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.

Data Availability Statement: Not applicable.

Acknowledgments: Authors gratefully acknowledge the PHENOTIC platform node of the french national infrastructure on plant phenotyping ANR PHENOME 11-INBS-0012. The authors would like also to thank the anonymous reviewers for their helpful comments valuable comments and suggestions.

To further encourage the dissemination of these results, we provide a code available as attached file to this paper. This is given in Matlab with a specific case of p = 3. This can be easily extended to any value of p, thanks to the general closed-form expression established in this paper. ---------------------------------J --------------------------------- -------------------------------------------H ------------------------------------------ -------------------------------------------J -------------------------------------------

Conflicts of Interest:

The authors declare no conflict of interest.

Appendix A. Lauricella Function

In 1893, G. Lauricella [START_REF] Lauricella | Sulle funzioni ipergeometriche a piu variabili[END_REF] investigated the properties of four series

of n variables. When n = 2, these functions coincide with Appell's F 2 , F 3 , F 4 , F 1 , respectively. When n = 1, they all coincide with Gauss' 2 F 1 . We present here only the Lauricella series F (n) D given as follows

where |x 1 |, . . . , |x | < 1. The Pochhammer symbol (q) i indicates the i-th rising factorial of q, i.e., (q) i = q(q + 1) . . .

) can be expressed in terms of multiple integrals as follows [START_REF] Hattori | On the Euler integral representations of hypergeometric functions in several variables[END_REF]

where

Lauricella's F D can be written as a one-dimensional Euler-type integral for any number n of variables. The integral form of F

(n) D (.) is given as follows when Real(a) > 0 and Real(ca) > 0

We use the following notation α = ∑ p i=1 m i to alleviate the writing of equations.

Using the fact that

we can state that

By developing the previous expression we can state that

As a consequence, ∂ ∂a Let f be a function of λ defined as

The multiplication of the derivative of f with respect to λ by (1λ) is given as follows

As a consequence,

where f is a function of λ i . The multiplication of the derivative of f with respect to λ i by (1

we can deduce an expression of

where f is a function of λ. The multiplication of the derivative of f with respect to λ by (1λ) is given as follows

Knowing that 2 F 1 1, 1;

we can state that

As a consequence,


Ayman Beghdadi
email: aymanaymar.beghdadi@univ-evry.fr

Malik Mallem
email: malik.mallem@univ-evry.fr

Lotfi Beji
email: lotfi.beji@univ-evry.fr

BENCHMARKING PERFORMANCE OF OBJECT DETECTION UNDER IMAGES DISTORTIONS OF AN UNCONTROLLED ENVIRONMENT

Keywords: Deep learning, Object detection, Distortion, Robustness, Benchmark

The robustness of object detection algorithms plays a prominent role in real-world applications, especially in the uncontrolled environments due to distortions during image acquisition. It has been proved that object detection methods suffer from in-capture distortions to perform a reliable detection. In this study, we present a performance evaluation framework of the state-of-the-art object detection methods on a dedicated dataset containing images with various distortions at different levels of severity. Furthermore, we propose an original strategy of image distortion generation applied to the MS-COCO dataset that combines some local and global distortions to reach a better realism. We have shown that training with this distorted dataset improves the robustness of models by 31.5%. Finally, we provided a custom dataset including the natural images distorted from MS-COCO to perform a more relevant evaluation of the robustness concerning distortions. The database and the generation source codes of the different distortions are publicly available 1,2 .

INTRODUCTION

Object detection and recognition is one of the most challenging and widely studied problems in computer vision and especially in uncontrolled environment. Deep Neural Networks (DNN) based approaches [START_REF] Zou | Object detection in 20 years: A survey[END_REF][START_REF] Zhao | Object detection with deep learning: A review[END_REF][START_REF] Liu | Deep learning for generic object detection: A survey[END_REF] have been proven very promising paradigms to solve such complex problems in various real-world applications. However, generally, the proposed DNN-based architectures do not consider the sensitivity of the learning models to common data distortions and corruption effects. Indeed, it has been shown by Szegedy et al. [START_REF] Szegedy | Intriguing properties of neural networks[END_REF] that a small perturbation invisible to the human visual system can have significant repercussions on the performance of object recognition models. It is therefore important to conduct ad hoc studies to understand the effects of such perturbations and to be able to generate them at different levels in order to develop more robust object detection networks. The limitations of object detection methods, in terms of robustness against some common signal degradation, data corruptions, or adversarial examples, have been highlighted in several works [START_REF] Michaelis | Benchmarking robustness in object detection: Autonomous driving when winter is coming[END_REF][START_REF] Hendrycks | Benchmarking neural network robustness to common corruptions and perturbations[END_REF][START_REF] Li | Universal adversarial perturbations against object detection[END_REF][START_REF] Benz | Revisiting batch normalization for improving corruption robustness[END_REF]. A solution would be to perform pre-processing before the high-level phases. But this implies being able to identify the distortion automatically in order to apply the most appropriate solution. A more elegant and efficient solution present in many research [START_REF] Khodabandeh | A robust learning approach to domain adaptive object detection[END_REF][START_REF] Geirhos | Imagenettrained cnns are biased towards texture; increasing shape bias improves accuracy and robustness[END_REF][START_REF] Cygert | Robustness in compressed neural networks for object detection[END_REF][START_REF] Poursaeed | Robustness and generalization via generative adversarial training[END_REF] is to enrich the database used for learning by simulating the most relevant distortions. It is in this last framework that the proposal described in this 1 https://github.com/Aymanbegh/Benchmarking-performance 2 https://github.com/Aymanbegh/Distorted-Natural-COCO work is made. Thereby, we propose a full framework evaluation of robustness for a set of object detection methods [START_REF] Bochkovskiy | Yolov4: Optimal speed and accuracy of object detection[END_REF][START_REF] Tan | Efficientdet: Scalable and efficient object detection[END_REF][START_REF] He | Mask r-cnn[END_REF] through several distortions applied on the MS-COCO dataset [START_REF] Lin | Microsoft coco: Common objects in context[END_REF]. This data augmentation is performed through some common global distortions (noise, motion blur, defocus blur, haze, rain, contrast change, and compression artefacts) in the whole image and some local distortions (local blur motion, BackLight illumination and local defocus blur) in specific areas that include the possible dynamic objects or scene conditions. To the best of our knowledge, our work is the first to consider the local distortions to assess the robustness of object detection models. Finally, we proceeded to a manual selection of the distorted images present in the validation set of the MS-COCO dataset and retained only the annotations of the distorted objects to perform a more relevant evaluation of the impact of distortions on models. The main contributions of our study are summarized in the following:

• A comprehensive evaluation of the robustness of the stateof-the-art object detection methods against global and local distortions at 10 levels of distortion is provided (see section 4.1).

• A dataset dedicated to the study of the impact of local and global distortions on the robustness of object detection is built from MS-COCO dataset [START_REF] Beghdadi | Benchmarking performance of object detection[END_REF].

• It is shown that the robustness of object detection process is improved by enriching the training process by using data augmentation based on generated distorted images (see table 3).

• An additional new benchmark dataset constructed for natural in-capture distortions from MS-COCO dataset that allows to better assess the reliability of the object detection models in case of real distortions (see tables 2 and 4) [START_REF] Beghdadi | Natural distorted ms-coco dataset for object detection[END_REF].

The remainder of the paper is organized as follows. Section 2 summarizes previous related literature. Section 3 is devoted to detail the methods of evaluation and the distortion generations. Then, section 4 is dedicated to show and discuss results. Finally, concluding remarks and perspectives are provided in section 5.

RELATED WORK

The evaluation of object detection methods from images acquired in an uncontrolled environment has been the subject of a few studies [START_REF] Cygert | Toward robust pedestrian detection with data augmentation[END_REF][START_REF] Chen | Robust and accurate object detection via adversarial learning[END_REF]. But most of the time, the databases used contains images with limited number of degradation types and therefore cannot be generalized to practical cases where one has to face several distortions. Some papers deal with the impact of various distortions on the detection performance such as in the study [START_REF] Azulay | Why do deep convolutional networks generalize so poorly to small image transformations?[END_REF] that attempts to evaluate the robustness of some backbones of CNN architectures against geometric transformations. A comparative study between DNN-based architectures (ResNet-152, VGG-19, GoogLeNet) and human observer performance was conducted in [START_REF] Geirhos | Generalisation in humans and deep neural networks[END_REF]. In this study twelve common distortions are considered in the performance evaluation of object recognition methods. It has been shown that these DNN-based architectures perform remarkably well in the case of the distortions on which the training was done. However, this performance decreases, compared to that of the human observer, for unseen distortions. Another interesting study on the evaluation of the robustness of classification models of images from the imagenet database and subjected to fifteen artificial distortions at five levels of severity, revealed the effect of image quality degradation on the stability of the training process for various DNN-based architectures [START_REF] Hendrycks | Benchmarking neural network robustness to common corruptions and perturbations[END_REF]. Recently, a benchmark study on the robustness of many object detection models against 15 types of global distortions and stylized images through data augmentation process has been done in [START_REF] Michaelis | Benchmarking robustness in object detection: Autonomous driving when winter is coming[END_REF].

PROPOSED BENCHMARKING FRAMEWORK

In order to evaluate the robustness of the considered models against image distortions in an uncontrolled environment, we propose 10 types of distortions (global and local) at 10 linearly increasing levels of severity [START_REF] Beghdadi | Benchmarking performance of object detection[END_REF]. In figures 1 (b), (c), and (d), distortion levels are 8, 7, and 10 respectively.

Image distortions

Global distortions affect the image as a whole and come from different sources related in general to the acquisition conditions. Some are directly dependent on the physical characteristics of the camera and are of photometric or geometric origin. Among the most common distortions that affect the quality of the signal are defocus blur (Defocus-Blur), photon noise, geometric or chromatic aberrations, and blur (Motion-Blur) due to the movement of the camera. The other types of degradation are related to the environment and, more particularly, the lighting and atmospheric disturbances in the case of outdoor scenes. Compression and image transmission artifacts are another source of degradation that is difficult to control. These common distortions have been already considered in benchmarking the performance of some models [START_REF] Michaelis | Benchmarking robustness in object detection: Autonomous driving when winter is coming[END_REF][START_REF] Hendrycks | Benchmarking neural network robustness to common corruptions and perturbations[END_REF][START_REF] Benz | Revisiting batch normalization for improving corruption robustness[END_REF].

Local distortions are undesirable signals affecting one or more localized areas in the image (see figure 1). A typical case is the Motion Blur (Loc. MBlur) due to the movement of an object of relatively high speed. Another photometric distortion is the appearance of a halo around the object contours due to the limited sensitivity of the sensors or backlight illumination (BackLight). The artistic blur (Loc. Defoc.) affecting a particular part of the targeted scene, the object to be highlighted by the pro-shooter, is another type of local distortion. Thus, integrating the local distortions in the database increases its size and makes it richer and more representative of scenarios close to real applications. It is worth noticing that these local distortions are made possible thanks to the annotations provided in the MS-COCO database and, in particular, the details of the shape and location of the object of interest.

Description of the considered models

In this study we focus on few models of state-of-the-art object detection methods which have the best performance or notoriety, namely YOLOv4 [START_REF] Bochkovskiy | Yolov4: Optimal speed and accuracy of object detection[END_REF], Mask R-CNN [START_REF] He | Mask r-cnn[END_REF] and EfficientDet [START_REF] Tan | Efficientdet: Scalable and efficient object detection[END_REF] with distinct backbone architectures that highlights the correlation between the backbone and their robustness. The table 1 summarizes the different models and their components which are used in this benchmarking.

Robustness Evaluation against Natural Distortions

We manually identified the natural distortions present in the MS-COCO validation set to better assess the robustness in real scenarios. Indeed, in the existing methods, robustness evaluations were performed on synthetic distortions or specific datasets, limiting those studies' extension in real-world applications. It is worth noticing that the selected sub-sets from MS-COCO contain images with global distortions and the associated object annotations. In contrast, only affected objects are considered in the case of local distortions (see figure 2). This benchmarking framework for robustness evaluation uses one of the most widespread databases for object detection, restricted to a limited set of seven natural distortions types (see table 2). In this way, we are laying the foundations for a future common evaluation method for evaluating object detection robustness against real-world distortions. It is important to notice that some natural distortions such as Defocus, Object Blur, and BackLight have low ratios (see table 2). This results from the nature of local distortions, which only affect specific objects in the images (see figure 2).

RESULTS AND DISCUSSIONS

Evaluation of the robustness

We evaluated the robustness of models against distortions on the MS-COCO 2017 validation set, which contains 5K annotated images that we have corrupted through our 10 types of distortions [START_REF] Beghdadi | Benchmarking performance of object detection[END_REF]. Note that level 0 of distortion represents the COCO AP score for non-distorted images. We assessed the object detection performance of models through the AP (Average Precision) from the COCO metrics that provides a relevant performance indicator. Figures 3 and4 Fig. 4: Robustness of object detection models against distortions.

represent the evolution of the AP score of models according to each distortion type and level. These results show that distortions reduce object detection performance of models by 20% to 89.2%. Moreover, they highlight the weaknesses of YOLOv4 methods against distortions (strongly decreasing curves in figures 4 and 3), which need a deeper study. The evaluations in figure 3 indicated that global distortions have more impact than local distortions on the models. While local and atmospheric disturbances have a roughly uniform impact on methods despite their different architectures, the other distortions produce effects separate and independent. Indeed, overall models performance are greatly affected by distorted images but not uniformly. Each type of distortion has a different rate of impact on models performance depending on the distortion level and object detection method. In order to evaluate the impact of distortions on models relevantly, we computed the robustness rate that represents the ratio between mAP scores of models through each distortion type and level, and non-distorted images (level 0). This robustness rate is illustrated (see figure 5) with a data distribution visualization including boxplots and violin plots that contribute to extract richer information such as the mean, median, and probability density. The diversity of architectures with different network depths made it possible to estimate a possible correlation between this depth and the robustness of the models. Results from figures 3, 4 and 5 indicate that the EfficientDet method's robustness decreases despite its network depth and input size increasing through its different backbone architectures (D0 to D3). Furthermore, the YOLOv4-tiny model was the least robust to distortions with an average decrease in its AP score of */-: Mean/Median of mAP for all distortion type and level, Light/dark surrounding shape: probability density/percentiles. 36.7% between the distortion levels 0 and 10 (lowest average robustness rate at 63.3%, see mAP score in figure 5). Therefore, it seemed relevant to us to choose this model to evaluate the improvement of robustness through a data augmentation with our distorted dataset.

Training on distorted dataset

This section presents experimental results on MS-COCO 2017 evaluation set with 5K images. YOLOv4-tiny model is trained following the protocol in MS-COCO challenge, i.e., using the trainval35k set for training that includes 80k images from the train set and 35k images from the validation set. We first performed this training with the original MS-COCO trainval35k set to obtain reference values, then training with the distorted dataset, which allowed us to evaluate the contribution of the data augmentation on the robustness of the model. Each distortions type represents 5% of the database used for the training process, i.e., 5.9K images for each distortion. In our case, the hyper-parameters are taken from the YOLOv4 method but with some tuning. The training step is 500,050, with the initial learning rate multiplied by 0.1 at 400,040 step then 450,045 step. Evaluation of the training improvement is achieved by comparing the evaluation scores of both models trained on original and distorted train sets, respectively, on the distorted MS-COCO validation sets according to each distortion level. These results are formulated in the table 3 (Columns 2-6) by computing the ratio between AP COCO scores of the distorted model on the natural model. We observed that training has a significant impact on the robustness improvement with an average increase of 31.5% of performance (see last column in table 3). Furthermore, the study highlighted that the training improves more robustness for high-level distortions (see last row in table 3). According to the results of the table 4, training has a weak impact on robustness to natural distortions with an average improvement of 2.5% for the mAP score. However, these results need perspective due to the distortion levels in the distorted natural datasets [START_REF] Beghdadi | Natural distorted ms-coco dataset for object detection[END_REF], which are primarily low, and, therefore, not impacted enough by data augmentation from training according to the previously highlighted observation (see average per level of column 2 in table 3).

CONCLUSIONS AND FUTURE WORK

This study revealed the impact of the quality of the image frames in the development of object detection methods. Indeed, through this study we have shown that increasing the training database with complex scenarios containing different distortions improves the performance of object detection models. We have also shown that this is also true even for the best state-of-the-art methods. We, particularly, propose an original strategy of image distortion generation applied to the MS-COCO dataset that combines some local and global distortions to reach a better realism. We have shown that training with this distorted dataset improves the robustness of models by up to 31,5%. We, also, provide a custom dataset including the natural images distorted from MS-COCO dataset to perform a more relevant evaluation of the robustness concerning distortions. The database and the generation source codes of the different distortions is made publicly available on github plateform. Considering dataset with local and global distortions is an interesting way for future research.

Fig. 1 :

 1 Fig. 1: Illustration of some local distortions.

Fig. 2 :

 2 Fig. 2: Natural distortions selection in the MS-COCO dataset.

Fig. 3 :

 3 Fig. 3: Robustness of object detection models against global(left)/local (right) distortions.

Fig. 5 :

 5 Fig. 5: Average robustness rate of object detection models against the distortions.

Table 1 :

 1 State-of-the-art object detection models.

	Methods	Architecture	Size Add-blocks
	Mask R-CNN [15] Resnet-101	512	FPN, RPN [23, 24]
	YOLOv4 [13]	CSPDarknet53 YOLOv4-tiny	512 PAN [25] 416
		EfficientNet-B0	512
		EfficientNet-B1	640
	EfficientDet [14]	EfficientNet-B2	768	Bi-FPN [14]
		EfficientNet-B3	896
		EfficientNet-B4	1024

FPN: Feature Pyramid Network, RPN: Region Proposal Network, PAN: Path Aggregation Network, Bi-FPN: Bi-directional FPN.

Table 2 :

 2 Features of distorted natural sub-sets.

		Noise Contrast Blur Defoc. Rain O.Blur BackLight
	Images 44	42	32 201	21	127	128
	Objects 289 312	224 1299 110 464	934
	Ratio	1.0 0.83	0.97 0.72 0.95 0.34 0.68

Ratio: Number of retained objects/annotated objects in images.

Table 3 :

 3 Object detection percentage performance of the YOLOv4-tiny trained model on our distorted dataset.

	Distortion		Distortion level (%)		mAP
	Type	2	4	6	8	10	Average
	Noise	9.58	22.3	37.9 80.3 114	47.2%
	Contrast	-0.54 1.8	2.47 3.23 5.71	1.99%
	Compression -2.36 -0.48 -0.5	4.84 35.4	4.26%
	Rain	12.7	40.4	71.4 101.8 114.3 61.7%
	Haze	-0.49 3.14	10.5 29.8 51.7	15.8%
	Motion-Blur 6.11	39.5	95.3 151.4 183.3 86.3%
	Defoc-Blur	1.02	14.2	22.4 26.3 26.7	16.5%
	Loc. MBlur	10.3	31.6	46.9 59.3 67.0	39.8%
	Loc. Defoc.	3.08	16.0	23.3 25.5 25.9	17.2%
	BackLight	1.45	5.85	19.1 40.3 76.8	24.1%
	Average per level	4.11% 17.9% 32.9% 52.3% 70.1% 31.5%

Table 4 :

 4 Object detection performance of the YOLOv4-tiny trained model on our distorted natural sub-sets.

	Natural distortions	AP	rel. AP	mIoU	rel.mIoU
	Noise	0.003	0.99%	0.005	0.65%
	Contrast	0.011	5.47%	0.007	0.95%
	Motion-Blur	-0.012	-5.17%	-0.004	-0.54%
	Defocus	0.015	12.9%	0.0	0.0%
	Rain	-0.009	-2.81%	0.006	0.79%
	Loc. MBlur	-0.006	-3.13%	0.0	0.0%
	BackLight	0.011	8.94%	0.005	0.68%
	rel.: ratio between scores from distorted and original images, mIoU:
	mean Intersection over Union.