ARXIV PREPUBLICATION

1

Regularized maximum likelihood estimation of
covariance matrices of elliptical distributions

arXiv:1611.10266v1 [math.ST] 30 Nov 2016

Christophe Culan , Claude Adnet∗ ,

Advanced Radar Concepts division, Thales, Limours 91338, France

Abstract—The maximum likelihood principle is widely used in
statistics, and the associated estimators often display good properties. indeed maximum likelihood estimators are guaranteed to
be asymptotically efficient under mild conditions.
However in some settings, one has too few samples to get a good
estimation. It then becomes desirable to take into account prior
information about the distribution one wants to estimate.
One possible approach is to extend the maximum likelihood
principle in a bayesian context, which then becomes a maximum
a posteriori estimate; however this requires a distribution model
on the distribution parameters.
We shall therefore concentrate on the alternative approach of
regularized estimators in this paper; we will show how they
can be naturally introduced in the framework of maximum
likelihood estimation, and how they can be extended to form
robust estimators which can reject outliers.

A. Notations and conventions
In the following development, the following notations and
conventions shall be observed:
•

•

Index Terms—Maximum likelihood estimation, KL-divergence,
entropy, elliptical distributions, complex elliptical distributions,
adaptive detection, iterative algorithm, outlier detection, outlier
rejection, regularized estimators

I. I NTRODUCTION : REGULARIZED ESTIMATION
Many statistical estimators are shown to be asymptotically
well behaved; in particular maximum likelihood estimators are
consistent and efficient under relatively mild conditions, and
often have good convergence properties, thus making them
a popular choice for parametric estimation [21][25][12][1].
However in a few sample regime, their performances can be
greatly degraded, and in some cases one may not even be able
to compute these estimators [20]; this is notably the case for
covariance estimations when the number of samples is of the
same order or lower than the dimension of the space [13][29].
In such situations, one has to add additional information
in the estimation procedure, other than available samples.
Hence a popular approach is to use regularized estimators,
which are estimators purposely biased towards a prior distribution representing an initial knowledge of the model; the
effect generally is a smoothing of the estimated parameter,
although this largely depends upon the shape of the prior
distribution; moreover this prior is not necessarily explicitly
known. Various regularization methods have been proposed
[30][4][2][9][16][6][11]. A popular approach is to derive
such estimates by a penalized maximum likelihood approach
[31][22][29][17].
The applied penalties are often chosen for convenience of
computation; we introduce in this article a way to define such
penalties which are in agreement with the theory of maximum
likelihood.
Manuscript received October XX, 2016; revised October XX, 2016. Corresponding author: M. Culan (email: christophe.culan@thalesgroup.com).

•
•

•
•
•
•

•

For any topological space R1 and R2 , the topological
space R1 × R2 is the product of R1 and R2 , whereas
R1 ∨ R2 is the disjoint union of R1 and R2 . One also
notes, for a topological space R, RN to be the product
of N copies of R and ∨N R to be the disjoint union of
N copies of R.
Let R1 , ..., RN be topological spaces, and let µn be
a measure defined onQ Rn for 1 ≤ n ≤ N . The
N
measure µ1 . . . µn = n=1 µn is the product measure
QN
of (µn )1≤n≤N defined on n=1 Rn , whereas
WN
WN
µ1 ∨ · · · ∨ µn = n=1 µn is the measure of n=1 Rn
such that its restriction to Rn is µn for 1 ≤ n ≤ N and
is the joint measure of (µn )1≤n≤N .
Moreover for a topological space R and a measure µ
defined on R, µN is the product measure of N copies of
µ defined on RN , whereas ∨N µ is the joint measure of
N copies of µ defined on N · R.
The vector space Cd is canonically identified to R2d [23].
Hd (R) is the set of symmetric matrices of size (d, d);
Hd (C) is the set of hermitian matrices of size (d, d).
Hd+ (K) is the subset of matrices of Hd+ (K) which are
+
positive; HP +
d (K) is the subset of matrices of Hd (K)
of unit determinant.
The vector space Cd is canonically identified to R2d [23].
X † is the transpose conjugate of any matrix or vector X.
µRd is the Lebesgues measure of the space Rd . Similarly,
µCd is the Lebesgues measure of Cd ' R2d .
Sd−1 is the (d − 1)-sphere, which is identified as the following part of Rd : Sd−1 ' x ∈ Rd ; x† x = 1 . Similarly
S2d−1 is identified to x ∈ Cd ; x† x = 1 in Cd .
sd−1 shall denote the probability distribution of Sd−1
isotropic for the canonical scalar product of Cd , defined
by:
d−1
X
hy, xi =
yk xk
(1)
k=0

