Kernel Density Estimation on the Siegel Space with an Application to Radar Processing ^†

Abstract

This paper studies probability density estimation on the Siegel space. The Siegel space is a generalization of the hyperbolic space. Its Riemannian metric provides an interesting structure to the Toeplitz block Toeplitz matrices that appear in the covariance estimation of radar signals. The main techniques of probability density estimation on Riemannian manifolds are reviewed. For computational reasons, we chose to focus on the kernel density estimation. The main result of the paper is the expression of Pelletier’s kernel density estimator. The computation of the kernels is made possible by the symmetric structure of the Siegel space. The method is applied to density estimation of reflection coefficients from radar observations.

1. Introduction

Various techniques can be used to estimate the density of probability measure in the Euclidean spaces, such as histograms, kernel methods, or orthogonal series. These methods can sometimes be adapted to densities in Riemannian manifolds. The computational cost of the density estimation depends on the isometry group of the manifold. In this paper, we study the special case of the Siegel space. The Siegel space is a generalization of the hyperbolic space. It has a structure of symmetric Riemannian manifold, which enables the adaptation of different density estimation methods at a reasonable cost. Convergence rates of the density estimation using kernels and orthogonal series were gradually generalized to Riemannian manifolds (see [1,2,3]).

The Siegel space appears in radar processing in the study of Toeplitz block Toeplitz matrices, whose blocks represent covariance matrices of a radar signal (see [4,5,6]). The Siegel also appears in statistical mechanics, see [7] and was recently used in image processing (see [8]). Information geometry is now a standard framework in radar processing (see [4,5,6,9,10,11,12,13]). The information geometry on positive definite Teoplitz block Teoplitz matrices is directly related to the metric on the Siegel space (see [14]). Indeed, Toeplitz block Toeplitz matrices can be represented by a symmetric positive definite matrix and a point laying in a product of Siegel disks. The metric considered on Toeplitz block Toeplitz matrices is induced by the product metric between a metric on the symmetric positive definite matrices and the Siegel disks metrics (see [4,5,6,9,14]).

One already encounters the problem of density estimation in the hyperbolic space for electrical impedance [15], networks [16] and radar signals [17]. In [18], a generalization of the Gaussian law on the hyperbolic space was proposed. Apart from [19], where authors propose a generalization of the Gaussian law, probability density estimation on the Siegel space has not yet been addressed.

The contributions of the paper are the following. We review the main non parametric density estimation techniques on the Siegel disk. We provide some rather simple explicit expressions of the kernels defined by Pelletier in [1]. These expressions make the kernel density estimation the most adapted method. We present visual results of estimated densities in the simple case where the Siegel disk reduces to the Poincaré disk.

The paper begins with an introduction to the Siegel space in Section 2. Section 3 reviews the main non-parametric density estimation techniques on the Siegel space. Section 3.3 contains the original results of the paper. Section 4 presents an application to radar data estimation.

2. The Siegel Space

This section presents facts about the Siegel space. The interested reader can find more details in [20,21]. The necessary background on Lie groups and symmetric space can be found in [22].

2.1. The Siegel Upper Half Space

The Siegel upper half space is a generalization of the Poincaré upper half space (see [23]) for a description of the hyperbolic space. Let $Sym\left(n\right)$ be the space of real symmetric matrices of size $n\times n$ and $Sy{m}_{+}\left(n\right)$ the set of real symmetric positive definite matrices of size $n\times n$. The Siegel upper half space is defined by

$${\mathcal{H}}_n=\left\{Z=X+iY|X\in Sym\left(n\right),Y\in Sy{m}_{+}\left(n\right)\right\}.$$

${\mathcal{H}}_n$ is equipped with the following metric:

$$ds=2tr\left(Y^{-1}dZ{Y}^{-1}d\overline{Z}\right).$$

The set of real symplectic matrices $Sp(n,\mathbb{R})$ is defined by

$$g\in Sp(n,\mathbb{R})\iff g^{t}Jg=J,$$

where

$$J=\left(\begin{array}{cc}0& I_n\\ -I_n& 0\end{array}\right),$$

and $I_n$ is the $n\times n$ identity matrix. $Sp(n,\mathbb{R})$ is a subgroup of $S{L}_{2n}\left(\mathbb{R}\right)$, the set of $2n\times 2n$ invertible matrices of determinant 1. Let $g=\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\in Sp(n,\mathbb{R})$. The metric $ds$ is invariant under the following action of $Sp(n,\mathbb{R})$,

$$g.Z=(AZ+B){(CZ+D)}^{-1}.$$

This action is transitive, i.e.,

$$\forall Z\in {\mathcal{H}}_n,\exists g\in Sp(n,\mathbb{R}),g.iI=Z.$$

The stabilizer K of $iI$ is the set of elements g of $Sp(n,\mathbb{R})$ whose action leaves $iI$ fixed. K is a subgroup of $Sp(n,\mathbb{R})$ called the isotropy group. We can verify that

