*8.3. Case λ<sup>p</sup>* > 1 *and λ*<sup>1</sup> < 1

This case guarantees that 0 ≤ 1 − *λj*/*λ<sup>p</sup>* < 1, *j* = 1, ... , *p* − 1 and 0 ≤ 1 − 1/*λ<sup>p</sup>* < 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 *<sup>∂</sup> <sup>∂</sup><sup>a</sup>* {*F*(*p*) *<sup>D</sup>* (.)}|*a*=<sup>0</sup> is possible. The following case where *λ*<sup>1</sup> = ... = *λ<sup>p</sup>* = *λ* allows the Lauricella series to be equivalent to the Gauss hypergeometric function given as follows

$$F\_D^{(p)}\left(a,\underbrace{1\ldots1}\_{p};\underbrace{1\ldots1}\_{p};a+\frac{1+p}{2};1-\lambda,\ldots,1-\lambda\right) = \,\_2F\_1\left(a,\underbrace{p}\_{2};a+\frac{1+p}{2};1-\lambda\right).\tag{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 *<sup>F</sup>*(*p*) *<sup>D</sup>* (.) and <sup>2</sup>*F*1(.), along with the absolute value of error ||, where

*<sup>p</sup>* <sup>=</sup> 2, *<sup>N</sup>* <sup>=</sup> {20, 30, <sup>40</sup>}. The exact expression of *<sup>∂</sup> <sup>∂</sup><sup>a</sup>* {2*F*1(.)}|*a*=<sup>0</sup> 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.

**Table 2.** Computation of A = *<sup>∂</sup> <sup>∂</sup><sup>a</sup>* {2*F*1(.)}|*a*=<sup>0</sup> and <sup>B</sup> <sup>=</sup> *<sup>∂</sup> <sup>∂</sup><sup>a</sup>* {*F*(*p*) *<sup>D</sup>* (.)}|*a*=<sup>0</sup> when *p* = 2 and *λ*<sup>1</sup> = ... = *λ<sup>p</sup>* = *λ*.


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 **Σ**<sup>1</sup> and **Σ**<sup>2</sup> 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.

**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*<sup>3</sup> + 1.5 × *IQR* or smaller than *Q*<sup>1</sup> − 1.5 × *IQR*, where *Q*1, *Q*3, and *IQR* are the 25th, 75th percentiles, and the interquartile range, respectively.

**Table 3.** Parameters **Σ**<sup>1</sup> and **Σ**<sup>2</sup> used to compute KLD for central MCD.


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.

```
%********************************************************************************
% Compute the KL dive rgence and dis t ance between two ce n t r al mul tiv a ri a te Cauchy
% distribution .
%% Inpu t :
% + Sigma1 : Symmetric positive definite (p*p) scale matrix
% + Sigma2 : Symmetric positive definite (p*p) scale matrix
% + nb : indi ce s used to compute the KL and di s ; nb = { 2 0 , 3 0 , 4 0 , e t c } .
% Inc re a se nb means inc re a se the p reci sion and also the computation time .
%% Output :
% + KL_12 : KL dive rgence between X1 and X2 : KL ( X1||X2 )
% + KL_21 : KL dive rgence between X2 and X1 : KL ( X2||X1 )
% + Esp_12 : expe c t a ti o n E_ {X } \ { ln [1+ X^T* Sigma2^{ −1}*X ] } } where X~MCD( Sigma1 , p=3)
% + Esp_21 : expe c t a ti o n E_ {X } \ { ln [1+ X^T* Sigma1^{ −1}*X ] } } where X~MCD( Sigma2 , p=3)
% + di s : di s t a n c e between X1 and X2 : di s = KL ( X1||X2 ) + KL ( X2||X1 )
%% Example :
% Sigma1 = [1 0.6 0.2; 0.6 1 0.3; 0.2 0.3 1];
% Sigma2 = [1 0.3 0.1; 0.3 1 0.4; 0.1 0.4 1];
% [ KL_12 , KL_21 , Esp_12 , Esp_21 , dis ] = fonction_KL_MCD_final ( Sigma1 , Sigma2 , 2 0 ) ;
%********************************************************************************
function [ KL_12 , KL_21 , Esp_12 , Esp_21 , dis ] = fonction_KL_MCD_final ( Sigma1 , Sigma2 , nb )
format long ;
p = 3;
vpr = real ( eig ( Sigma1 * inv ( Sigma2 ) ) ) ;
vpr = sort ( vpr , ' ascend ' ) ;
nbre = nb ;
[N,M, L ] = ndgrid ( 0 : nbre , 0 : nbre , 0 : nbre ) ;
i f vpr (p)< 1
%**************************************************************************
%% D e ri v a ti v e o f Fd ( a , 1 / 2 , 1 / 2 , 1 / 2; a+(1+p)/2;1 −vpr (1) ,1 − vpr (2) ,1 − vpr ( 3 ) )| a=0
%**************************************************************************
  H = N+M+L ;
  H(H==0)= i n f ;
  commun = (1 − vpr ( 1 ) ) . ^N./ f a c t o r i a l (N) .*( 1 − vpr ( 2 ) ) . ^M./ f a c t o r i a l (M) .*(1 − vpr ( 3 ) ).^L./ f a c t o ri al (L ) . * ...
      pochhammer ( 1/ 2 ,N ) . * pochhammer ( 1/ 2 ,M) . * pochhammer ( 1/ 2 ,L ) ;
  h1 = commun .* pochhammer ( 1 ,N+M+L ) . / pochhammer ( ( 1 + p )/ 2 , N+M+L )* 1./H;
  derive1 = sum(sum(sum( h1 ) ) ) % Eq . ( 1 0 2 ) and (A1 )
  %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
  J = N+M+L− 1;
  A= [0 , cumsum(1./((p+1)/2 +(0:p* nbre −1)))];
  for i =1: nbre+1
       for j =1: nbre+1
            for l =1: nbre+1
              G( i , j , l )= −A( J ( i , j , l ) + 2 ) ;
            end
       end
  end
  h2 = commun .*G;
  derive2 = sum(sum( sum( h2 ) ) ) % Eq . ( 1 0 4 ) and (A1 )
  %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
  Esp_12 = psi (1/2 + p/2)−psi (1/2) − derive1
  Esp_21 = psi (1/2 + p/2)−psi (1/2) − prod ( vpr ).^ (1/2 )* derive2
  KL_12 = −1/2*log ( prod ( vpr )) −(1+p)/2* derive1 % Eq . ( 1 0 9 )
```

```
KL_21 = 1/2*log ( prod ( vpr )) −(1+p)/2* prod ( vpr ).^ (1/2 )* derive2 % Eq . ( 1 1 5 )
elseif vpr (1 ) > 1
%*************************************************************************************
%% De ri v a ti ve o f Fd ( ( 1 + p )/ 2 , 1/ 2 , 1/ 2 , 1/ 2; a+(1+p)/2;1−1/vpr (1) ,1 −1/vpr (2) ,1 −1/vpr ( 3 ) )| a=0
%*************************************************************************************
 J = N+M+L− 1;
 A= [0 , cumsum(1./((p+1)/2 +(0:p* nbre −1)))];
 for i =1: nbre+1
      for j =1: nbre+1
           for l =1: nbre+1
              G( i , j , l )= −A( J ( i , j , l ) + 2 ) ;
           end
      end
 end
 commun = (1 −1/ vpr ( 1 ) ) . ^N./ f a c t o r i a l (N) .*( 1 − 1 / vpr ( 2 ) ) . ^M./ f a c t o r i a l (M) .*(1 −1/ vpr ( 3 ) ).^L./ f a c t o ri al (L ) . * ...
      pochhammer ( 1/ 2 ,N ) . * pochhammer ( 1/ 2 ,M) . * pochhammer ( 1/ 2 ,L ) ;
 h1 = commun .*G;
 derive1 = sum( sum(sum( h1 ) ) ) % Eq . ( 1 0 4 ) and (A1 )
 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
 H = N+M+L ;
 H(H==0)= i n f ;
 h2 = commun .* pochhammer ( 1 ,N+M+L ) . / pochhammer ( ( 1 + p )/ 2 , N+M+L )* 1./H;
 derive2 = sum( sum(sum( h2 ) ) ) % Eq . ( 1 0 2 ) and (A1 )
 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
 Esp_12 = psi (1/2 + p/2)−psi (1/2) − prod ( vpr ).^( −1/2)* derive1
 Esp_21 = psi (1/2 + p/2)−psi (1/2) − derive2
 KL_12 = −1/2*log ( prod ( vpr )) −(1+p)/2* prod ( vpr ).^( −1/2)* derive1 % Eq . ( 1 1 0 )
 KL_21 = 1/2*log ( prod ( vpr )) −(1+p)/2* derive2 % Eq . ( 1 1 4 )
else
%******************************************************************************************
%% D e ri v a ti v e o f Fd ( a , 1/ 2 , 1/ 2 , a +1/2; a+(1+p)/2;1 −vpr (1)/ vpr (3) ,1 − vpr (2)/ vpr (3) ,1 −1/vpr ( 3 ) )| a=0
%*******************************************************************************************
 H = N+M+L ;
 H(H==0)= i n f ;
 commun = (1 − vpr ( 1 ) / vpr ( 3 ) ) . ^N./ f a c t o r i a l (N) .*( 1 − vpr ( 2 ) / vpr ( 3 ) ) . ^M./ f a c t o r i a l (M ) . * ...
            (1−1/vpr ( 3 ) ).^L./ f a c t o ri al (L ) . * ...
            pochhammer ( 1/ 2 ,N ) . * pochhammer ( 1/ 2 ,M) . * pochhammer ( 1/ 2 ,L ) ;
 h1 = commun .* pochhammer ( 1 ,N+M+L ) . / pochhammer ( ( 1 + p )/ 2 ,N+M+L ) . /H;
 derive1 = sum(sum( sum( h1 ) ) ) % Eq . ( 1 0 1 ) and (A1 )
 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
 J = N+M+L− 1;
 JJ = L−1;
 A = [0 , cumsum(1./((p+1)/2 +(0:p* nbre −1)))];
 B = [0 , cumsum(1./(1/2 +(0: nbre −1)))];
 for i =1: nbre+1
      for j =1: nbre+1
           for l =1: nbre+1
                G( i , j , l )= B ( J J ( i , j , l )+2) − A( J ( i , j , l )+ 2 );
           end
      end
 end
 h2 = commun .*G;
 derive2 = sum( sum( sum( h2 ) ) ) % Eq . ( 1 0 3 ) and (A1 )
 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
 Esp_12 = psi (1/2 + p/2)−psi (1/2) + log ( vpr (p ) ) − derive1
 Esp_21 = psi (1/2 + p/2)−psi (1/2) − vpr (p)^(−p/2)*( vpr ( 1 )* vpr (2 ) )^ (1/2 )* derive2
 KL_12 = −1/2*log ( prod ( vpr ) ) −(1+p)/2*( − log ( vpr (p ) ) + derive1 ) % Eq . ( 1 0 6 )
 KL_21 = 1/2*log ( prod ( vpr ) ) −(1+p)/2* vpr (p)^(−p/2)*( vpr ( 1 )* vpr (2 ) )^ (1/2 )* derive2 % Eq . ( 1 1 2 )
end
```
dis = KL\_12 + KL\_21

#### **9. Conclusions**

Since the MCDs have various applications in signal and image processing, the KLD between central MCDs tackles an important problem for future work on statistics, machine learning and other related fields in computer science. In this paper, we derived a closedform expression of the KLD and distance between two central MCDs. The similarity measure can be expressed as function of the Lauricella D-hypergeometric series *<sup>F</sup>*(*p*) *<sup>D</sup>* . 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.

**Author Contributions:** Conceptualization, N.B.; methodology, N.B.; software, N.B.; writing original draft preparation, N.B.; writing review and editing, N.B. and D.R.; supervision, D.R. All authors have read and agreed to the published version of the manuscript.

**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.

**Conflicts of Interest:** The authors declare no conflict of interest.