R
• δx
•

is the Dirac delta distribution centered in x in space
R.
One shall note:
 R 
Z

dδ0
(t) dδ0R (t)
(2)
hδ = H δ0R |µR =
log
dt
t∈R
It is the entropy of the Dirac distribution relative to
the Lebesgues measure of R due to its translational

ARXIV PREPUBLICATION

2

symmetry, which is positive and infinite. It can be used
to express different other entropies in higher dimensional
settings:
 d

H δxR |µRd = dhδ
(3)

 d
H δxC |µCd = 2dhδ
II. M AXIMUM LIKELIHOOD WITH PRIOR INFORMATION
A. General theory
Let us first recall that likelihood can be understood in terms
of information; indeed it can be expressed as the relative
entropy (or KL-divergence [15]) of the sampling distribution
S relative to the model distribution P [7]:


Z
l(P |S) = −H(S|P ) =

log
x∈R


dP
(x) dS(x)
dS

(4)

Under the usual hypothesis in which one has N i.i.d samples
(xn )1≤n≤N , the usual sampling distribution is defined on a
space ∨N B, with B being some base space, and is given by the
disjoint union of Dirac distributions centered on each sample
[7]:
N
1 _ B
S=
δ
N n=1 xn

The model distribution is given by
being some distribution of B.

(5)
1
N

∨N PB , with PB

Thus the likelihood can be expressed under this hypothesis
by:
l(P |S) = −

N
1 X
H(δxBn |PB )
N n=1

(6)

Obviously this particular instantiation of the maximum
likelihood problem is not well suited for our current problem,
as we would like to take into account some prior information
in addition to the information provided by samples.
We shall focus on a situation in which we have an existing
prior estimate, and we want to take into account this prior in
our estimation procedure [29][17].
Fortunately the entropic formulation of the maximum
likelihood is quite flexible; in particular one can use different
models for the sampling distribution, as is done for example
in [7] to take into account the circular symmetry of the
distributions.
We therefore propose to take into account the prior knowledge
as part of the sampling distribution. The sampling distribution
under the i.i.d hypothesis, with additional prior information
can then be modeled as:

S=
with:

N
α _
Sn
N n=1

!
∨ ((1 − α)Pprior )

(7)

Sn the sample distribution associated to the n-th sample,
defined on a base space B; it would generally be given
by δxBn .
• Pprior the prior distribution. It is not necessarily defined
on the same space as the sample distributions depending
on the information it carries, and might be defined on a
quotiented image of the base space, which shall be noted
as [B].
• α ∈ [0; 1] the integration factor, which represents how
much ones trusts the prior estimate compared to the newly
acquired samples.
The distribution model one wants to fit can on the other
hand be expressed as:

α
∨N P ∨ ((1 − α)[P ])
(8)
N
with:
• P the distribution model
• [P ] the distribution model resulting from P on the
quotiented space [B]
The corresponding likelihood can therefore be expressed as:
•

N
α X
H(Sn |P ) − (1 − α)H(Pprior |[P ])
l(P |S) = −
N n=1

(9)

Therefore this results in a penalized version of maximum
likelihood , the penalty being defined from the relative entropy
of the prior distribution for the model distribution [17][29].
III. R EGULARIZED COVARIANCE ESTIMATORS
We shall now apply the regularized maximum likelihood
equation to different covariance estimation problems. We’ll
focus on the two following cases:
• The distribution which we try to estimate is elliptically
symmetric (ES/CES) [5][18]. No hypothesis is made on
the radial distribution of this distribution; this corresponds
to the conditions for which Tyler’s estimator is the maximum likelihood estimator for i.i.d samples [7][27][19].
• The distribution is gaussian, or circular gaussian for complex distributions. Note that the corresponding maximum
likelihood estimators under the i.i.d hypothesis are given
in [7].
A. Regularized Tyler’s estimator
Let us first recall that an elliptical distribution can be
expressed in general by [7]:
√

1
dP (x) = p dQ
x† R−1 x dsd−1
|R|



L(R)−1 x
√
x† R−1 x


(10)