$$K=\left\{\left(\begin{array}{cc}A& B\\ -B& A\end{array}\right),A+iB\in SU\left(n\right)\right\}.$$

A symmetric space is a Riemannian manifold, where the reversal of the geodesics is well defined and is an isometry. Formally, $ex{p}_p\left(u\right)\mapsto ex{p}_p(-u)$ is an isometry for each p on the manifold, where u is a vector in the tangent space at p, and $ex{p}_p$ the Riemannian exponential application at p. In other words, the symmetry around each point is an isometry. ${\mathcal{H}}_n$ is a symmetric space (see [20]). The structure of a symmetric space can be studied through its isometry group and the Lie algebra of its isometry group. The present work will make use of the Cartan and Iwasawa decompositions of the Lie algebra of $Sp(n,\mathbb{R})$ (see [22]). Let $\mathfrak{sp}(n,\mathbb{R})$ be the Lie algebra of $Sp(n,\mathbb{R})$. Given A, B and C three real $n\times n$ matrices, let denote $\left(\begin{array}{cc}A& B\\ C& -A^t\end{array}\right)=(A,B,C)$. We have

$$\mathfrak{sp}(n,\mathbb{R})=\left\{(A,B,C)|B\mathrm{and}C\mathrm{symmetric}\right\}.$$

The Cartan decomposition of $\mathfrak{sp}(n,\mathbb{R})$ is given by

$$\mathfrak{sp}(n,\mathbb{R})=\mathfrak{t}\oplus \mathfrak{p},$$

where

$$\mathfrak{t}=\left\{(A,B,-B)|B\mathrm{symmetric}\mathrm{and}A\text{skew-symmetric}\right\},$$

$$\mathfrak{p}=\left\{(A,B,B)|A,B,\mathrm{symmetric}\right\}.$$

(1)

The Iwasawa decomposition is given by

$$\mathfrak{sp}(n,\mathbb{R})=\mathfrak{t}\oplus \mathfrak{a}\oplus \mathfrak{n},$$

where

$$\mathfrak{a}=\left\{(H,0,0)|H\mathrm{diagonal}\right\},$$

$$\mathfrak{n}=\left\{(A,B,0)|A\mathrm{upper}\mathrm{triangular}\mathrm{with}0\mathrm{on}\mathrm{the}\mathrm{diagonal},B\mathrm{symmetric}\right\}.$$

It can be shown that

$$\mathfrak{p}={\cup }_{k\in K}A{d}_k\left(\mathfrak{a}\right),$$

(2)

where $Ad$ is the adjoint representation of $Sp(n,\mathbb{R})$.

2.2. The Siegel Disk

The Siegel disk ${\mathcal{D}}_n$ is the set of complex matrices $\left\{Z|I-Z^{*}Z\ge 0\right\}$, where ≥ stands for the Loewner order (see [24] for details on the Loewner order). Recall that for A and B two Hermitian matrices, $A\ge B$ with respect to the Loewner order means that $A-B$ is positive definite. The transformation

$$Z\in {\mathcal{H}}_n\mapsto (Z-iI){(Z+iI)}^{-1}\in {\mathcal{D}}_n$$

is an isometry between the Siegel upper half space and the Siegel disk. Let $C=\left(\begin{array}{cc}I& -iI\\ I& iI\end{array}\right)$. The application $g\in Sp(n,\mathbb{R})\mapsto Cg{C}^{-1}$ identifies the set of isometries of ${\mathcal{H}}_n$ and of ${\mathcal{D}}_n$. Thus, it can be shown that a matrix $g=\left(\begin{array}{cc}A& B\\ \overline{A}& \overline{B}\end{array}\right)\in Sp(n,\mathbb{C})$ acts isometrically on ${\mathcal{D}}_n$ by

$$g.Z=(AZ+B){(\overline{A}Z+\overline{B})}^{-1},$$

where $\overline{A}$ stands for the conjugate of A. The point $iI$ in ${\mathcal{H}}_n$ is identified with the null matrix noted 0 in ${\mathcal{D}}_n$. Let $Z\in {\mathcal{D}}_n$. There exists P a diagonal matrix with decreasing positive real entries and U a unitary matrix such that $Z=UP{U}^t$. Let ${\tau }_1\ge ...\ge {\tau }_n$ be such that

$$P=\left(\begin{array}{cc}th\left({\tau }_1\right)& \\ & \ddots & \\ & & th\left({\tau }_n\right)\end{array}\right).$$

Let

$$A_0=\left(\begin{array}{cc}ch\left({\tau }_1\right)& \\ & \ddots & \\ & & ch\left({\tau }_n\right)\end{array}\right),B_0=\left(\begin{array}{cc}sh\left({\tau }_1\right)& \\ & \ddots & \\ & & sh\left({\tau }_n\right)\end{array}\right)$$

and

$$g_Z=\left(\begin{array}{cc}U& 0\\ 0& \overline{U}\end{array}\right).\left(\begin{array}{cc}{A}_0& B_0\\ A_0& B_0\end{array}\right).$$

