Riemannian Laplace Distribution on the Space of Symmetric Positive Definite Matrices

Abstract

The Riemannian geometry of the space ${\mathcal{P}}_m$, of $m\times m$ symmetric positive definite matrices, has provided effective tools to the fields of medical imaging, computer vision and radar signal processing. Still, an open challenge remains, which consists of extending these tools to correctly handle the presence of outliers (or abnormal data), arising from excessive noise or faulty measurements. The present paper tackles this challenge by introducing new probability distributions, called Riemannian Laplace distributions on the space ${\mathcal{P}}_m$. First, it shows that these distributions provide a statistical foundation for the concept of the Riemannian median, which offers improved robustness in dealing with outliers (in comparison to the more popular concept of the Riemannian center of mass). Second, it describes an original expectation-maximization algorithm, for estimating mixtures of Riemannian Laplace distributions. This algorithm is applied to the problem of texture classification, in computer vision, which is considered in the presence of outliers. It is shown to give significantly better performance with respect to other recently-proposed approaches.

1. Introduction

Data with values in the space ${\mathcal{P}}_m$, of $m\times m$ symmetric positive definite matrices, play an essential role in many applications, including medical imaging [1,2], computer vision [3,4,5,6,7] and radar signal processing [8,9]. In these applications, the location where a dataset is centered has a special interest. While several definitions of this location are possible, its meaning as a representative of the set should be clear. Perhaps, the most known and well-used quantity to represent a center of a dataset is the Fréchet mean. Given a set of points $Y_1,\cdots ,Y_n$ in ${\mathcal{P}}_m$, their Fréchet mean is defined to be:

$$\text{Mean}(Y_1,\cdots ,Y_n)={\text{argmin}}_{Y\in {\mathcal{P}}_m}\sum _{i=1}^n{d}^2(Y,Y_i)$$

(1)

where d is Rao’s Riemannian distance on ${\mathcal{P}}_m$ [10,11].

Statistics on general Riemannian manifolds have been powered by the development of different tools for geometric measurements and new probability distributions on manifolds [12,13]. On the manifold $({\mathcal{P}}_m,d)$, the major advances in this field have been achieved by the recent papers [14,15], which introduce the Riemannian Gaussian distribution on $({\mathcal{P}}_m,d)$. This distribution depends on two parameters $\overline{Y}\in {\mathcal{P}}_m$ and $\sigma >0$, and its density with respect to the Riemannian volume form $dv(Y)$ of ${\mathcal{P}}_m$ (see Formula (13) in Section 2) is:

$$\frac{1}{Z_m(\sigma )}exp\left[-\frac{d^2(Y,\overline{Y})}{2{\sigma }^2}\right]$$

(2)

where $Z_m(\sigma )$ is a normalizing factor depending only on σ (and not on $\overline{Y}$).

For the Gaussian distribution Equation (2), the maximum likelihood estimate (MLE) for the parameter $\overline{Y}$ based on observations $Y_1,\cdots ,Y_n$ corresponds to the mean Equation (1). In [15], a detailed study of statistical inference for this distribution was given and then applied to the classification of data in ${\mathcal{P}}_m$, showing that it yields better performance, in comparison to recent approaches [2].

When a dataset contains extreme values (or outliers), because of the impact of these values on $d^2$, the mean becomes less useful. It is usually replaced with the Riemannian median:

$$\text{Median}(Y_1,\cdots ,Y_n)={\text{argmin}}_{Y\in {\mathcal{P}}_m}\sum _{i=1}^{n}d(Y,Y_i)$$

(3)

Definition Equation (3) corresponds to that of the median in statistics based on ordering of the values of a sequence. However, this interpretation does not continue to hold on ${\mathcal{P}}_m$. In fact, the Riemannian distance on ${\mathcal{P}}_m$ is not associated with any norm, and it is therefore only possible to compare distances of a set of matrices to a reference matrix.

In the presence of outliers, the Gaussian distribution on ${\mathcal{P}}_m$ also loses its robustness properties. The main contribution of the present paper is to remedy this problem by introducing the Riemannian Laplace distribution while maintaining the same one-to-one relation between MLE and the Riemannian median. This will be shown to offer considerable improvement in dealing with outliers.

This paper is organized as follows.

Section 2 reviews the Riemannian geometry of ${\mathcal{P}}_m$, when this manifold is equipped with the Riemannian metric known as the Rao–Fisher or affine invariant metric [10,11]. In particular, it gives analytic expressions for geodesic curves, Riemannian distance and recalls the invariance of Rao’s distance under affine transformations.

Section 3 introduces the Laplace distribution $\mathcal{L}(\overline{Y},\sigma )$ through its probability density function with respect to the volume form $dv(Y)$:

$$p(Y|\overline{Y},\sigma )=\frac{1}{{\zeta }_m(\sigma )}exp\left[-\frac{d(Y,\overline{Y})}{\sigma }\right]$$