with the radial distribution Q being any probability measure
on R+ , and the correlation matrix R being any matrix of
HP +
d (K) [7].
Supposing here that one has no a priori on the radial
distribution, the prior distribution one wants to use should
be defined on the unit sphere, and corresponds to the unique

ARXIV PREPUBLICATION

3

probability distribution of the unit sphere which is isotropic
Rprior
for the corresponding correlation matrix Rprior , noted as sd−1
.
The entropy between two such distributions isotropic for
correlation matrices R1 and R2 shall be noted H(R1 |R2 ).
Thus the concentrated likelihood for samples (xn )1≤n≤N in
Rd under i.i.d hypothesis is given, after maximization over
the radial distribution, by [7]:

l(R|S) = − dhδ −

N

α(d − 1) X
log xn † R−1 xn
2N
n=1

1:

function REG T YLER(Rprior ,(xn )1≤n≤N ,α,)

2:

R ← Rprior

3:

repeat

√

R 12 ←

4:
5:

R− 12 ← R 12 −1

6:



S ← (1 − α)anscm R− 21 Rprior R− 21

(11)

+

− (1 − α)H(Rprior |R)
H(Rprior |R) however represents a challenge to compute.
Fortunately, one can characterize its differential with respect
to R (see A for more details):



1
1
1
d − 1  1
tr R 2 anscm(R− 2 Rprior R− 2 )R 2 d R−1
2
(12)

d
tr(S) S

S←

8:

R ← R 21 exp(S − I)R 12
R
R ← tr(R)


2
until tr (R−1 R − I) ≤ 

10:

=

N
X
α
xn xn †
R− 21
R 1
†
N
x R− 21 R− 21 xn − 2
n=1 n

7:

9:

dH(Rprior |R)

R

return R− 12

11:
12:

end function

In the case of a distribution in Cd under the circularity
hypothesis, the maximum likelihood equation under i.i.d hypothesis is slightly different; indeed the sampling distribution
is then given by:

with anscm(Σ) being the expected value of the normalized
sample covariance estimator (NSCM) for an elliptical distribution of correlation matrix proportional to Σ.
Unfortunately there is no explicit way of computing this matrix
in the general case to the authors’ knowledge and one has
to resort to numerical methods to extract its eigenvalues.
However it has been characterized for complex circularly
symmetric distribution in [3].

By deriving the concentrated likelihood in a manner similar
to [7] and differentiating it, on gets the following equation:

In any case this leads to the following maximum likelihood
equation:


tr λR

S=

−

Without any structural constrain on R, this simplifies to:


 1
1
1
1
1
R=
(1 − α)R 2 anscm R− 2 Rprior R− 2 R 2
2λ(d − 1)

N
α X xn xn †
+
N n=1 xn † R−1 xn
(14)
This can be solved by using the following numerical
procedure:

(15)




 1
1
1
1
1
− (1 − α) d −
R 2 anscm R− 2 Rprior R− 2 R 2
2


N

α(d − 1) X xn xn †
−1
−
d R
=0
N
x † R−1 xn
n=1 n
(16)


tr λR

 1
1
1
(1 − α)(d − 1) 1
R 2 anscm R− 2 Rprior R− 2 R 2
2


N

α(d − 1) X xn xn †
−1
−
d
R
=0
2N
x R−1 xn
n=1 n
(13)

N
α _ C
Rprior
δ
+ (1 − α)sd−1
N n=1 [xn ]

This is almost the same equation, although each sample
actually carries less information compared to the prior
distribution for a same α parameter.
The equation can thus be solved numerically by calling the
same algorithm REG T YLER with a different α parameter:
REG

T YLER(Rprior ,(xn )1≤n≤N , d−α(d−1)
,)
1
(1+α)
2

Note however that the computation of the ansm function
differs in the real and complex circular cases (see annex).
These regularized Tyler estimators are more complicated
to compute than those found in the existing litterature by
applying different penalties (specifically ones maintaining the
convexity of the problem) and offer no theoretical guaranty

ARXIV PREPUBLICATION

4

of convergence [28][17][26]. However they use a penalty
which is consistent with the definition of the likelihood as an
information theoretic quantity, which is the aim of this paper.
B. Regularized estimators for gaussian distributions
A subclass of interest of elliptical models are estimations
constrained on a subclass of models corresponding to a fixed
scaled radial distribution, of the from:


L(R)−1 x
(17)
dP (x) = σdQ0 (σr)dsd−1 √
x† R−1 x
Indeed most M-estimators can be expressed as maximum
likelihood estimators for one of these subclasses of models
[7].
However the existence of an easy to compute regularized
estimator then depends on the expression of the gradient
of H(Σ1 |Σ2 ) with respect to Σ2 for two distributions with
the same base radial distribution Q0 . Unfortunately there is
no analytical expression of this gradient except in specific
cases [14]. We shall therefore focus on one of these specific
cases: that of gaussian distributions, which corresponds to Q0
following a χ2 (cK d) distribution, with cK = 1 for K = R,
and cK = 2 for K = C.
The relative entropy of two such gaussian models is given
by [10]:


cK
tr Σ2 −1 Σ1 + log |Σ2 | − log |Σ1 | − d
2
(18)
with Σ1 and Σ2 being the covariance matrices of said
models. In the standard case of gaussian distributions in Rd ,
the resulting likelihood is thus given up to constant additive
terms by:
H(Σ1 |Σ2 ) =

N
α X † −1
xn Σ xn
2N n=1
(19)
 1
1−α
−1
tr Σ Σprior − log |Σ|
−
2
2
The corresponding maximum likelihood estimator is then
given by a weighted mean of the covariance of the prior
distribution and that of the sample covariance matrix of the
new samples (SCM):

This can be differentiated, leading to the following maximum likelihood equation:
Σ=

1
α ((1 − α)Σprior
1 − 2d

N 
α X
1
+
1−
xn xn † )
N n=1
2xn † Σ−1 xn

(22)

A solution can be found numerically by adapting the
algorithm described in [7]:

CG COV

1:
2:
3:

function REG CG COV(Σprior , (xn )1≤n≤N , α, , Kmax )
Σ ← (1 − α)Σprior +
repeat

√

4:

Σ 21 ←

5:

Σ− 12 ← Σ 21 −1

6:

S←

7:
8:
9:
10:

α
N

N
P

xn xn †

n=1

Σ

1
α
(1− 2d
)

((1 − α)Σprior


N 
1
α X
1−
+
xn xn † )
N n=1
2xn † Σ−1 xn


R ← Σ 12 exp Σ− 12 SΣ− 12 − I Σ 21

2 
≤
until tr Σ− 21 ΣΣ− 12 − I
return Σ− 21
end function

Thus in the specific case of gaussian and complex circular
gaussian variables, the regularized maximum likelihood estimators can be computed efficiently.

l(Σ|S) = −

Σ = (1 − α)Σprior +

N
α X
xn xn †
N n=1

(20)

In the complex circular case, one can take into account the
phase symmetry of the sampling distribution, in a manner
similar to Tyler’s estimator; this results in the following
likelihood, up to constant additive terms:

N 

α X
1
† −1
† −1
l(Σ|S) = −
xn Σ xn − log xn Σ xn
N n=1
2


α
− (1 − α)tr Σ−1 Σprior − 1 −
log |Σ|
2d
(21)

IV. PARTIAL ESTIMATION AND REGULARIZED ESTIMATORS
We have thus successfully extended maximum likelihood
estimators of covariance matrices to a regularized estimation
setting. We shall now discuss a further extension of these
regularized estimators in order to improve their robustness to
cases of contamination by outliers. In particular, we would
like to extend the partial maximum likelihood estimation
method developed in [8] to these estimators.
Let us recall that partial estimation as introduced in [8]
consists in the maximization of a partial likelihood:

dP
(x) dS(x)
l|X (P |S) = −H|X (P |S) =
log
dS
x∈X
(23)
under the constrain that X can be any measurable set such
that S(X) ≥ p with p ∈ [0; 1] being the order of the partial
estimation.
Z



Unfortunately this definition is impractical in the case of
regularized estimators for which R = ∨N B ∨ [B]. Indeed one

ARXIV PREPUBLICATION

5

would have to consider partial entropy integrals between the
prior and model distribution of the form:



dPprior
H|X∩[B] (Pprior |P ) =
log
(x) dPprior (x)
dP
x∈X∩[B]
(24)
We therefore propose to restrict the possible partial domains
to subsets X of the form X ∨ [B], with X ⊂ ∨N B and such
that S(X) ≥ (1 − α) + αp. Thus the corresponding restricted
partial likelihood, maximized on the partial domain X is now
given by [8]:
Z

11:

R ← R 21 exp(S − I)R 12