It can be shown that

$$g_Z\in Sp(n,\mathbb{C})\mathrm{and}{g}_Z.0=Z.$$

(3)

We provide now a correspondence between the elements of ${\mathcal{D}}_n$ and the elements of $\mathfrak{p}$ defined in Equation (1). Let

$$H_Z=\left(\begin{array}{cc}{\tau }_1& \\ & \ddots & \\ & & {\tau }_n\\ & & & -{\tau }_1& \\ & & & & \ddots & \\ & & & & & -{\tau }_n\end{array}\right)\in \mathfrak{a},$$

(4)

and

$$a_Z=\left(\begin{array}{cc}{e}^{{\tau }_1}& \\ & \ddots & \\ & & e^{{\tau }_n}\\ & & & e^{-{\tau }_1}& \\ & & & & \ddots & \\ & & & & & e^{-{\tau }_n}\end{array}\right)\in A=exp\left(\mathfrak{a}\right).$$

It can be shown that there exists $k\in K$ such that

$$Cexp\left(A{d}_k\left(H_Z\right)\right)C^{-1}.0=Z,$$

or equivalently

$$Ck{a}_{Z}k{C}^{-1}.0=Z.$$

Recall that Equation (2) gives $A{d}_k\left(H\right)\in \mathfrak{p}$ and $kak\in exp\left(\mathfrak{p}\right)$. The distance between Z and 0 in ${\mathcal{D}}_n$ is given by

$$d(0,Z)={\left(2\sum {\tau }_i^2\right)}^{1/2}$$

(5)

(see p. 292 in [20] ).

3. Non Parametric Density Estimation on the Siegel Space

Let Ω be a space, endowed with a σ-algebra and a probability measure p. Let X be a random variable $\Omega \to {\mathcal{D}}_n$. The Riemannian measure of ${\mathcal{D}}_n$ is called $vol$ and the measure on ${\mathcal{D}}_n$ induced by X is noted ${\mu }_X$. We assume that ${\mu }_X$ has a density, noted f, with respect to $vol$, and that the support of X is a compact set noted $Supp$. Let $(x_1,...,x_k)\in {\mathcal{D}}_n^k$ be a set of draws of X.

The Dirac measure at a point $a\in {\mathcal{D}}_n$ is denoted ${\delta }_a$. Let ${\mu }_k=\frac{1}{k}{\sum }_{i=1}^k{\delta }_{x_i}$ denotes the empirical measure of the set of draws. This section presents four non-parametric techniques of estimation of the density f from the set of draws $(x_1,...,x_k)$. The estimated density at x in ${\mathcal{D}}_n$ is noted ${\widehat{f}}_k\left(x\right)=\widehat{f}(x,x_1,...,x_k)$. The relevance of a density estimation technique depends on several aspects. When the space allows it, the estimation technique should equally consider each direction and location. This leads to an isotropy and a homogeneity condition. In the kernel method, for instance, a kernel density function $K_{x_i}$ is placed at each observation $x_i$. Firstly, in order to treat directions equally, the function $K_{x_i}$ should be invariant under the isotropy group of $x_i$; Secondly, for another observation $x_j$, functions $K_{x_i}$ and $K_{x_j}$ should be similar up to the isometries that send $x_i$ on $x_j$. These considerations strongly depend on the geometry of the space: if the space is not homogeneous and the isotropy group is empty, these indifference principles have no meaning. Since the Siegel space is symmetric, it is homogeneous and has a non empty isotropy group. Thus, the density estimation technique should be chosen accordingly.

The convergence of the different estimation techniques is widely studied. Results were first obtained in the Euclidean case, and are gradually extended to the probability densities on manifold (see [1,2,15,25]).

The last relevant aspect is computational. Each estimation technique has its own computational framework that presents pros and cons given the different applications. For instance, the estimation by orthogonal series needs an initial pre-processing, but provides a fast evaluation of the estimated density in compact manifolds.

3.1. Histograms

The histogram is the simplest density estimation method. Given a partition of the space ${\mathcal{D}}_n={\cup }_i{A}_i$, the estimated density is given by

$$\widehat{f}(x\in A_i)=\frac{1}{vol\left(A_i\right)}\sum _{j=1}^k{1}_{A_i}\left(x_j\right),$$

where $1_{A_i}$ stands for the indicator function of $A_i$. Following the considerations of the previous sections, the elements of the partition should firstly be as isotropic as possible, and secondly as similar as possible to each other. Regarding the problem of histograms, the case of the Siegel space is similar to the case of the hyperbolic space. There exist various uniform polygonal tilings on the Siegel space that could be used to compute histograms. However, there are ratio $\lambda \in \mathbb{R}$ for which there is no homothety. Thus, it is not always possible to adapt the size of the bins to a given set of draws of the random variable. Modifying the size of the bins can require a change of the structure of the tiling. This is why the study of histograms has not been deepened.

3.2. Orthogonal Series