here, σ lies in an interval $]0,{\sigma }_{\text{max}}[$ with ${\sigma }_{\text{max}}<\infty$. This is because the normalizing constant ${\zeta }_m(\sigma )$ becomes infinite for $\sigma \ge {\sigma }_{\text{max}}$. It will be shown that ${\zeta }_m(\sigma )$ depend only on σ (and not on $\overline{Y}$) for all $\sigma <{\sigma }_{\text{max}}$. This important fact leads to simple expressions of MLEs of $\overline{Y}$ and σ. In particular, the MLE of $\overline{Y}$ based on a family of observations $Y_1,\cdots ,Y_N$ sampled from $\mathcal{L}(\overline{Y},\sigma )$ is given by the median of $Y_1,\cdots ,Y_N$ defined by Equation (3) where d is Rao’s distance.

Section 4 focuses on mixtures of Riemannian Laplace distributions on ${\mathcal{P}}_m$. A distribution of this kind has a density:

$$p(Y|{({\omega }_{\mu },{\overline{Y}}_{\mu },{\sigma }_{\mu })}_{1\le \mu \le M})=\sum _{\mu =1}^M{\varpi }_{\mu }p(Y|{\overline{Y}}_{\mu },{\sigma }_{\mu })$$

(4)

with respect to the volume form $dv(Y)$. Here, M is the number of mixture components, ${\varpi }_{\mu }>0,{\overline{Y}}_{\mu }\in {\mathcal{P}}_m,{\sigma }_{\mu }>0$ for all $1\le \mu \le M$ and ${\sum }_{\mu =1}^M{\varpi }_{\mu }=1$. A new EM (expectation-maximization) algorithm that computes maximum likelihood estimates of the mixture parameters ${({\varpi }_{\mu },{\overline{Y}}_{\mu },{\sigma }_{\mu })}_{1\le \mu \le M}$ is provided. The problem of the order selection of the number M in Equation (4) is also discussed and performed using the Bayesian information criterion (BIC) [16].

Section 5 is an application of the previous material to the classification of data with values in ${\mathcal{P}}_m$, which contain outliers (abnormal data points). Assume to be given a training sequence $Y_1,\cdots ,Y_n\in {\mathcal{P}}_m$. Using the EM algorithm developed in Section 4, it is possible to subdivide this sequence into disjoint classes. To classify new data points, a classification rule is proposed. The robustness of this rule lies in the fact that it is based on the distances between new observations and the respective medians of classes instead of the means [15]. This rule will be illustrated by an application to the problem of texture classification in computer vision. The obtained results show improved performance with respect to recent approaches which use the Riemannian Gaussian distribution [15] and the Wishart distribution [17].

2. Riemannian Geometry of ${\mathcal{P}}_m$

The geometry of Siegel homogeneous bounded domains, such as Kähler homogeneous manifolds, have been studied by Felix A. Berezin [18] and P. Malliavin [19]. The structure of Kähler homogeneous manifolds has been used in [20,21] to parameterize (Toeplitz–) Block–Toeplitz matrices. This led to a Hessian metric from information geometry theory with a Kähler potential given by entropy and to an algorithm to compute medians of (Toeplitz–)Block–Toeplitz matrices by Karcher flow on Mostow/Berger fibration of a Siegel disk. Optimal numerical schemes of this algorithm in a Siegel disk have been studied, developed and validated in [22,23,24].

This section introduces the necessary background on the Riemannian geometry of ${\mathcal{P}}_m$, the space of symmetric positive definite matrices of size $m\times m$. Precisely, ${\mathcal{P}}_m$ is equipped with the Riemannian metric known as the affine-invariant metric. First, analytic expressions are recalled for geodesic curves and Riemannian distance. Then, two properties are stated, which are fundamental to the following. These are affine-invariance of the Riemannian distance and the existence and uniqueness of Riemannian medians.

The affine-invariant metric, called the Rao–Fisher metric in information geometry, has the following expression:

$$g_Y(A,B)=\text{tr}(Y^{-1}A{Y}^{-1}B)$$

(5)

where $Y\in {\mathcal{P}}_m$ and $A,B\in T_Y{\mathcal{P}}_m$, the tangent space to ${\mathcal{P}}_m$ at Y, which is identified with the vector space of $m\times m$ symmetric matrices. The Riemannian metric Equation (5) induces a Riemannian distance on ${\mathcal{P}}_m$ as follows. The length of a smooth curve $c:[0,1]\to {\mathcal{P}}_m$ is given by:

$$L(c)={\int }_0^1\sqrt{g_{c(t)}(\dot{c}(t),\dot{c}(t))}dt$$

(6)

where $\dot{c}(t)=\frac{dc}{dt}$. For $Y,Z\in {\mathcal{P}}_m$, the Riemannian distance $d(Y,Z)$, called Rao’s distance in information geometry, is defined to be:

$$d(Y,Z)=inf\left\{L(c),c:[0,1]\to {\mathcal{P}}_m\text{is}\text{a}\text{smooth}\text{curve}\text{with}c(0)=Y,c(1)=Z\right\}.$$

This infimum is achieved by a unique curve $c=\gamma$, called the geodesic connecting Y and Z, which has the following equation [10,25]:

$$\gamma (t)=Y^{1/2}{(Y^{-1/2}Z{Y}^{-1/2})}^t{Y}^{1/2}$$

(7)

Here, and throughout the following, all matrix functions (for example, square root, logarithm or power) are understood as symmetric matrix functions [26]. By definition, $d(Y,Z)$ coincides with $L(\gamma )$, which turns out to be:

$$d^2(Y,Z)=\mathrm{tr}{[log(Y^{-1/2}Z{Y}^{-1/2})]}^2$$

(8)

Equipped with the affine-invariant metric Equation (5), the space ${\mathcal{P}}_m$ enjoys two useful properties, which are the following. The first property is invariance under affine transformations [10,25]. Recall that an affine transformation of ${\mathcal{P}}_m$ is a mapping $Y\mapsto Y·A$, where A is an invertible real matrix of size $m\times m$,

$$Y·A=A^{†}YA$$

(9)

and ${}^{†}$ denotes the transpose. Denote by $\mathrm{GL}(m)$ the group of $m\times m$ invertible real matrices on ${\mathcal{P}}_m$. Then, the action of $\mathrm{GL}(m)$ on ${\mathcal{P}}_m$ is transitive. This means that for any $Y,Z\in {\mathcal{P}}_m$, there exists $A\in \mathrm{GL}(m)$, such that $Y.A=Z$. Moreover, the Riemannian distance Equation (8) is invariant by affine transformations in the sense that for all $Y,Z\in {\mathcal{P}}_m$:

$$d(Y,Z)=d(Y·A,Z·A)$$

(10)

where $Y·A$ and $Z·A$ are defined by Equation (9). The transitivity of the action Equation (9) and the isometry property Equation (10) make ${\mathcal{P}}_m$ a Riemannian homogeneous space.

The affine-invariant metric Equation (5) turns ${\mathcal{P}}_m$ into a Riemannian manifold of negative sectional curvature [10,27]. As a result, ${\mathcal{P}}_m$ enjoys the property of the existence and uniqueness of Riemannian medians. The Riemannian median of N points $Y_1,\cdots ,Y_N\in {\mathcal{P}}_m$ is defined to be:

$${\widehat{Y}}_N={\text{argmin}}_Y\sum _{n=1}^{N}d(Y,Y_n)$$

(11)

where $d(Y,Y_n)$ is the Riemannian distance Equation (8). If $Y_1,\cdots ,Y_N$ do not belong to the same geodesic, then ${\widehat{Y}}_N$ exists and is unique [28]. More generally, for any probability measure π on ${\mathcal{P}}_m$, the median of π is defined to be:

$${\widehat{Y}}_{\pi }={\text{argmin}}_Y{\int }_{{\mathcal{P}}_m}d(Y,Z)d\pi (Z)$$

(12)

Note that Equation (12) reduces to Equation (11) for $\pi =\frac{1}{N}{\sum }_{n=1}^N{\delta }_{Y_n}$. If the support of π is not carried by a single geodesic, then again, ${\widehat{Y}}_{\pi }$ exists and is unique by the main result of [28].

To end this section, consider the Riemannian volume associated with the affine-invariant Riemannian metric [10]:

$$dv(Y)=det{(Y)}^{-\frac{m+1}{2}}\prod _{i\le j}d{Y}_{ij}$$

(13)

where the indices denote matrix elements. The Riemannian volume is used to define the integral of a function $f:{\mathcal{P}}_m\to \mathbb{R}$ as:

$${\int }_{{\mathcal{P}}_m}f(Y)dv(Y)=\int \dots \int f(Y)det{(Y)}^{-\frac{m+1}{2}}\prod _{i\le j}d{Y}_{ij}$$

(14)

where the integral on the right-hand side is a multiple integral over the $m(m+1)/2$ variables, $Y_{ij}$ with $i\le j$. The integral Equation (14) is invariant under affine transformations. Precisely:

$${\int }_{{\mathcal{P}}_m}f(Y·A)dv(Y)={\int }_{{\mathcal{P}}_m}f(Y)dv(Y)$$

(15)

where $Y·A$ is the affine transformation given by Equation (9). It takes on a simplified form when $f(Y)$ only depends on the eigenvalues of Y. Precisely, let the spectral decomposition of Y be given by $Y=U^{†}\mathrm{diag}(e^{r_1},\cdots ,e^{r_m})U$, where U is an orthogonal matrix and $e^{r_1},\cdots ,e^{r_m}$ are the eigenvalues of Y. Assume that $f(Y)=f(r_1,\dots ,r_m)$, then the invariant integral Equation (14) reduces to:

$${\int }_{{\mathcal{P}}_m}f(Y)dv(Y)=c_m\times {\int }_{{\mathbb{R}}^m}f(r_1,\cdots ,r_m)\prod _{i<j}sinh\left(\frac{|r_i-r_j|}{2}\right)d{r}_1\cdots d{r}_m$$

(16)

where the constant $c_m$ is given by $c_m=\frac{1}{m!}\times {\omega }_m\times 8^{\frac{m(m-1)}{4}},{\omega }_m=\frac{{\pi }^{m^2/2}}{{\Gamma }_m(m/2)}$ and ${\Gamma }_m$ is the multivariate gamma function given in [29]. See Appendix A for the derivation of Equation (16) from Equation (14).

3. Riemannian Laplace Distribution on ${\mathcal{P}}_m$

3.1. Definition of $\mathcal{L}(\overline{Y},\sigma )$

The Riemannian Laplace distribution on ${\mathcal{P}}_m$ is defined by analogy with the well-known Laplace distribution on $\mathbb{R}$. Recall the density of the Laplace distribution on $\mathbb{R}$,

$$p(x|\overline{x},\sigma )=\frac{1}{2\sigma }{e}^{-|x-\overline{x}|/\sigma }$$

where $\overline{x}\in \mathbb{R}$ and $\sigma >0$. This is a density with respect to the length element $dx$ on $\mathbb{R}$. The density of the Riemannian Laplace distribution on ${\mathcal{P}}_m$ will be given by:

$$p(Y|\overline{Y},\sigma )=\frac{1}{{\zeta }_m(\sigma )}exp\left[-\frac{d(Y,\overline{Y})}{\sigma }\right]$$

(17)

here, $\overline{Y}\in {\mathcal{P}}_m$, $\sigma >0$, and the density is with respect to the Riemannian volume element Equation (13) on ${\mathcal{P}}_m$. The normalizing factor ${\zeta }_m(\sigma )$ appearing in Equation (17) is given by the integral:

$${\int }_{{\mathcal{P}}_m}exp\left[-\frac{d(Y,\overline{Y})}{\sigma }\right]dv(Y)$$

Assume for now that this integral is finite for some choice of $\overline{Y}$ and σ. It is possible to show that its value does not depend on $\overline{Y}$. To do so, recall that the action of $\mathrm{GL}(m)$ on ${\mathcal{P}}_m$ is transitive. As a consequence, there exists $A\in {\mathcal{P}}_m$, such that $\overline{Y}=I.A$, where $I.A$ is defined as in Equation (9). From Equation (10), it follows that $d(Y,\overline{Y})=d(Y,I.A)=d(Y.A^{-1},I)$. From the invariance property Equation (15):

$${\int }_{{\mathcal{P}}_m}exp\left[-\frac{d(Y,\overline{Y})}{\sigma }\right]dv(Y)={\int }_{{\mathcal{P}}_m}exp\left[-\frac{d(Y,I)}{\sigma }\right]dv(Y)$$

(18)

The integral on the right does not depend on $\overline{Y}$, which proves the above claim.

The last integral representation and formula Equation (16) lead to the following explicit expression:

$${\zeta }_m(\sigma )=c_m\times {\int }_{{\mathbb{R}}^m}{e}^{-\frac{|r|}{\sigma }}\prod _{i<j}sinh\left(\frac{|r_i-r_j|}{2}\right)d{r}_1\cdots d{r}_m$$

(19)

where $|r|={(r_1^2+\cdots +r_2^m)}^{\frac{1}{2}}$ and $c_m$ is the same constant as in Equation (16) (see Appendix B for more details on the derivation of Equation (19)).

A distinctive feature of the Riemannian Laplace distribution on ${\mathcal{P}}_m$, in comparison to the Laplace distribution on $\mathbb{R}$ is that there exist certain values of σ for which it cannot be defined. This is because the integral Equation (19) diverges for certain values of this parameter. This leads to the following definition.

Definition 1. 

Set ${\sigma }_m=sup\{\sigma >0:{\zeta }_m(\sigma )<\infty \}$. Then, for $\overline{Y}\in {\mathcal{P}}_m$ and $\sigma \in (0,{\sigma }_m)$, the Riemannian Laplace distribution on ${\mathcal{P}}_m$, denoted by $\mathcal{L}(\overline{Y},\sigma )$, is defined as the probability distribution on ${\mathcal{P}}_m$, whose density with respect to $dv(Y)$ is given by Equation (17), where ${\zeta }_m(\sigma )$ is defined by Equation (19).

The constant ${\sigma }_m$ in this definition satisfies $0<{\sigma }_m<\infty$ for all m and takes the value $\sqrt{2}$ for $m=2$ (see Appendix C for proofs).

3.2. Sampling from $\mathcal{L}(\overline{Y},\sigma )$

The current section presents a general method for sampling from the Laplace distribution $\mathcal{L}(\overline{Y},\sigma )$. This method relies in part on the following transformation property.

Proposition 2. 

Let Y be a random variable in ${\mathcal{P}}_m$. For all $A\in \mathrm{GL}(m)$,

$$Y\sim \mathcal{L}(\overline{Y},\sigma )⟹Y·A\sim \mathcal{L}(\overline{Y}·A,\sigma )$$

where $Y·A$ is given by Equation (9).

Proof. 

Let $\phi :{\mathcal{P}}_m\to \mathbb{R}$ be a test function. If $Y\sim \mathcal{L}(\overline{Y},\sigma )$ and $Z=Y·A$, then the expectation of $\phi (Z)$ is given by:

$$\begin{array}{c}\hfill E[\phi (Z)]={\int }_{{\mathcal{P}}_m}\phi (X·A)p(X|\overline{Y},\sigma )dv(X)={\int }_{{\mathcal{P}}_m}\phi (X)p(X·A^{-1}|\overline{Y},\sigma )dv(X)\end{array}$$

where the equality is a result of Equation (15). However, $p(X·A^{-1}|\overline{Y},\sigma )=p(X|\overline{Y}·A,\sigma )$ by Equation (10), which proves the proposition.

 ☐

The following algorithm describes how to sample from $\mathcal{L}(\overline{Y},\sigma )$ where $0<\sigma <{\sigma }_m$. For this, it is first required to sample from the density p on ${\mathbb{R}}^m$ defined by:

$$p(r)=\frac{c_m}{{\zeta }_m(\sigma )}{e}^{-\frac{|r|}{\sigma }}\prod _{i<j}sinh\left(\frac{|r_i-r_j|}{2}\right),r=(r_1,\cdots ,r_m).$$

This can be done by a usual Metropolis algorithm [30].

It is also required to sample from the uniform distribution on $\mathrm{O}(m)$, the group of real orthogonal $m\times m$ matrices. This can be done by generating A, an $m\times m$ matrix, whose entries are i.i.d. with normal distribution $\mathcal{N}(0,1)$, then the orthogonal matrix U, in the decomposition $A=UT$ with T upper triangular, is uniformly distributed on $\mathrm{O}(m)$ [29] (p. 70). Sampling from $\mathcal{L}(\overline{Y},\sigma )$ can now be described as follows.

Algorithm 1 Sampling from $\mathcal{L}(\overline{Y},\sigma )$.

1: Generate i.i.d. samples $(r_1,\cdots ,r_m)\in {\mathbb{R}}^m$ with density p    2: Generate U from a uniform distribution on $\mathrm{O}(m)$    3: $X\leftarrow U^{†}\text{diag}(e^{r_1},\cdots ,e^{r_m})U$    4: $Y\leftarrow X.{\overline{Y}}^{\frac{1}{2}}$

Note that the law of X in Step 3 is $\mathcal{L}(I,\sigma )$; the proof of this fact is given in Appendix D. Finally, the law of Y in Step 4 is $\mathcal{L}(I.{\overline{Y}}^{\frac{1}{2}}=\overline{Y},\sigma )$ by proposition Equation (2).

3.3. Estimation of $\overline{Y}$ and σ

The current section considers maximum likelihood estimation of the parameters $\overline{Y}$ and σ, based on independent observations $Y_1,\dots ,Y_N$ from the Riemannian Laplace distribution $\mathcal{L}(\overline{Y},\sigma )$. The main results are contained in Propositions 3 and 5 below.

Proposition 3 states the existence and uniqueness of the maximum likelihood estimates ${\widehat{Y}}_N$ and ${\widehat{\sigma }}_N$ of $\overline{Y}$ and σ. In particular, the maximum likelihood estimate ${\widehat{Y}}_N$ of $\overline{Y}$ is the Riemannian median of $Y_1,\dots ,Y_N$, defined by Equation (11). Numerical computation of ${\widehat{Y}}_N$ will be considered and carried out using a Riemannian sub-gradient descent algorithm [8].

Proposition 5 states the convergence of the maximum likelihood estimate ${\widehat{Y}}_N$ to the true value of the parameter $\overline{Y}$. It is based on Lemma 4, which states that the parameter $\overline{Y}$ is the Riemannian median of the distribution $\mathcal{L}(\overline{Y},\sigma )$ in the sense of definition Equation (12).

Proposition 3 

(MLE and median). The maximum likelihood estimate of the parameter $\overline{Y}$ is the Riemannian median ${\widehat{Y}}_N$ of $Y_1,\dots ,Y_N$. Moreover, the maximum likelihood estimate of the parameter σ is the solution ${\widehat{\sigma }}_N$ of:

$${\sigma }^2\times \frac{d}{\phantom{(}d\sigma }log{\zeta }_m(\sigma )=\frac{1}{N}\sum _{n=1}^{N}d(\overline{Y},Y_n)$$

(20)

Both ${\widehat{Y}}_N$ and ${\widehat{\sigma }}_N$ exist and are unique for any realization of the samples $Y_1,\dots ,Y_N$.

 

Proof of Proposition 3. The log-likelihood function, of the parameters $\overline{Y}$ and σ, can be written as:

$$\begin{array}{ccc}\hfill \sum _{n=1}^{N}logp(Y_n|\overline{Y},\sigma )& =& \sum _{n=1}^{N}log\left(\frac{1}{{\zeta }_m(\sigma )}{e}^{-\frac{d(\overline{Y},Y_n)}{\sigma }}\right)\hfill \\ & =& -Nlog{\zeta }_m(\sigma )-\frac{1}{\sigma }\sum _{n=1}^{N}d(\overline{Y},Y_n)\hfill \end{array}$$

As the first term in the last expression does not contain $\overline{Y}$,

$${\text{argmax}}_{\overline{Y}}\sum _{n=1}^{N}logp(Y_n|\overline{Y},\sigma )={\text{argmin}}_{\overline{Y}}\sum _{n=1}^{N}d(\overline{Y},Y_n)$$

The quantity on the right is exactly ${\widehat{Y}}_N$ by Equation (11). This proves the first claim. Now, consider the function:

$$F(\eta )=-Nlog({\zeta }_m(\frac{-1}{\eta }))+\eta \sum _{n=1}^{N}d({\widehat{Y}}_N,Y_n),\eta <\frac{-1}{{\sigma }_m}$$

This function is strictly concave, since it is the logarithm of the moment generating function of a positive measure. Note that ${lim}_{\eta \to \frac{-1}{{\sigma }_m}}F(\eta )=-\infty$, and admit for a moment that ${lim}_{\eta \to -\infty }F(\eta )=-\infty$. By the strict concavity of F, there exists a unique ${\widehat{\eta }}_N<\frac{-1}{{\sigma }_m}$ (which is the maximum of F), such that $F^{\prime }({\widehat{\eta }}_N)=0$. It follows that ${\widehat{\sigma }}_N=\frac{-1}{{\widehat{\eta }}_N}$ lies in $(0,{\sigma }_m)$ and satisfies Equation (20). The uniqueness of ${\widehat{\sigma }}_N$ is a consequence of the uniqueness of ${\widehat{\eta }}_N$. Thus, the proof is complete. Now, it remains to check that ${lim}_{\eta \to -\infty }F(\eta )=-\infty$ or just ${lim}_{\sigma \to +\infty }\frac{1}{\sigma }log({\zeta }_m(\frac{1}{\sigma }))=0$. Clearly:

$$\prod _{i<j}sinh\left(\frac{|r_i-r_j|}{2}\right)\le A_m{e}^{B_m|r|}$$

where $A_m$ and $B_m$ are two constants only depending on m. Using this, it follows that:

$$\frac{1}{\sigma }log({\zeta }_m(\frac{1}{\sigma }))\le \frac{1}{\sigma }log(c_m{A}_m)+\frac{1}{\sigma }log\left({\int }_{{\mathbb{R}}^m}exp((-\sigma +B_m)|r|)d{r}_1\cdots d{r}_m\right)$$

(21)

However, for some constant $C_m$ only depending on m,

$${\int }_{{\mathbb{R}}^m}exp((-\sigma +B_m)|r|)d{r}_1\cdots d{r}_m=C_m{\int }_0^{\infty }exp((-\sigma +B_m)u)u^{m-1}du$$

$$\le (m-1)!C_m{\int }_0^{\infty }exp((-\sigma +B_m+1)u)du=\frac{(m-1)!C_m}{\sigma -B_m-1}$$

Combining this bound and Equation (21) yields ${lim}_{\sigma \to +\infty }\frac{1}{\sigma }log({\zeta }_m(\frac{1}{\sigma }))=0$. ☐

Remark 1. 

Replacing F in the previous proof with $F(\eta )=-log({\zeta }_m(\frac{-1}{\eta }))+\eta c$ where $c>0$ shows that the equation:

$${\sigma }^2\times \frac{d}{d\sigma }log{\zeta }_m(\sigma )=c$$

has a unique solution $\sigma \in (0,{\sigma }_m)$. This shows in particular that $\sigma \mapsto {\sigma }^2\times \frac{d}{\phantom{(}d\sigma }log{\zeta }_m(\sigma )$ is a bijection from $(0,{\sigma }_m)$ to $(0,\infty )$.

Consider now the numerical computation of the maximum likelihood estimates ${\widehat{Y}}_N$ and ${\widehat{\sigma }}_N$ given by Proposition 3. Computation of ${\widehat{Y}}_N$ consists in finding the Riemannian median of $Y_1,\dots ,Y_N$, defined by Equation (11). This can be done using the Riemannian sub-gradient descent algorithm of [8]. The k-th iteration of this algorithm produces an approximation ${\widehat{Y}}_N^{{}^k}$ of ${\widehat{Y}}_N$ in the following way.

For $k=1,2,\dots$, let ${\Delta }_k$ be the symmetric matrix:

$${\Delta }_k=\frac{1}{N}\sum _{n=1}^N\frac{{\mathrm{Log}}_{{\widehat{Y}}_N^{{}^{k-1}}}(Y_n)}{||{\mathrm{Log}}_{{\widehat{Y}}_N^{{}^{k-1}}}(Y_n)||}$$

(22)

Here, $\mathrm{Log}$ is the Riemannian logarithm mapping inverse to the the Riemannian exponential mapping:

$${\mathrm{Exp}}_Y(\Delta )=Y^{1/2}exp\left(Y^{-1/2}\Delta Y^{-1/2}\right)Y^{1/2}$$

(23)

and $||{\mathrm{Log}}_a(b)||=\sqrt{g_a(b,b)}$. Then, ${\widehat{Y}}_N^{{}^k}$ is defined to be:

$${\widehat{Y}}_N^{{}^k}={\mathrm{Exp}}_{{\widehat{Y}}_N^{{}^{k-1}}}({\tau }_k{\Delta }_k)$$

(24)

where ${\tau }_k>0$ is a step size, which can be determined using a backtracking procedure.

Computation of ${\widehat{\sigma }}_N$ requires solving a non-linear equation in one variable. This is readily done using Newton’s method.

It is shown now that the empirical Riemannian median ${\widehat{Y}}_N$ converges almost surely to the true median $\overline{Y}$. This means that ${\widehat{Y}}_N$ is a consistent estimator of $\overline{Y}$. The proof of this fact requires few notations and a preparatory lemma.

For $\overline{Y}\in {\mathcal{P}}_m$ and $\sigma \in (0,{\sigma }_m)$, let:

$$\mathcal{E}(Y|\overline{Y},\sigma )={\int }_{{\mathcal{P}}_m}d(Y,Z)p(Z|\overline{Y},\sigma )dv(Z)$$

The following lemma shows how to find $\overline{Y}$ and σ from the function $Y\mapsto \mathcal{E}(Y|\overline{Y},\sigma )$.

Lemma 4. 

For any $\overline{Y}\in {\mathcal{P}}_m$ and $\sigma \in (0,{\sigma }_m)$, the following properties hold

(i) 

$\overline{Y}$ is given by:

$$\overline{Y}={\mathrm{argmin}}_Y\mathcal{E}(Y|\overline{Y},\sigma )$$

(25a)

That is, $\overline{Y}$ is the Riemannian median of $\mathcal{L}(\overline{Y},\sigma )$.

(ii) 

σ is given by:

$$\sigma =\Phi \left(\mathcal{E}(\overline{Y}|\overline{Y},\sigma )\right)$$

(25b)

where the function Φ is the inverse function of $\sigma \mapsto {\sigma }^2\times dlog{\zeta }_m(\sigma )/d\sigma$.

Proof of Lemma 4. 

(i) Let $\mathcal{E}(Y)=\mathcal{E}(Y|\overline{Y},\sigma )$. According to Theorem 2.1 in [28], this function has a unique global minimum, which is also a unique stationary point. Thus, to prove that $\overline{Y}$ is the minimum point of $\mathcal{E}$, it will suffice to check that for any geodesic γ starting from $\overline{Y}$, $\frac{d}{dt}{|}_{t=0}\mathcal{E}(\gamma (t))=0$ [31] (p. 76). Note that:

$$\frac{d}{dt}{|}_{t=0}\mathcal{E}(\gamma (t))={\int }_{{\mathcal{P}}_m}\frac{d}{dt}{|}_{t=0}d(\gamma (t),Z)p(Z|\overline{Y},\sigma )dv(Z)$$

(26)

where for all $Z\ne \overline{Y}$ [32]:

$$\frac{d}{dt}{|}_{t=0}d(\gamma (t),Z)=-g_{\overline{Y}}({log}_{\overline{Y}}(Z),{\gamma }^{\prime }(0))d{(\overline{Y},Z)}^{-1}$$

The integral in Equation (26) is, up to a constant,

$$\frac{d}{dt}{|}_{t=0}{\int }_{{\mathcal{P}}_m}p(Z|\gamma (t),\sigma )dv(Z)=0$$

since ${\int }_{{\mathcal{P}}_m}p(Z|\gamma (t),\sigma )dv(Z)=1$.

(ii) Differentiating ${\int }_{{\mathcal{P}}_m}exp(-\frac{d(Z,\overline{Y})}{\sigma })dv(Z)={\zeta }_m(\sigma )$ with respect to σ, it comes that:

$${\sigma }^2\times dlog{\zeta }_m(\sigma )/d\sigma ={\sigma }^2\frac{{\zeta }_m^{\prime }(\sigma )}{{\zeta }_m(\sigma )}={\int }_{{\mathcal{P}}_m}d(Z,\overline{Y})p(Z|\overline{Y},\sigma )dv(Z)=\mathcal{E}(\overline{Y}|\overline{Y},\sigma )$$

which proves (ii). ☐

Proposition 5 

(Consistency of ${\widehat{Y}}_N$). Let $Y_1,Y_2,\cdots$ be independent samples from a Laplace distribution $G(\overline{Y},\sigma )$. The empirical median ${\widehat{Y}}_N$ of $Y_1,\dots ,Y_N$ converges almost surely to $\overline{Y}$, as $N\to \infty$.

Proof of Proposition 5. 

Corollary 3.5 in [33] (p. 49) states that if $(Y_n)$ is a sequence of i.i.d. random variables on ${\mathcal{P}}_m$ with law π, then the Riemannian median ${\widehat{Y}}_N$ of $Y_1,\cdots ,Y_N$ converges almost surely as $N\to \infty$ to ${\widehat{Y}}_{\pi }$, the Riemannian median of π defined by Equation (12). Applying this result to $\pi =\mathcal{L}(\overline{Y},\sigma )$ and using ${\widehat{Y}}_{\pi }=\overline{Y}$, which follows from item (i) of Lemma 4, shows that ${\widehat{Y}}_N$ converges almost surely to $\overline{Y}$. ☐

4. Mixtures of Laplace Distributions

There are several motivations for considering mixtures of distributions in general. The most natural approach is to envisage a dataset as constituted of several subpopulations. Another approach is the fact that there is a support for the argument that mixtures of distributions provide a good approximation to most distributions in a spirit similar to wavelets.

The present section introduces the class of probability distributions that are finite mixtures of Riemannian Laplace distributions on ${\mathcal{P}}_m$. These constitute the main theoretical tool, to be used for the target application of the present paper, namely the problem of texture classification in computer vision, which will be treated in Section 5.

A mixture of Riemannian Laplace distributions is a probability distribution on ${\mathcal{P}}_m$, whose density with respect to the Riemannian volume element Equation (13) has the following expression:

$$p(Y|{({\varpi }_{\mu },{\overline{Y}}_{\mu },{\sigma }_{\mu })}_{1\le \mu \le M})=\sum _{\mu =1}^M{\varpi }_{\mu }\times p(Y|{\overline{Y}}_{\mu },{\sigma }_{\mu })$$

(27)

where ${\varpi }_{\mu }$ are nonzero weights, whose sum is equal to one, ${\overline{Y}}_{\mu }\in {\mathcal{P}}_m$ and ${\sigma }_{\mu }\in (0,{\sigma }_m)$ for all $1\le \mu \le M$, and the parameter M is called the number of mixture components.

Section 4.1 describes a new EM algorithm, which computes the maximum likelihood estimates of the mixture parameters ${({\varpi }_{\mu },{\overline{Y}}_{\mu },{\sigma }_{\mu })}_{1\le \mu \le M}$, based on independent observations $Y_1,\dots ,Y_N$ from the mixture distribution Equation (27).

Section 4.2 considers the problem of order selection for mixtures of Riemannian Laplace distributions. Precisely, this consists of finding the number M of mixture components in Equation (27) that realizes the best representation of a given set of data $Y_1,\dots ,Y_N$. This problem is solved by computing the BIC criterion, which is here found in explicit form for the case of mixtures of Riemannian Laplace distributions on ${\mathcal{P}}_m$.

4.1. Estimation of the Mixture Parameters

In this section, $Y_1,\dots ,Y_N$ are i.i.d. samples from Equation (27). Based on these observations, an EM algorithm is proposed to estimate ${({\varpi }_{\mu },{\overline{Y}}_{\mu },{\sigma }_{\mu })}_{1\le \mu \le M}$. The derivation of this algorithm can be carried out similarly to [15].

To explain how this algorithm works, define for all $\vartheta =\{({\varpi }_{\mu },{\overline{Y}}_{\mu },{\sigma }_{\mu })\}$,

$${\omega }_{\mu }(Y_n,\vartheta )=\frac{{\varpi }_{\mu }\times p(Y_n|{\overline{Y}}_{\mu },{\sigma }_{\mu })}{{\sum }_{s=1}^M{\varpi }_s\times p(Y_n|{\overline{Y}}_s,{\sigma }_s)},N_{\mu }(\vartheta )=\sum _{n=1}^N{\omega }_{\mu }(Y_n)$$

(28)

The algorithm iteratively updates $\widehat{\vartheta }=\{({\widehat{\varpi }}_{\mu },{\widehat{Y}}_{\mu },{\widehat{\sigma }}_{\mu })\}$, which is an approximation of the maximum likelihood estimate of the mixture parameters $\vartheta =({\varpi }_{\mu },{\overline{Y}}_{\mu },{\sigma }_{\mu })$ as follows.

• Update for ${\widehat{\varpi }}_{\mu }$: Based on the current value of $\widehat{\vartheta }$, assign to ${\widehat{\varpi }}_{\mu }$ the new value ${\widehat{\varpi }}_{\mu }=N_{\mu }(\widehat{\vartheta })/N.$
• Update for ${\widehat{Y}}_{\mu }$: Based on the current value of $\widehat{\vartheta }$, assign to ${\widehat{Y}}_{\mu }$ the value:

  $${\widehat{Y}}_{\mu }={\text{argmin}}_Y\sum _{n=1}^N{\omega }_{\mu }(Y_n,\widehat{\vartheta })d(Y,Y_n)$$

  (29)
• Update for ${\widehat{\sigma }}_{\mu }$: Based on the current value of $\widehat{\vartheta }$, assign to ${\widehat{\sigma }}_{\mu }$ the new value:

  $${\widehat{\sigma }}_{\mu }=\Phi (N_{\mu }^{-1}(\widehat{\vartheta })\times {\sum }_{n=1}^N{\omega }_{\mu }(Y_n,\widehat{\vartheta })d({\widehat{Y}}_{\mu },Y_n))$$

  (30)

where the function Φ is defined in Proposition 4.

These three update rules should be performed in the above order. Realization of the update rules for ${\widehat{\varpi }}_{\mu }$ and ${\widehat{\sigma }}_{\mu }$ is straightforward. The update rule for ${\widehat{Y}}_{\mu }$ is realized using a slight modification of the sub-gradient descent algorithm described in Section 3.2. More precisely, the factor $1/N$ appearing in Equation (22) is only replaced with ${\omega }_{\mu }(Y_n,\widehat{\vartheta })$ at each iteration.

In practice, the initial conditions $({\widehat{\varpi }}_{{\mu }_0},{\widehat{Y}}_{{\mu }_0},{\widehat{\sigma }}_{{\mu }_0})$ in this algorithm were chosen in the following way. The weights $({\varpi }_{{\mu }_0})$ are uniform and equal to $1/M$; $({\widehat{Y}}_{{\mu }_0})$ are M different observations from the set $\{Y_1,..,Y_N\}$ chosen randomly; and $({\widehat{\sigma }}_{{\mu }_0})$ is computed from $({\varpi }_{{\mu }_0})$ and $({\widehat{Y}}_{{\mu }_0})$ according to the rule Equation (30). Since the convergence of the algorithm depends on the initial conditions, the EM algorithm is run several times, and the best result is retained, i.e., the one maximizing the log-likelihood function.

4.2. The Bayesian Information Criterion

The BIC was introduced by Schwarz to find the appropriate dimension of a model that will fit a given set of observations [16]. Since then, BIC has been used in many Bayesian modeling problems where priors are hard to set precisely. In large sample settings, the fitted model favored by BIC ideally corresponds to the candidate model that is a posteriori most probable; i.e., the model that is rendered most plausible by the data at hand. One of the main features of the BIC is its easy computation, since it is only based on the empirical log-likelihood function.

Given a set of observations $\{Y_1,\cdots ,Y_N\}$ arising from Equation (27) where M is unknown, the BIC consists of choosing the parameter:

$$\overline{M}={\text{argmax}}_{M}BIC(M)$$

where:

$$BIC(M)=LL-\frac{1}{2}\times DF\times log(N)$$

(31)

Here, $LL$ is the log-likelihood given by:

$$LL=\sum _{n=1}^{N}log\left(\sum _{k=1}^M{\widehat{\varpi }}_{k}p(Y_n|{\widehat{Y}}_k,{\widehat{\sigma }}_k)\right)$$

(32)

and $DF$ is the number of degrees of freedom of the statistical model:

$$DF=M\times \frac{m(m+1)}{2}+M+M-1$$

(33)

In Formula Equation (32), ${({\widehat{\varpi }}_k,{\widehat{Y}}_k,{\widehat{\sigma }}_k)}_{1\le k\le M}$ are obtained from an EM algorithm as stated in Section 4.1 assuming the exact dimension is M. Finally, note that in Formula Equation (33), $M\times \frac{m(m+1)}{2}$ (respectively M and $M-1$) corresponds to the number of degrees of freedom associated with ${({\widehat{Y}}_k)}_{1\le k\le M}$ (respectively ${({\widehat{\sigma }}_k)}_{1\le k\le M}$ and ${({\widehat{\varpi }}_k)}_{1\le k\le M}$).

5. Application to Classification of Data on ${\mathcal{P}}_m$

Recently, several approaches have used the Riemannian distance in general as the main innovation in image or signal classification problems [2,15,34]. It turns out that the use of this distance leads to more accurate results (in comparison, for example, with the Euclidean distance). This section proposes an application that follows a similar approach, but in addition to the Riemannian distance, it also relies on a statistical approach. It considers the application of the Riemannian Laplace distribution (RLD) to the classification of data in ${\mathcal{P}}_m$ and gives an original Laplace classification rule, which can be used to carry out the task of classification, even in the presence of outliers. It also applies this classification rule to the problem of texture classification in computer vision, showing that it leads to improved results in comparison with recent literature.

Section 5.1 considers, from the point of view of statistical learning, the classification of data with values in ${\mathcal{P}}_m$. Given data points $Y_1,\cdots ,Y_N\in {\mathcal{P}}_m$, this proceeds in two steps, called the learning phase and the classification phase, respectively. The learning phase uncovers the class structure of the data, by estimating a mixture model using the EM algorithm developed in Section 4.1. Once training is accomplished, data points are subdivided into disjoint classes. Classification consists of associating each new data point to the most suitable class. For this, a new classification rule will be established and shown to be optimal.

Section 5.2 is the implementation of the Laplace classification rule together with the BIC criterion to texture classification in computer vision. It highlights the advantage of the Laplace distribution in the presence of outliers and shows its better performance compared to recent approaches.

5.1. Classification Using Mixtures of Laplace Distributions

Assume to be given a set of training data $Y_1,\cdots ,Y_N$. These are now modeled as a realization of a mixture of Laplace distributions:

$$p(Y)=\sum _{\mu =1}^M{\varpi }_{\mu }\times p(Y|{\overline{Y}}_{\mu },{\sigma }_{\mu })$$

(34)

In this section, the order M in Equation (34) is considered as known. The training phase of these data consists of learning its structure as a family of M disjoint classes $C_{\mu },\mu =1,\cdots ,M$. To be more precise, depending on the family $({\varpi }_{\mu })$, some of these classes may be empty. Training is done by applying the EM algorithm described in Section 4.1. As a result, each class $C_{\mu }$ is represented by a triple $({\widehat{\varpi }}_{\mu },{\widehat{Y}}_{\mu },{\widehat{\sigma }}_{\mu })$ corresponding to maximum likelihood estimates of $({\varpi }_{\mu },Y_{\mu },{\sigma }_{\mu })$. Each observation $Y_n$ is now associated with the class $C_{{\mu }^{*}}$ where ${\mu }^{*}={\text{argmax}}_{\mu }\omega (Y_n,\widehat{\nu })$ (recall the definition from Equation (28)). In this way, $\{Y_1,\cdots ,Y_N\}$ is subdivided into M disjoint classes.

The classification phase requires a classification rule. Following [15], the optimal rule (in the sense of a Bayesian risk criterion given in [35]) consists of associating any new data $Y_t$ to the class $C_{{\mu }^{*}}$ where:

$${\mu }^{*}={\text{argmax}}_{\mu }\left\{{\widehat{N}}_{\mu }\times p(Y_t|{\widehat{Y}}_{\mu },{\widehat{\sigma }}_{\mu })\right\}$$

(35)

Here, ${\widehat{N}}_{\mu }$ is the number of elements in $C_{\mu }$. Replacing ${\widehat{N}}_{\mu }$ with $N\times {\widehat{\varpi }}_{\mu }$, Equation (35) becomes ${\text{argmax}}_{\mu }{\widehat{\varpi }}_{\mu }\times p(Y_t|{\widehat{Y}}_{\mu },{\widehat{\sigma }}_{\mu })$. Note that when the weights ${\varpi }_{\mu }$ in Equation (34) are assumed to be equal, this rule reduces to a maximum likelihood classification rule ${max}_{\mu }p(Y_t|{\widehat{Y}}_{\mu },{\widehat{\sigma }}_{\mu })$. A quick look at the expression Equation (17) shows that Equation (35) can also be expressed as:

$${\mu }^{*}={\text{argmin}}_{\mu }\left\{-log{\widehat{\varpi }}_{\mu }+log\zeta ({\widehat{\sigma }}_{\mu })+\frac{d(Y_t,{\widehat{Y}}_{\mu })}{{\widehat{\sigma }}_{\mu }}\right\}$$

(36)

The rule Equation (36) will be called the Laplace classification rule. It favors clusters $C_{\mu }$ having a larger number of data points (the minimum contains $-log{\widehat{\varpi }}_{\mu }$) or a smaller dispersion away from the median (the minimum contains $log\zeta ({\widehat{\sigma }}_{\mu })$). When choosing between two clusters with the same number of points and the same dispersion, this rule favors the one whose median is closer to $Y_t$. If the number of data points inside clusters and the respective dispersions are neglected, then Equation (36) reduces to the nearest neighbor rule involving only the Riemannian distance introduced in [2].

The analogous rules of Equation (36) for the Riemannian Gaussian distribution (RGD) [15] and the Wishart distribution (WD) [17] on ${\mathcal{P}}_m$ can be established by replacing $p(Y_t|{\widehat{Y}}_{\mu },{\widehat{\sigma }}_{\mu })$ in Equation (35) with the RGD and the WD and then following the same reasoning as before. Recall that a WD depends on an expectation $\Sigma \in {\mathcal{P}}_m$ and a number of degrees of freedom n [29]. For the WD, Equation (36) becomes:

$${\mu }^{*}={\text{argmin}}_{\mu }\left\{-2log\widehat{\varpi }(\mu )-\widehat{n}(\mu )(logdet({\widehat{\Sigma }}^{-1}(\mu )Y_t)-\mathrm{tr}({\widehat{\Sigma }}^{-1}(\mu )Y_t))\right\}$$

Here, $\widehat{\varpi }(\mu )$, $\widehat{\Sigma }(\mu )$ and $\widehat{n}(\mu )$ denote maximum likelihood estimates of the true parameters $\varpi (\mu )$, $\Sigma (\mu )$ and $n(\mu )$, which define the mixture model (these estimates can be computed as in [36,37]).

5.2. Application to Texture Classification

This section presents an application of the mixture of Laplace distributions to the context of texture classification on the MIT Vision Texture (VisTex) database [38]. The purpose of this experiment is to classify the textures, by taking into consideration the within-class diversity. In addition, the influence of outliers on the classification performances is analyzed. The obtained results for the RLD are compared to those given by the RGD [15] and the WD [17].

The VisTex database contains 40 images, considered as being 40 different texture classes. The database used for the experiment is obtained after several steps. First of all, each texture is decomposed into 169 patches of $128\times 128$ pixels, with an overlap of 32 pixels, giving a total number of 6760 textured patches. Next, some patches are corrupted, in order to introduce abnormal data into the dataset. Therefore, their intensity is modified by applying a gradient of luminosity. For each class, between zero and 60 patches are modified in order to become outliers. An example of a VisTex texture with one of its patches and an outlier patch are shown in Figure 1.

Once the database is built, it is 15-times equally and randomly divided in order to obtain the training and the testing sets that are further used in the supervised classification algorithm. Then, for each patch in both databases, a feature vector has to be computed. The luminance channel is first extracted and then normalized in intensity. The grayscale patches are filtered using the stationary wavelet transform Daubechies db4 filter (see [39]), with two scales and three orientations. To model the wavelet sub-bands, various stochastic models have been proposed in the literature. Among them, the univariate generalized Gaussian distribution has been found to accurately model the empirical histogram of wavelet sub-bands [40]. Recently, it has been proposed to model the spatial dependency of wavelet coefficients. To this aim, the wavelet coefficients located in a $p\times q$ spatial neighborhood of the current spatial position are clustered in a random vector. The realizations of these vectors can be further modeled by elliptical distributions [41,42], copula-based models [43,44], etc. In this paper, the wavelet coefficients are considered as being realizations of zero-mean multivariate Gaussian distributions. In addition, for this experiment the spatial information is captured by using a vertical ($2\times 1$) and a horizontal ($1\times 2$) neighborhood. Next, the $2\times 2$ sample covariance matrices are estimated for each wavelet sub-band and each neighborhood. Finally, each patch is represented by a set of $F=12$ covariance matrices (2 scales × 3 orientations × 2 neighborhoods) denoted $Y=[Y_1,\cdots ,Y_F]$.

The estimated covariance matrices are elements of ${\mathcal{P}}_m$, with $m=2$, and therefore, they can be modeled by Riemannian Laplace distributions. More precisely, in order to take into consideration the within-class diversity, each class in the training set is viewed as a realization of a mixture of Riemannian Laplace distributions (Equation (27)) with M mixture components, characterized by $({\varpi }_{\mu },{\overline{Y}}_{\mu ,f},{\sigma }_{\mu ,f})$, having ${\overline{Y}}_{\mu ,f}\in {\mathcal{P}}_2$, with $\mu =1,\cdots ,M$ and $f=1,\cdots ,F$. Since the sub-bands are assumed to be independent, the probability density function is given by:

$$p(Y|{({\varpi }_{\mu },{\overline{Y}}_{\mu ,f},{\sigma }_{\mu ,f})}_{1\le \mu \le M,1\le f\le F})=\sum _{\mu =1}^M{\varpi }_{\mu }\prod _{f=1}^{F}p(Y_f|{\overline{Y}}_{\mu ,f},{\sigma }_{\mu ,f})$$

(37)

The learning step of the classification is performed using the EM algorithm presented in Section 4, and the number of mixture components is determined using the BIC criterion recalled in Equation (31). Note that for the considered model given in Equation (37), the degree of freedom is expressed as:

$$DF=M-1+M\times F\times \left(\frac{m(m+1)}{2}+1\right)$$

(38)

since one centroid and one dispersion parameter should be estimated per feature and per component of the mixture model. In practice, the number of mixture components M varies between two and five, and the M yielding to the highest BIC criterion is retained. As mentioned earlier, the EM algorithm is sensitive to the initial conditions. In order to minimize this influence, for this experiment, the EM algorithm is repeated 10 times, and the result maximizing the log-likelihood function is retained. Finally, the classification is performed by assigning each element $Y_t\in {\mathcal{P}}_2$ in the testing set to the class of the closest cluster ${\mu }^{*}$, given by:

$${\mu }^{*}={\text{argmin}}_{\mu }\left\{-log{\widehat{\varpi }}_{\mu }+\sum _{f=1}^{F}log\zeta ({\widehat{\sigma }}_{\mu ,f})+\sum _{f=1}^F\frac{d(Y_t,{\widehat{Y}}_{\mu ,f})}{{\widehat{\sigma }}_{\mu ,f}}\right\}$$

(39)

This expression is obtained starting from Equations (36) and (37), knowing that F features are extracted for each patch.

The classification results of the proposed model (solid red line), expressed in terms of overall accuracy, shown in Figure 2, are compared to those given by a fixed number of mixture components (that is, for $M=3$, dashed red line) and with those given when the within-class diversity is not considered (that is, for $M=1$, dotted red line). In addition, the classification performances given by the RGD model (displayed in black) proposed in [15] and the WD model (displayed in blue) proposed in [17] are also considered. For each of these models, the number of mixture components is first computed using the BIC, and next, it is fixed to $M=3$ and $M=1$. For all of the considered models, the classification rate is given as a function of the number of outliers, which varies between zero and 60 for each class.

It is shown that, as the number of outliers increases, the RLD gives progressively better results than the RGD and the WD. The results are improved by using the BIC criterion for choosing the suitable number of clusters. In conclusion, the mixture of RLDs combined with the BIC criterion to estimate the best number of mixture components can minimize the influence of abnormal samples present in the dataset, illustrating the relevance of the proposed method.

6. Conclusions

Motivated by the problem of outliers in statistical data, this paper introduces a new distribution on the space ${\mathcal{P}}_m$ of $m\times m$ symmetric positive definite matrices, called the Riemannian Laplace distribution. Denoted throughout the paper by $\mathcal{L}(\overline{Y},\sigma )$, where $\overline{Y}\in {\mathcal{P}}_m$ and $\sigma >0$ are the indexing parameters, this distribution may be thought of as specifying the law of a family of observations on ${\mathcal{P}}_m$ concentrated around the location $\overline{Y}$ and having dispersion σ. If d denotes Rao’s distance on ${\mathcal{P}}_m$ and $dv(Y)$ its associated volume form, the density of $\mathcal{L}(\overline{Y},\sigma )$ with respect to $dv(Y)$ is proportional to $exp(-\frac{d(Y,\overline{Y}}{\sigma }))$. Interestingly, the normalizing constant depends only on σ (and not on $\overline{Y}$). This allows us to deduce exact expressions for maximum likelihood estimates of $\overline{Y}$ and σ relying on the Riemannian median on ${\mathcal{P}}_m$. These estimates are also computed numerically by means of sub-gradient algorithms. The estimation of parameters in mixture models of Laplace distributions are also considered and performed using a new expectation-maximization algorithm. Finally, the main theoretical results are illustrated by an application to texture classification. The proposed experiment consists of introducing abnormal data (outliers) into a set of images from the VisTex database and analyzing their influences on the classification performances. Each image is characterized by a set of $2\times 2$ covariance matrices modeled as mixtures of Riemannian Laplace distributions in the space ${\mathcal{P}}_2$. The number of mixtures is estimated using the BIC criterion. The obtained results are compared to those given by the Riemannian Gaussian distribution, showing the better performance of the proposed method.