12:

R
R ← tr(R)


2
if tr (R−1 R − I) ≤  then

13:
14:
15:
16:
17:
18:
19:

lp (P |S) = − Hp (P |S)


Z
=

sup
S(X)≥αp+(1−α)

log
x∈X


dP
(25)
(x) dS(x)
dS

This is the restricted partial likelihood of order p, conditional to the fact that the prior distribution Pprior contains no
outliers.
This can be recast as:
dpN e

α X  B
H δxo(n) |P
lp (P |S) = −
dpN e n=1

This is a penalized form of the partial likelihood as expressed in [8]. Thus the partial estimation procedure introduced in [8] can be directly extended to such regularized
estimators.
A. Restricted partial maximum likelihood principle for regularized estimators
We shall now apply the partial estimation procedure to the
different regularized estimators introduced so far.
Outlined below is the algorithmic procedures corresponding
to the partial version of the regularized Tyler estimator for a
distribution in Rd :
function REG P T YLER(Rprior , (xn )1≤n≤N , α, p, , Kmax )
p
2:
R 21 ← Rprior

3:

R− 12 ← R 21 −1

4:

for k from 1 to Kmax do

5:
6:
7:
8:
9:

for n from 1 to N do
τn ← xn † R− 12 R− 12 xn
end for


o ← argsort↑ (τn )1≤n≤N


S ← (1 − α)anscm R− 12 Rprior R− 12
dpN e
X xo(n) xo(n) †
α
1
+
R− 2
R− 21
dpN e
τo(n)
n=1

10:

S←

d
tr(S) S

else
R 12 ←

√

R

R− 21 ← R 12 −1

end if
end for

20:

return R− 21

21:

end function

The corresponding partial procedure for a distribution in
Cd is the same as Tyler but using a different α parameter:

(26)

− (1 − α)H(Pprior |[P ])

1:

break

REG P T YLER (Rprior ,(xn )1≤n≤N ,

α(d−1)
,p,)
d− 21 (1+α)

As stated in III-B, it is in general quite difficult to derive
regularized estimators for M-estimators which are consistent
with maximum likelihood theory as stated in II-Aexcept
for some specific cases. The same remains valid for partial
regularized estimates, for which we shall again concentrate
on the specific case of gaussian models. The corresponding
partial regularized estimators, corresponding to the partial
versions of the regularized sample covariance in the real case
and the REG CG COV algorithm in the complex circular case
are outlined below:
1:

function REG PARTIAL SCM(Σprior , (xn )1≤n≤N , α, Kmax )

2:

Σ ← Σprior

3:

for n from 1 to N do

4:

τn ← xn † Σ−1 xn

5:

end for

6:



o0 ← argsort↑ (xn )1≤n≤N

7:

for k from 1 to Kmax do
α
dpN e

8:

Σ ← (1 − α)Σprior +

9:

for n from 1 to N do

10:

dpN
Pe
n=1

τn ← xn † Σ−1 xn

11:

end for

12:



o1 ← argsort↑ (τn )1≤n≤N

xo0 (n) xo0 (n) †

ARXIV PREPUBLICATION

13:
14:
15:
16:
17:

6

if o0 = o1 then

the background noise is generated as a white gaussian noise
of unit variance.
The target signal is generated as a complex centered circular
1-dimensional gaussian signal aligned with the test signal s,
whose variance σ is such that:

break
else
o0 ← o1
end if

18:

end for

19:

return Σ

20:

end function

10 log10 (σ) = SiNR

In the complex circular case the corresponding algorithm
is given by:
1:
2:
3:

function REG PCG COV(Σprior (xn )1≤n≤N , α, , Kmax )
Σ ← (1 − α)Σprior +
√
Σ 12 ← Σ

α
N

N
P

Σ− 12 ← Σ 12 −1

5:

for k from 1 to Kmax do

7:

for n from 1 to N do
end for

9:

o ← argsort↑

11:
12:
13:
14:
15:
16:
17:

S←

1
α
(1− 2d
)

τn −

1
2



log (τn ) 1≤n≤N

This covariance/correlation matrix is then used in an
adaptive detection scheme.


(1 − α)Σprior



dpN e 
α X
1
†
+
xo(n) xo(n)
1−
dpN e n=1
2τo(n)


R ← Σ 12 exp Σ− 12 SΣ− 12 Σ 12

2 
if tr Σ− 21 ΣΣ− 12 − I
≤  then