The estimation of the density f can be made out of the estimation of the scalar product between f and a set of “orthonormal” functions $\left\{e_j\right\}$. The most standard choice for $\left\{e_j\right\}$ is the eigenfunctions of the Laplacian. When the variable X takes its values in ${\mathbb{R}}^n$, this estimation technique becomes the characteristic function method. When the underlying space is compact, the spectrum of the Laplacian operator is countable, while when the space is non-compact, the spectrum is uncountable. In the first case, the estimation of the density f is made through the estimation of a sum, while in the second case is made through the estimation of an integral. In practice, the second situation presents a larger computational complexity. Unfortunately, eigenfunctions of the Laplacian operator are known on ${\mathcal{D}}_n$ but not on compact sub-domains. This is why the study of this method has not been deepened.

3.3. Kernels

Let $\mathcal{K}:{\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$ be a map which verifies the following properties:

(i)

${\int }_{{\mathbb{R}}^d}\mathcal{K}\left(\right|\left|x\right|\left|\right)dx=1$;

(ii)

${\int }_{{\mathbb{R}}^d}x\mathcal{K}\left(\right|\left|x\right|\left|\right)dx=0$;

(iii)

$\mathcal{K}(x>1)=0$;

(iv)

$sup\left(\mathcal{K}\right(x\left)\right)=\mathcal{K}\left(0\right)$.

Let $p\in {\mathcal{D}}_n$. Generally, given a point p on a Riemannian manifold, $ex{p}_p$ defines an injective application only on a neighborhood of 0. On the Siegel space, $ex{p}_p$ is injective on the whole space. When the tangent space $T_p{\mathcal{D}}_n$ is endowed with the local scalar product,

$$\left|\right|u\left|\right|=d(p,ex{p}_p\left(u\right)),$$

where $\left|\right|.\left|\right|$ is the Euclidean distance associated with the local scalar product and $d(.,.)$ is the Riemannian distance. The corresponding Lebesgue measure on $T_p{\mathcal{D}}_n$ is noted $Le{b}_p$. Let $ex{p}_p^{*}\left(Le{b}_p\right)$ denote the push-forward measure of $Le{p}_p$ by $ex{p}_p$. The function ${\theta }_p$ defined by:

$${\theta }_p:q\mapsto {\theta }_p\left(q\right)=\frac{dvol}{dex{p}_p^{*}\left(Le{b}_p\right)}\left(q\right)$$

(6)

is the density of the Riemannian measure on ${\mathcal{D}}_n$ with respect to the Lebesgue measure $Le{b}_p$ after the identification of ${\mathcal{D}}_n$ and $T_p{\mathcal{D}}_n$ induced by $ex{p}_p$ (see Figure 1).

Given $\mathcal{K}$ and a positive radius r, the estimator of f proposed by [1] is defined by:

$${\widehat{f}}_k=\frac{1}{k}\sum _i\frac{1}{r^n}\frac{1}{{\theta }_{x_i}\left(x\right)}\mathcal{K}\left(\frac{d(x,x_i)}{r}\right).$$

(7)

The corrective factor ${\theta }_{x_i}{\left(x\right)}^{-1}$ is necessary since the kernel $\mathcal{K}$ originally integrates to one with respect to the Lebesgue measure and not with respect to the Riemannian measure. It can be noticed that this estimator is the usual kernel estimator in the case of Euclidean space. When the curvature of the space is negative, which is the case of the Siegel space, the distribution placed over each sample $x_i$ has $x_i$ as intrinsic mean. The following theorem provides convergence rate of the estimator. It is a direct adaptation of Theorem 3.1 of [1].

Theorem 1.

Let $(\mathcal{M},g)$ be a Riemannian manifold of dimension n and μ its Riemannian volume measure. Let X be a random variable taking its values in a compact subset C of $(\mathcal{M},g)$. Let $0<r\le r_{inj}$, where $r_{inj}$ is the infimum of the injectivity radius on C. Assume the law of X has a twice differentiable density f with respect to the Riemannian volume measure. Let ${\widehat{f}}_k$ be the estimator defined in Equation (7). There exists a constant $C_f$ such that

$${\int }_{x\in \mathcal{M}}{E}_{x_1,...,x_k}[{(f\left(x\right)-{\widehat{f}}_k\left(x\right))}^2]d\mu \le C_f(\frac{1}{k{r}^n}+r^4).$$

(8)

If $r\sim k^{\frac{-1}{n+4}}$,

$${\int }_{x\in \mathcal{M}}{E}_{x_1,...,x_k}[{(f\left(x\right)-{\widehat{f}}_k\left(x\right))}^2]d\mu =O(k^{-\frac{4}{n+4}}).$$

(9)

Proof. 

See Appendix A. ☐

It can be checked that on the Siegel space $r_{inj}=+\infty$ and that, for an isometry α, we have:

$${\widehat{f}}_k(x,x_1,...,x_k)={\widehat{f}}_k(\alpha \left(x\right),\alpha \left(x_1\right),...,\alpha \left(x_k\right)).$$

Each location and direction are processed as similarly as possible. This density estimator can be used for data classification on Riemannian manifolds, see [26].

In order to obtain the explicit expression of the estimator, one must have the explicit expression of the Riemannian exponential, of its inverse, and of the function ${\theta }_p$ (see Equations (6) and (7)). These expressions are difficult and sometimes impossible to obtain for general Riemannian manifolds. In the case of the Siegel space, the symmetric structure makes the computation possible. Since the space is homogeneous, the computation can be made at the origin $iI\in {\mathcal{H}}_n$ or $0\in {\mathcal{D}}_n$ and transported to the whole space. In the present work, the random variable lays in ${\mathcal{D}}_n$. However, in the literature, the Cartan and Iwasawa decompositions are usually given for the isometry group of ${\mathcal{H}}_n$. Thus, our computation starts in ${\mathcal{H}}_n$ before moving to ${\mathcal{D}}_n$.

The Killing form on the Lie algebra $\mathfrak{sp}(n,\mathbb{R})$ of the isometry group of ${\mathcal{H}}_n$ induces a scalar product on $\mathfrak{p}$. This scalar product can be transported on $exp\left(\mathfrak{p}\right)$ by left multiplication. This operation gives $exp\left(\mathfrak{p}\right)$ a Riemannian structure. It can be shown that on this Riemannian manifold, the Riemannian exponential at the identity coincides with the group exponential. Furthermore,

$$\begin{array}{ccccc}\varphi & :& exp\left(\mathfrak{p}\right)& \to & {\mathcal{H}}_n\\ & & g& \mapsto & g.iI\end{array}$$

(10)

is a bijective isometry, up to a scaling factor. Since the volume change ${\theta }_p$ is invariant under rescaling of the metric, this scaling factor has no impact. Thus, ${\mathcal{H}}_n$ can be identified with $exp\left(\mathfrak{p}\right)$ and $ex{p}_{iI\in {\mathcal{H}}_n}$ with $ex{p}_{|\mathfrak{p}}$. The expression of the Riemannian exponential is difficult to obtain in general; however, it boils down to the group exponential in the case of symmetric spaces. This is the main element of the computation of ${\theta }_p$. The Riemannian volume measure on $exp\left(\mathfrak{p}\right)$ is noted $vo{l}^{\prime }$. Let

$$\begin{array}{ccccc}\psi & :& K\times \mathfrak{a}& \to & \mathfrak{p}\\ & & (k,H)& \mapsto & A{d}_k\left(H\right).\end{array}$$

Let ${\mathfrak{a}}^{+}$ be the diagonal matrices with strictly decreasing positive eigenvalues. Let ${\Lambda }^{+}$ be the set of positive roots of $\mathfrak{sp}(n,\mathbb{R})$ as described in p. 282 in [20],

$${\Lambda }^{+}=\{e_i+e_j,i\le j\}\cup \{e_i-e_j,i<j\},$$

where $e_i\left(H\right)$ is the i-th diagonal term of the diagonal matrix H. Let $C_c\left(E\right)$ be the set of continuous compactly supported functions on the space E. In [27], at page 73, it is given that for all $t\in C_c\left(\mathfrak{p}\right)$, there exists $c_1>0$ such that

$${\int }_{\mathfrak{p}}t\left(Y\right)dY=c_1{\int }_K{\int }_{{\mathfrak{a}}^{+}}t\left(\psi (k,H)\right)\prod _{\lambda \in {\Lambda }^{+}}\lambda \left(H\right)dkdH,$$

(11)

where $dY$ is a Lebesgue measure on the coefficients of Y. Let $\tilde{\mathfrak{p}}=\psi (K\times {\mathfrak{a}}^{+})$. $\lambda \in {\Lambda }^{+}$ never vanishes on ${\mathfrak{a}}^{+}$ and $\mathfrak{p}\backslash \tilde{\mathfrak{p}}$ has a null measure. Thus,

$${\int }_{\tilde{\mathfrak{p}}}t\left(Y\right)\prod _{\lambda \in {\Lambda }^{+}}\frac{1}{\lambda \left(H_Y\right)}dY=c_1{\int }_K,{\int }_{{\mathfrak{a}}^{+}}t\left(A{d}_k\left(H\right)\right)dkdH,$$

(12)

where $H_Y$ is the point in ${\mathfrak{a}}^{+}$ such that there exists k in K such that $\psi (k,H_Y)=Y$. Calculation in p. 73 in [27] also gives that for all $t\in C_c\left(\mathfrak{p}\right)$, there exists $c_2>0$, such that

$${\int }_{Sp(n,\mathbb{R})}t\left(g\right)dg=c_2{\int }_K{\int }_{{\mathfrak{a}}^{+}}{\int }_{K}t(k_2.exp\left(A{d}_{k_1}\left(H\right)\right))J\left(H\right)d{k}_{1}dHd{k}_2,$$

(13)

where $dg$ is the Haar measure on $Sp(n,\mathbb{R})$ and

$$\begin{array}{ccc}\hfill J\left(H\right)& =& \prod _{\lambda \in {\Lambda }^{+}}{e}^{\lambda \left(H\right)}-e^{-\lambda \left(H\right)}\hfill \\ \hfill & =& 2^{|{\Lambda }^{+}|}\prod _{\lambda \in {\Lambda }^{+}}sinh\left(\lambda \left(H\right)\right).\hfill \end{array}$$

Thus, for all $t\in C_c(Sp(n,\mathbb{R})/K)$,

$${\int }_{Sp(n,\mathbb{R})/K}t\left(x\right)dx=c_2{\int }_K{\int }_{{\mathfrak{a}}^{+}}t\left(exp\left(A{d}_k\left(H\right)\right)\right)J\left(H\right)dkdH,$$

(14)

where $dx$ is the invariant measure on $Sp(n,\mathbb{R})/K$. After identifying $Sp(n,\mathbb{R})/K$ and $exp\left(\mathfrak{p}\right)$, the Riemannian measure on $exp\left(\mathfrak{p}\right)$ coincides with the invariant measure on $Sp(n,\mathbb{R})/K$. Thus, for all $t\in C_c\left(exp\left(\mathfrak{p}\right)\right)$,

$${\int }_{exp\left(\mathfrak{p}\right)}t\left(x\right)dvo{l}^{\prime }=c_2{\int }_K{\int }_{{\mathfrak{a}}^{+}}t\left(exp\left(A{d}_k\left(H\right)\right)\right)J\left(H\right)dkdH.$$

(15)

Using the notation $H_Y$ of Equation (12),

$${\int }_{\tilde{\mathfrak{p}}}t\left(exp\left(Y\right)\right)J\left(H_Y\right)\prod _{\lambda \in {\Lambda }^{+}}\frac{1}{\lambda \left(H_Y\right)}dY=c_1{\int }_K{\int }_{{\mathfrak{a}}^{+}}t\left(exp\left(A{d}_k\left(H\right)\right)\right)J\left(H\right)dkdH.$$

(16)

Combining Equations (15) and (16), we obtain that there exists $c_3$ such that

$${\int }_{\tilde{\mathfrak{p}}}t\left(exp\left(Y\right)\right)\prod _{\lambda \in {\Lambda }^{+}}\frac{sinh\left(\lambda \left(H_Y\right)\right)}{\lambda \left(H_Y\right)}dY=c_3{\int }_{exp\left(\mathfrak{p}\right)}t\left(x\right)dvo{l}^{\prime }.$$

(17)

The term $\frac{sinh\left(\lambda \right(H\left)\right)}{\lambda \left(H\right)}$ can be extended by continuity on $\mathfrak{a}$; thus,

$${\int }_{\mathfrak{p}}t\left(exp\left(Y\right)\right)\prod _{\lambda \in {\Lambda }^{+}}\frac{sinh\left(\lambda \left(H_Y\right)\right)}{\lambda \left(H_Y\right)}dY=c_3{\int }_{exp\left(\mathfrak{p}\right)}t\left(x\right)dvo{l}^{\prime }.$$

(18)

Let $dY$ be the Lebesgue measure corresponding to the metric. Then, the exponential application does not introduce a volume change at $0\in \mathfrak{p}$. Since $H_0=0$ and $\frac{sinh\left(\lambda \right(H\left)\right)}{\lambda \left(H\right)}\underset{H\to 0}{⟶}1$, we have $c_3=1$. Let $log$ denote the inverse of the exponential application. We have

$$\frac{dlo{g}^{*}\left(vo{l}^{\prime }\right)}{dY}=\prod _{\lambda \in {\Lambda }^{+}}\frac{sinh\left(\lambda \left(H_Y\right)\right)}{\lambda \left(H_Y\right)}.$$

Since ϕ from Equation (10) is an isometry up to a scaling factor, if $Y\in \mathfrak{p}$ and $C\varphi \left(exp\left(Y\right)\right)C^{-1}=ex{p}_0(u\in T_0{\mathcal{D}}_n)$, then

$$\frac{dlo{g}^{*}\left(vol\right)}{dLe{b}_0}\left(u\right)=\frac{dlo{g}^{*}\left(vo{l}^{\prime }\right)}{dY}\left(Y\right),$$

where $Le{b}_0$ refers to the Lebesgue measure on the tangent space $T_0{\mathcal{D}}_n$ as in Equation (6). Given $Z\in {\mathcal{D}}_n$, $H_Z$ from Equation (4) verifies $C\varphi \left(exp\left(A{d}_k\left(H_Z\right)\right)\right)C^{-1}=Z$ for some k in K. Thus,

$${\theta }_0\left(Z\right)=\frac{dlo{g}^{*}\left(vo{l}^{\prime }\right)}{dY}\left(A{d}_k\left(H_Z\right)\right)=\prod _{\lambda \in {\Lambda }^{+}}\frac{sinh\left(\lambda \left(H_Z\right)\right)}{\lambda \left(H_Z\right)}.$$

We have then

$${\theta }_0\left(Z\right)=\prod _{i<j}\frac{sinh({\tau }_i-{\tau }_j)}{{\tau }_i-{\tau }_j}\prod _{i\le j}\frac{sinh({\tau }_i+{\tau }_j)}{{\tau }_i+{\tau }_j},$$

where the $\left({\tau }_i\right)$ are described in Section 2.2. Given $Z_1,Z_2\in {\mathcal{D}}_n$,

$${\theta }_{Z_1}\left(Z_2\right)={\theta }_0(g_{Z_1}^{-1}.Z_2),$$

where $g_{Z_1}^{-1}$ is defined in Equation (3). It is thus possible to use the density estimator defined in Equation (7). Indeed,

$$\frac{1}{{\theta }_{Z_1}\left(Z_2\right)}\mathcal{K}\left(\frac{d(Z_1,Z_2)}{r}\right)=\prod _{i<j}\frac{{\tau }_i-{\tau }_j}{sinh({\tau }_i-{\tau }_j)}\prod _{i\le j}\frac{{\tau }_i+{\tau }_j}{sinh({\tau }_i+{\tau }_j)}\mathcal{K}\left(\frac{{(2\sum {\tau }_i^2)}^{1/2}}{r}\right),$$

(19)

where the $\left({\tau }_i\right)$ are the diagonal elements of $H_{g_{Z_1}^{-1}.Z_2}$. Recall that when $n=1$, the Siegel disk corresponds to the Poincaré disk. Thus, we retrieve the expression of the kernel for the hyperbolic space,

$$\frac{1}{{\theta }_{Z_1}\left(Z_2\right)}\mathcal{K}\left(\frac{d(Z_1,Z_2)}{r}\right)=\frac{2\tau }{sinh\left(2\tau \right)}\mathcal{K}\left(\frac{{\left(2{\tau }^2\right)}^{1/2}}{r}\right).$$

(20)

4. Application to Radar Processing

4.1. Radar Data

In space time adaptative radar processing (STAP), the signal is formed by a succession of matrices X representing the realization of a temporal and spatial process. Let ${\mathcal{B}}_{n,m}^{+}$ be the set of positive definite block Teoplitz matrices composed of $n\times n$ blocks of $m\times m$ matrices (PD BT). For a stationary signal, the autocorrelation matrix R is PD BT (see [5,6,14]). Authors of [5,6,14] proposed a generalization of Verblunsky coefficients and defined a parametrization of PD BT matrices,

$$\begin{array}{ccc}{\mathcal{B}}_{n,m}^{+}& \to & Sy{m}^{+}\times {\mathbb{D}}_n^{m-1}\\ \mathcal{R}& \mapsto & (P_0,{\Omega }_1,...,{\Omega }_{m-1}),\end{array}$$

(21)

in which the metric induced by the Kähler potential is the product metric of an affine invariant metric on $Sy{m}^{+}$ and $m-1$ times the metric of the Siegel disk, up to a scaling factor. When the signal is not Gaussian, reflection/Verblunsky coefficients in Poincaré or Siegel Disks should be normalized as described in [28] by a normalized Burg algorithm. Among other references, positive definite block Teoplitz matrices have been studied in the context of STAP-radar processing in [4,5,6].

4.2. Marginal Densities of Reflection Coefficients

In this section, we show density estimation results of the marginal parameters ${\Omega }_k$. For the sake of visualization, only the Siegel disk ${\mathbb{D}}_1$ is considered. Recall that ${\mathbb{D}}_1$ coincides with the Poincaré disk. The results are partly extracted from the conference paper [17]. Data used in the experimental tests are radar observations from THALES X-band Radar, recorded during 2014 field trials campaign at Toulouse Blagnac Airport for European FP7 UFO study (Ultra-Fast wind sensOrs for wake-vortex hazards mitigation) (see [29,30]). Data are representative of Turbulent atmosphere monitored by radar. Figure 2 illustrates the density estimation of six coefficients on the Poincaré unit disk under a rainy environment. The densities are individually re-scaled for visualization purposes. For each environment, the dataset is composed of 120 draws. The densities of the coefficients ${\Omega }_k$ are representative of the background. This information on the background is expected to ease the detection of interesting targets.

4.3. Radar Clutter Segmentation

Clutter refers to background Doppler signal related to meteorological conditions (e.g., wind in wooded areas, currents and breaking waves on water), which hinders detection of small and slow targets. At each range, a set of reflection coefficients are computed from the Doppler spectrum (see [31]). This set of coefficients is a point in the Poincaré poly-disk. From this set of points in the poly-disk, it is possible to estimate the underlying density. Segmenting clutter, i.e., determining zones of homogeneous Doppler characteristics (see Figure 3), enables the improvement of detection algorithms on each zone. The mean-shift algorithm enables segmentation of the space according to the kernel density estimation of a set of points. It was introduced by Fukunaga and Hostetler in 1975 (see [32]). It corresponds to a gradient ascent of the density estimator (see [33]) for a study of the statistical consistency of the gradient lines estimation. Each data point moves to a local mode of the density estimator, which yields as many clusters as modes. This algorithm has been generalized on manifolds in [34], and applied to radar images in [35]. It can thus be used to segment the set of points in the Poincaré poly-disk. Unfortunately, the mean-shift algorithm requires working with a kernel depending only on the distance to its barycenter, which is not the case of the kernel defined in Equation (19). Thus, the computations are performed without the use of the corrective term ${\theta }_p$. It is possible to solve this problem by replacing the corrective term by its average at a given radius, which leads to a kernel depending only on the distance to its barycenter. Our future work will focus on the computation of these averages. Let

$${\widehat{f}}_r^{\mathcal{K}}\left(x\right)=\frac{c_d}{k}\sum _{i=1}^k\frac{1}{r^n}\mathcal{K}\left(\frac{d{(x_i,x)}^2}{r^2}\right),$$

where $c_n$ is a normalization constant. Let $g=-k^{\prime }$.

The mean-shift is defined by

$$m\left(x\right)=\sum _{i=1}^k\frac{\frac{1}{r^{n+2}}g\left(\frac{d{(x,x_i)}^2}{r^2}\right)}{{\sum }_{i=1}^k\frac{1}{r^{n+2}}g\left(\frac{d{(x,x_i)}^2}{r^2}\right)}lo{g}_x\left(x_i\right)\propto \frac{\nabla f_r^{\mathcal{K}}}{f_r^g},$$

where $m\left(x\right)$ is in the tangent space at x. The algorithm moves from x to $ex{p}_x\left(m\left(x\right)\right)$ until convergence to a local maximum. The points of the space are segmented according to the local maxima to which they converge.

In order to assess the quality of unsupervised classification, we use the notion of Silhouette, see [36], which computes for each point a proximity criterion with respect to other points of the same cluster and other points of different clusters (see Figure 4). Let x be in the cluster A. We respectively define $a\left(x\right)=mi{n}_{y\in A}d(x,y)$ and $b\left(x\right)=mi{n}_{y\ne A}d(x,y)$, the minimum distance to points of the same (resp. other) class(es). The Silhouette of x is

$$\frac{a\left(x\right)-b\left(x\right)}{max\left\{a\right(x),b(x\left)\right\}},$$

which takes values between $-1$ and 1, respectively, when the data point is considered “badly” and “well” clustered. The average of all the silhouettes provides an indication of the relevance of the classification. One can represent graphically the silhouette profile by plotting for each class horizontal segments of the length of the silhouette value (see Figure 5).

In order to test the Riemannian Mean Shift performance, we generate simple synthetic radar clutter data. Given 250 range cells, we generate 125 cells of ground clutter (wind) centered at 0 m·s⁻¹, of spectral width 5 m·s⁻¹, to which we add 125 cells of rain clutter, centered at 5 m·s⁻¹, of spectral width 10 m·s⁻¹. This clutter is sampled 10 times and the segmentation is performed on each simulation (see Figure 6, Figure 7 and Figure 8).

It can be seen that, apart from a few outliers, the two clutters are well classified and that the algorithm was able to distinguish between two zones of different Doppler characteristics.

We then test our algorithm on real sea clutter data (see Figure 9, Figure 10 and Figure 11).

The results are more difficult to interpret in that case. The Doppler spectra are varying quite a lot along the range axis. Even though it looks over-segmented, the first classification (kernel size defined by the distance to the 10th closest neighbor point) displays the highest average silhouette value.

5. Conclusions

Three non parametric density estimation techniques have been considered. The main advantage of histograms in the Euclidean context is their simplicity of use. This makes histograms an interesting tool despite the fact that they do not present optimal convergence rates. On the Siegel space, histograms lose their simplicity advantage. They were thus not deeply studied. The orthogonal series density estimation also presents technical disadvantages on the Siegel space. Indeed, the series become integrals, which make the numerical computation of the estimator more difficult than in the Euclidean case. On the other hand, the use of the kernel density estimator does not present major differences with the Euclidean case. The convergence rate obtained in [1] can be extended to compactly supported random variables on non compact Riemannian manifolds. Furthermore, the corrective term whose computation is required to use Euclidean kernels on Riemannian manifolds turns out to have a reasonably simple expression. Our future efforts will concentrate on the use of kernel density estimation on the Siegel space in radar signal processing. As the experimental section suggests, we strongly believe that the estimation of the densities of the ${\Omega }_k$ will provide an interesting description of the different backgrounds. This non-parametric method of density estimation should be compared with parametric ones, as “Maximum Entropy Density” (Gibbs density) on homogenesous manifold as proposed in [37] based on the works of Jean-Marie Souriau. As proposed in [38], a median-shift approach might also be investigated.