The performances are considered in the asymptotic regime
(meaning that Σm follows the same distribution as Σm−1 for
any m ∈ Z). The performances are shown for N = 11, p = 34
for partial estimators, with various values of α; the results are
also compared to the non regularized versions of the estimators
for N = 22 [7][8].

break
else
√
Σ 21 ← Σ

Σ− 12 ← Σ 21 −1

end if

18:

end for

19:

return Σ− 21

20:



In this scenario, we have an infinite stream of learning
samples (xn,m )1≤n≤N,m∈Z . These learning samples are used
to learn the covariance/correlation matrix of the signal sequentially, using the following pattern:

Σm = reg(Σm−1 , (xn,m )1≤n≤N )

τn ← xn † Σ−1 xn

8:

10:

xn xn †

n=1

4:

6:

The detection thresholds are defined to have a false alarm
rate of 10−4 ; they are learned on clean training sets, that is
that they contain no outliers.

end function

V. S IMULATIONS
We shall now show some simulation results of adaptive
detectors using various estimators introduced in this article,
using the detection tests introduced in [7][24].
The simulation results are shown in a single channel scenario of dimension d = 8, as a function of the SiNR of the
target signal.

Fig. 1. Detection capabilities of Tyler-NMF and pTyler-NMF

ARXIV PREPUBLICATION

7

presence of outliers, by using the partial procedure defined in
[8].
Finally the potential gain in terms of detection performances
in adaptive detection schemes is shown on some simulations.
A PPENDIX A
O N THE RELATIVE ENTROPY OF ANGULAR MODELS
The relative entropy of two angular models of Rd , characterized by their correlation matrix R0 and R (with the convention
that |R0 | = |R| = 1), is given by:
H(R0 |R) =
 † −1 
Z
x R x
d−1
log
dsd−1
2
x† R0 −1 x
x∈Sd−1
Fig. 2. Detection capabilities of reg SCM-MF and reg pSCM-MF

L(R )−1 x
p 0
x† R0 −1 x

!

(27)
Although the exact computation of this integral represents
a challenge, what is of interest to us is its differential with
respect to the matrix R:
dH(R0 |R) =

Z
d−1 
tr 
2

xx

†

−1

x† R−1 x

dsd−1

x∈Sd−1

L(R ) x
p 0
x† R0 −1 x



!


d(R−1 )
(28)

Thus we have to compute:
xx†
dsd−1
x† R−1 x

Z
Σ=
x∈Sd−1

L(R )−1 x
p 0
x† R0 −1 x

!
(29)

1

By making a change of variable y = R− 2 x, this integral
can be recast as:




†
− 12 −1
− 21
yy
R
R
)
y
L(R
0
 R 12
Σ=R 
dsd−1  q
†
1
1
y∈Sd−1 y y
y † R− 2 R0 R− 2 y
(30)
Let us note for any symmetric positive matrix Σ:
1
2

As can be seen on figures 1, 2 and V, the use of a regularized
estimate in this case allows to greatly improve the detection
capabilities of adaptive tests. Indeed the performances quickly
approach the asymptotic limit in which the covariance matrix
is perfectly estimated when a slower integration factor is
chosen.
Moreover the degradation of performances of partial estimators can be almost completely negated for a slow enough
integration (α ≤ 14 ), thus resulting in estimators which offer
performances close to non-adaptive detectors, while retaining
the robustness to outliers offered by partial estimators.
VI. C ONCLUSION
We have presented a theoretical background for a likelihood
based approach to regularized estimations. This background
allows to define penalties used in penalized maximum likelihood schemes in a way which is consistent with the maximum
likelihood principle.
This maximum likelihood principle for regularized estimation
is then further extended to allow for robust estimators in the

Z

Z
anscm(Σ) =
y∈Sd−1

yy †
dsd−1
y† y

L(Σ)−1 y
p
y † Σ−1 y

!
(31)

†

This is the integral of yy
on the unit ellipse defined by the
y† y
correlation matrix Σ, and therefore corresponds to the expected
value of the normalized sample covariance estimator (NSCM)
for an elliptical distribution of correlation matrix proportional
to Σ. The structure of this matrix has been studied thoroughly
in [3] for complex circular distributions, in the case in which
all eigenvalues are distinct:
anscm(Σ) = U diag(µ)U †

(32)

with:
•

Σ = U diag(λ)U † with U U † = I being the eigen
decomposition of Σ

ARXIV PREPUBLICATION

•

µ = Ec with:









E
=








c=



8

λ
log( λj )
i
λj
λi −1

−



λi
λj

if i 6= j 
if i = j
1≤i,j≤d

0
!

1
λk
k6=j 1− λj

Q

1≤j≤d

Note that in a practical implementation, the limitation of
having distinct eigenvalues is not restrictive, as it is almost
surely the case. However one should take steps to insure a
decent numerical stability. We have on our end opted for
enforcing a minimum spacing between eigenvalues determined
by a numerical tolerance, as given by the following algorithm:
1:
2:

function R ESPACE
 ((λk )1≤k≤d
 , )
o ← argsort↑ (λk )1≤k≤d

3:

(µk )1≤k≤d ← (λk )1≤k≤d

4:

for k from 1 to d − 1 do
if

5:

µo(k+1)
µo(k)

− 1 <  then

for m from k + 1 to d do

6:

µo(m) ← (1 + )µo(m−1) − µo(m)

7:

end for

8:

10:

end if
end for

11:

return (µk )1≤k≤d

9:

12:

end function

Unfortunately we do not have a similar result in the general
case of distributions in Rd , in which case one has to resort
to numerical integration methods. Indeed the matrix can be
expressed as anscm U diag(λ)U † = U diag(µ)U † , with:
µk =
 ZZ
2 Γ d2

√
π Γ d−1
2

t∈[0;1]
x∈Sd−2

d

t2 (1 − t2 ) 2
P
t2 + (1 − t2 )

−1

m6=k

λm
2
λk xk

dsd−2 (x)dt

(33)
There is no known analytical expression of these eigenvalues to the authors’ knowledge.
R EFERENCES
[1] John Aldrich et al. Ra fisher and the making of maximum likelihood
1912-1922. Statistical Science, 12(3):162–176, 1997.
[2] Frédéric Barbaresco. Super-resolution spectrum analysis regularization:
Burg, capon & ago-antagonistic algorithms. In European Signal Processing Conference, 1996. EUSIPCO 1996. 8th, pages 1–4. IEEE, 1996.
[3] Sebastien Bausson, Frederic Pascal, Philippe Forster, Jean-Philippe
Ovarlez, and Pascal Larzabal. First-and second-order moments of the
normalized sample covariance matrix of spherically invariant random
vectors. IEEE Signal Processing Letters, 14(6):425–428, 2007.
[4] Peter J Bickel and Elizaveta Levina. Covariance regularization by
thresholding. The Annals of Statistics, pages 2577–2604, 2008.

[5] MA Chmielewski. Elliptically symmetric distributions: a review and
bibliography. International Statistical Review/Revue Internationale de
Statistique, pages 67–74, 1981.
[6] Philippe Ciuciu, Jérôme Idier, and Jean-François Giovannelli. Regularized estimation of mixed spectra using a circular Gibbs-Markov model.
IEEE Transactions on Signal Processing, 49(10):2201–2213, oct. 2001.
[7] Christophe Culan and Claude Adnet. Maximum likelihood estimation of covariances of elliptically symmetric distributions. eprint,
arXiv:1611.04365, 2016.
[8] Christophe Culan and Claude Adnet. Partial estimators and application
to covariance estimation of gaussian and elliptical distributions. eprint,
arXiv:1611.07418, 2016.
[9] Junbo Duan, Charles Soussen, David Brie, Jérôme Idier, Mingxi Wan,
and Yu-Ping Wang. Generalized LASSO with under-determined regularization matrices. Signal Processing, 127:239–246, 2016.
[10] John Duchi. Derivations for linear algebra and optimization. Berkeley,
California, 2007.
[11] Jean-François Giovannelli, Jérôme Idier, Guy Desodt, and Daniel Muller.
Regularized adaptive long autoregressive spectral analysis. IEEE Transactions on Geoscience and Remote Sensing, 39(10):2194–2202, oct.
2001.
[12] Anders Hald. On the history of maximum likelihood in relation to
inverse probability and least squares. Statistical Science, pages 214–
222, 1999.
[13] Joseph P Hoffbeck and David A Landgrebe. Covariance matrix estimation and classification with limited training data. IEEE Transactions on
Pattern Analysis and Machine Intelligence, 18(7):763–767, 1996.
[14] Reshad Hosseini, Suvrit Sra, Lucas Theis, and Matthias Bethge. Inference and mixture modeling with the elliptical gamma distribution.
Computational Statistics & Data Analysis, 101:29–43, 2016.
[15] Solomon Kullback and Richard A Leibler. On information and sufficiency. The annals of mathematical statistics, 22(1):79–86, 1951.
[16] Ali Mohammad-Djafari, Jean-François Giovannelli, Guy Demoment,
and Jérôme Idier. Regularization, maximum entropy and probabilistic
methods in mass spectrometry data processing problems. Int. Journal
of Mass Spectrometry, 215(1-3):175–193, avr. 2002.
[17] Esa Ollila and David E Tyler. Regularized-estimators of scatter matrix.
IEEE Transactions on Signal Processing, 62(22):6059–6070, 2014.
[18] Esa Ollila, David E Tyler, Visa Koivunen, and H Vincent Poor. Complex
elliptically symmetric distributions: Survey, new results and applications.
Signal Processing, IEEE Transactions on, 60(11):5597–5625, 2012.
[19] Frédéric Pascal, Yacine Chitour, Jean-Philippe Ovarlez, Philippe Forster,
and Pascal Larzabal. Covariance structure maximum-likelihood estimates in compound gaussian noise: Existence and algorithm analysis.
Signal Processing, IEEE Transactions on, 56(1):34–48, 2008.
[20] Johann Pfanzagl. Parametric statistical theory. Walter de Gruyter, 1994.
[21] John W Pratt. Fy edgeworth and ra fisher on the efficiency of maximum
likelihood estimation. The Annals of Statistics, pages 501–514, 1976.
[22] Pradeep Ravikumar, Martin J Wainwright, Garvesh Raskutti, Bin Yu,
et al. High-dimensional covariance estimation by minimizing 1penalized log-determinant divergence. Electronic Journal of Statistics,
5:935–980, 2011.
[23] Louis L Scharf. Statistical signal processing, volume 98. AddisonWesley Reading, MA, 1991.
[24] Louis L Scharf and Benjamin Friedlander. Matched subspace detectors.
Signal Processing, IEEE Transactions on, 42(8):2146–2157, 1994.
[25] Stephen M Stigler. Statistics on the table: The history of statistical
concepts and methods. Harvard University Press, 2002.
[26] Ying Sun, Prabhu Babu, and Daniel P Palomar. Regularized tyler’s scatter estimator: Existence, uniqueness, and algorithms. IEEE Transactions
on Signal Processing, 62(19):5143–5156, 2014.
[27] David E Tyler. Statistical analysis for the angular central gaussian
distribution on the sphere. Biometrika, 74(3):579–589, 1987.
[28] Ami Wiesel. Geodesic convexity and covariance estimation. Signal
Processing, IEEE Transactions on, 60(12):6182–6189, 2012.
[29] Ami Wiesel. Unified framework to regularized covariance estimation
in scaled gaussian models. IEEE Transactions on Signal Processing,
60(1):29–38, 2012.
[30] Daniela M Witten and Robert Tibshirani. Covariance-regularized regression and classification for high dimensional problems. Journal of the
Royal Statistical Society: Series B (Statistical Methodology), 71(3):615–
636, 2009.
[31] Lingzhou Xue, Shiqian Ma, and Hui Zou. Positive-definite 1-penalized
estimation of large covariance matrices. Journal of the American
Statistical Association, 107(500):1480–1491, 2012.

ARXIV PREPUBLICATION

Christophe Culan was born in Villeneuve-St-Georges, France, on November
1st , 1988. He received jointly the engineering degree from Ecole Centrale de
Paris (ECP), France and a master’s degree in Physico-informatics from Keio
University, Japan, in 2013.
He was a researcher in applied physics in Itoh laboratory from 2013 to
2014, specialized in quantum information and quantum computing, and has
contributed to several publications related to these subjects.
He currently holds a position as a research engineer in Thales Air Systems,
Limours, France, in the Advanced Radar Concepts division. His current research interests include statistical signal and data processing, robust statistics,
machine learning and information geometry.

Claude Adnet was born in A, France, in 1961. He received the DEA degree
in signal processing and Phd degree in 1988 and 1991 respectively, from the
Institut National Polytechnique de Grenoble (INPG), Grenoble France. Since
then, he has been working for THALES Group, where he is now Senior
Scientist. His research interests include radar signal processing and radar data
processing.

9

