# Towards a median signal detector through the total Bregman divergence and
its robustness analysis

Yusuke Ono and Linyu Peng
Department of Mechanical Engineering, Keio University,
Yokohama 223-8522, Japan E-mail<EMAIL_ADDRESS>author.
E-mail<EMAIL_ADDRESS>
###### Abstract

A novel family of geometric signal detectors are proposed through medians of
the total Bregman divergence (TBD), which are shown advantageous over the
conventional methods and their mean counterparts. By interpreting the
observation data as Hermitian positive-definite matrices, their mean or median
play an essential role in signal detection. As is difficult to be solved
analytically, we propose numerical solutions through Riemannian gradient
descent algorithms or fixed-point algorithms.

Beside detection performance, robustness of a detector to outliers is also of
vital importance, which can often be analyzed via the influence functions.
Introducing an orthogonal basis for Hermitian matrices, we are able to compute
the corresponding influence functions analytically and exactly by solving a
linear system, which is transformed from the governing matrix equation.
Numerical simulations show that the TBD medians are more robust than their
mean counterparts.

Keywords: Geometric median, matrix-CFAR, robustness analysis, nonhomogeneous
clutter

##  1 Introduction

In target detection, the constant false alarm rate (CFAR) is a well-known and
conventionally used clutter suppression method [36]. A popular CFAR is based
on Fourier transform and makes use of the Doppler spectral density to
distinguish targets from clutter. However, this method is rather weak for
short waveform of dense and nonhomogeneous clutter. One of the popular
statistical methods for discriminating between the target signal from clutter
is to estimate the covariance matrix and test the likelihood ratio. There has
been a well-known method for estimating the clutter covariance matrix, namely,
the sample covariance matrix (SCM) using the maximum likelihood estimation
[20], but this method performs well only when enough number of signals is
independent and identically distributed [35]. In order to reduce the number of
required observation data for achieving efficient detection performance within
nonhomogeneous clutter, various knowledge-based methods using a priori
information were proposed [14, 29, 30, 39, 41]. In [4, 8], the Bayesian
methods were leveraged to estimate the covariance matrix. In practice, it is
hard to accurately capture the statistical information of clutter in advance
which can be range-dependent and non-stationary and the observation data is
usually contaminated. Insufficient knowledge of clutter environment can lead
to the detection performance degradation.

Recently, the matrix-CFAR, which needs no a priori information but makes use
of the covariance matrix of the signals, was proposed and developed [5, 6]. It
is based on the fact that the covariance matrix that represents the
correlation of complex signals is Hermitian positive-definite (HPD). It was
applied to target detection [42] and drone detection [10]. Its advantages have
been demonstrated in detection of high frequency X-band radar clutter [28],
Burg estimation of radar scatter matrices [16], monitoring of wake eddy
turbulence [7], etc. In recent years, the matrix-CFAR with respect to various
divergence functions has been proposed, and its detection performance and
robustness to outliers were shown better compared with that of the geodesic
distance associated to the affine-invariant Riemannian metric (AIRM) [22, 23,
24, 26].

In [25], the total Bregman divergence (TBD) defined in HPD manifolds was
applied to signal detection in nonhomogeneous clutter. The TBD means
associated to the total square loss (TSL), the total log-determinant (TLD)
divergence and the total von-Neumann (TVN) divergence were investigated; it
was shown numerically that the corresponding detectors outperform the AIRM
detector and generalized likelihood ratio test (GLRT) detector. In the current
study, to improve the robustness of the matrix-CFAR, we extend the previous
study [25] on TBD means to TBD medians as median detectors are often more
robust to outliers than the mean counterparts. In addition, a novel analytical
method for robustness analysis is proposed. Main contributions of the current
paper are outlined as follows.

  * (1)

We define the TBD medians associated with the TSL, the TLD divergence and the
TVN divergence, and investigate their detection performance compared with the
TBD means, the Riemannian distance (RD) mean and median associated to the
AIRM, the GLRT detector as well as the adaptive normalized matched filter
(ANMF). Medians are often difficult to be calculated analytically; instead, we
introduce fixed-point algorithms or Riemannian gradient descent algorithms to
solve them numerically.

  * (2)

Influence functions can be used to describe the robustness of a detector about
outliers. As noticed in [25], the matrix equations for influence functions are
not commonly seen in matrix theory and to solve them directly can be very
challenging; approximate solutions were obtained and applied therein. By
defining an orthogonal basis for Hermitian matrix spaces with the induced
Frobenius metric in this study, we transform those matrix equations to linear
systems, which can immediately (at least in principle) be solved exactly. This
methodology allows us to precisely capture the robustness of all means and
medians discussed in the current paper.

The structure of this paper is as follows. Theory of the matrix-CFAR is
summarized in Section 2. In Section 3, after briefly reviewing the geometry of
HPD equipped with the AIRM, we introduce Riemannian gradient descent
algorithms for solving the corresponding RD mean and median. Riemannian
gradient, also called natural gradient, assigns the steepest direction in
Riemannian manifolds. Furthermore, the TBD medians are defined and solved
numerically through the fixed-point algorithms. Numerical detection
performance is conduced in Section 4 by comparing the RD mean and median, the
TBD means and medians, and the conventional GLRT. In Section 5, an orthogonal
basis for Hermitian matrices is defined and applied to the computation of
influence functions, allowing us to accurately illustrate robustness of all
means and medians. We conclude and address several lines of potential future
researches finally in Section 6.

##  2 Formulation of the matrix-CFAR

In this section, the method of matrix-CFAR will be recalled. Assume that the
observation data are collected from $N$ stationary channels and consider a
detection problem concerning a moving target under the nonhomogeneous clutter
environment. The problem of detection will be modeled as a binary hypothesis
testing. In this testing, the null hypothesis $H_{0}$ represents that the
observation data $\bm{x}$ is only the clutter $\bm{c}$ and the alternative
hypothesis $H_{0}$ represents that $\bm{x}$ contains both the target signal
$\bm{s}$ and the clutter $\bm{c}$, namely

$\displaystyle\begin{cases}H_{0}:\begin{cases}\mathrm{target\ is\ absent,}\\\
\bm{x}=\bm{c},\\\ \end{cases}\vspace{0.2cm}\\\
H_{1}:\begin{cases}\mathrm{target\ is\ present,}\\\
\bm{x}=\xi\bm{s}+\bm{c},\\\ \end{cases}\\\ \end{cases}$ (2.1)

where $\xi$ denotes the unknown amplitude coefficient which is complex scalar-
valued. The observation data $\bm{x}$ and the target signal $\bm{s}$ denoting
the steering vector are respectively given by

$\displaystyle\bm{x}=\left(x_{0},\ldots,x_{N-1}\right)^{\operatorname{T}}\in\mathbb{C}^{N},$
(2.2)
$\displaystyle\bm{s}=\dfrac{1}{\sqrt{N}}\left(1,\exp(-\operatorname{i}2\pi
f_{d}),\ldots,\exp(-\operatorname{i}2\pi
f_{d}(N-1))\right)^{\operatorname{T}},$ (2.3)

where $\mathbb{C}^{N}$ represents the $N$ dimensional complex space,
$\operatorname{i}$ denotes the imaginary unit, $f_{d}$ is the signal
normalized Doppler frequency, and $(\cdot)^{\operatorname{T}}$ denotes the
transpose of matrices or vectors. Conventional SCM estimator, derived by
secondary data set $\left\\{\bm{x}_{i}\right\\}_{i=1}^{m}$, is given by [20]

$\displaystyle\bm{R}_{\mathrm{SCM}}=\frac{1}{m}\sum_{i=1}^{m}\bm{x}_{i}\bm{x}_{i}^{\mathsf{H}},\
\ \bm{x}_{i}\in\mathbb{C}^{N},$ (2.4)

where $(\cdot)^{\mathsf{H}}$ represents the conjugate transpose of matrices or
vectors. This estimator has been widely utilized; for instance, the GLRT
detector reads [27]

$\displaystyle\frac{\left|\bm{x}^{\mathsf{H}}\bm{R}^{-1}_{\mathrm{SCM}}\bm{s}\right|^{2}}{\bm{s}^{\mathsf{H}}\bm{R}^{-1}_{\mathrm{SCM}}\bm{s}}\overset{H_{1}}{\underset{H_{0}}{\gtrless}}\gamma$
(2.5)

and the ANMF reads [13]

$\displaystyle\frac{\left|\bm{s}^{\mathsf{H}}\bm{R}^{-1}_{\mathrm{SCM}}\bm{x}\right|^{2}}{\left(\bm{x}^{\mathsf{H}}\bm{R}^{-1}_{\mathrm{SCM}}\bm{x}\right)\left(\bm{s}^{\mathsf{H}}\bm{R}^{-1}_{\mathrm{SCM}}\bm{s}\right)}\overset{H_{1}}{\underset{H_{0}}{\gtrless}}\gamma,$
(2.6)

where $\gamma$ denotes the detection threshold.

In the matrix-CFAR, to determine whether the received signal $\bm{x}$ includes
a target $\bm{s}$ or not, we model the observation data by the covariance
matrix as a Toeplitz HPD matrix defined as [3, 9]

$\displaystyle\bm{R}=\begin{pmatrix}r_{0}&\cdots&\overline{r}_{k}&\cdots&\overline{r}_{N-1}\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
r_{k}&\ddots&r_{0}&\ddots&\overline{r}_{k}\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\ r_{N-1}&\cdots&r_{k}&\cdots&r_{0}\\\
\end{pmatrix}\in\mathscr{P}(N,\mathbb{C}),$ (2.7)

where $\mathscr{P}(N,\mathbb{C})$ denotes the space of $N\times N$ HPD
matrices. Its component

$r_{k}=\mathbb{E}\left[x_{l}\overline{x}_{l+k}\right],\ \ 0\leq k\leq N-1,\
0\leq l\leq N-l-1$ (2.8)

is the $k$-th correlation coefficient. Here, $\overline{r}_{k}$ is the
conjugate of $r_{k}$, and $\mathbb{E}[\cdot]$ is the statistical expectation.
The Toeplitz covariance structure has been widely studied and applied as it
can often improve the estimate of the covariance matrix; see, e.g., [2, 18,
43].

Based on the ergodicity of a stationary process, each component $r_{l}$ can be
approximated by the observation data as

$\displaystyle r_{k}$
$\displaystyle=\dfrac{1}{N}\sum\limits_{l=0}^{N-1-l}x_{l}\overline{x}_{l+k},\
\ 0\leq k\leq N-1.$ (2.9)

The covariance matrix of clutter is estimated as $\bm{R}_{g}$ by the
observations $\\{\bm{R}_{i}\\}_{i=1}^{m}$. Consequently, the detection problem
can be modeled as

$\displaystyle\begin{cases}H_{0}:\bm{R}_{\mathrm{CUT}}=\bm{R}_{g},\vspace{0.15cm}\\\
H_{1}:\bm{R}_{\mathrm{CUT}}\neq\bm{R}_{g},\\\ \end{cases}$ (2.10)

where the matrix $\bm{R}_{\mathrm{CUT}}$ is the covariance matrix of the
observation in the cell under test (CUT). By defining a threshold $\gamma$,
the signal detection is modeled as discriminating $\bm{R}_{\mathrm{CUT}}$ from
$\bm{R}_{g}$ by

$\displaystyle
d\left(\bm{R}_{g},\bm{R}_{\mathrm{CUT}}\right)\overset{H_{1}}{\underset{H_{0}}{\gtrless}}\gamma,$
(2.11)

where $d\left(\bm{R}_{g},\bm{R}_{\mathrm{CUT}}\right)$ is the difference
between $\bm{R}_{\mathrm{CUT}}$ and $\bm{R}_{g}$, such as the Riemannian
distance or TBD divergences that we are going to introduce in the next
section.

The matrix-CFAR illustrated by Fig. 1 is summarized up as follows. Estimator
of covariance matrices $\\{\bm{R}_{i}\\}_{i=1}^{m}$ of the observation data
$\\{\bm{x}_{i}\\}_{i=1}^{m}$ is firstly calculated by (2.7) and (2.9). Mean or
median $\bm{R}_{g}$ of those covariance matrices will be obtained by using
either the distance or divergences. The matrix $\bm{R}_{g}$ can be considered
as an estimator of the covariance matrix of the clutter since we consider
almost all observation data $\\{\bm{x}_{i}\\}_{i=1}^{m}$ as the clutter, that
is essential in design of the matrix-CFAR. The detection decision is performed
by the comparison of the distance or difference of $\bm{R}_{\mathrm{CUT}}$ and
$\bm{R}_{g}$ and the threshold $\gamma$.

Figure 1: The process of matrix-CFAR.

##  3 Means and medians of HPD matrices

In this section, we briefly introduce geometric means and medians in HPD
spaces equipped with either the AIRM as Riemannian manifolds or the Frobenius
metric as metric spaces. In the latter, divergence functions can be defined,
playing the role of measuring the difference of two HPD matrices.

###  3.1 HPD manifolds and the RD mean and median

We firstly define the AIRM on the HPD matrix manifold
$\mathscr{P}(N,\mathbb{C})$, and introduce the corresponding geodesic
distance, also called the Riemannian distance. The tangent space at a point
$\bm{P}\in\mathscr{P}(N,\mathbb{C})$ is
$T_{\bm{P}}\mathscr{P}(N,\mathbb{C})\cong\mathscr{H}(N,\mathbb{C})$, which
denotes the space of $N\times N$ Hermitian matrices. The tangent bundle is
denoted by
$T\mathscr{P}(N,\mathbb{C})=\cup_{\bm{P}}T_{\bm{P}}\mathscr{P}(N,\mathbb{C})$.
The well-known AIRM is defined by

$\displaystyle\langle\bm{A},\bm{B}\rangle_{\bm{P}}:=\operatorname{tr}\left(\bm{P}^{-1}\bm{A}\bm{P}^{-1}\bm{B}\right)$
(3.1)

for $\bm{A},\bm{B}\in T_{\bm{P}}\mathscr{P}(N,\mathbb{C})$, where
$\operatorname{tr}$ denotes the trace of a matrix. The geodesic curve between
two points $\bm{P}_{0},\bm{P}_{1}\in\mathscr{P}(N,\mathbb{C})$ is given by

$\bm{P}(t)=\bm{P}_{0}^{\frac{1}{2}}\exp\left(t\bm{P}_{0}^{\frac{1}{2}}\operatorname{Log}\left(\bm{P}_{0}^{-\frac{1}{2}}\bm{P}_{1}\bm{P}_{0}^{-\frac{1}{2}}\right)\bm{P}_{0}^{\frac{1}{2}}\right)\bm{P}_{0}^{\frac{1}{2}},\
t\in[0,1],$ (3.2)

where $\exp$ denotes matrix exponential and $\operatorname{Log}$ denotes the
principal logarithm [21]. Obviously, $\bm{P}(0)=\bm{P}_{0},$ and
$\bm{P}(1)=\bm{P}_{1}$. Under the AIRM metric, the geodesic distance between
two points $\bm{P}_{0},\bm{P}_{1}\in\mathscr{P}(N,\mathbb{C})$ is given by

$\displaystyle d_{R}\left(\bm{P}_{0},\bm{P}_{1}\right)$ $\displaystyle=$
$\displaystyle\left\lVert\operatorname{Log}\left(\bm{P}_{0}^{-\frac{1}{2}}\bm{P}_{1}\bm{P}_{0}^{-\frac{1}{2}}\right)\right\rVert$
(3.3) $\displaystyle=$
$\displaystyle\left\lVert\operatorname{Log}\left(\bm{P}_{1}^{-1}\bm{P}_{0}\right)\right\rVert,$

where $\left\|\cdot\right\|$ is the Frobenius norm associated to the Frobenius
metric
$\left\langle\bm{A},\bm{B}\right\rangle:=\operatorname{tr}{\left(\bm{A}^{\mathsf{H}}\bm{B}\right)}$.
For more details, the reader may refer to [33].

In the Riemannian manifold $\mathscr{P}(N,\mathbb{C})$ equipped with the AIRM,
the exponential map
$\operatorname{Exp}:T\mathscr{P}(N,\mathbb{C})\rightarrow\mathscr{P}(N,\mathbb{C})$
is defined by making use of the geodesics. Locally, assuming $\bm{P}(t)$ is
the geodesic (3.2) connecting two arbitrary HPD matrices
$\bm{P}_{0}=\bm{P}(0),\bm{P}_{1}=\bm{P}(1)$, the exponential map gives
endpoint of the geodesic, namely

$\displaystyle\operatorname{Exp}\left(\bm{P}_{0},\dot{\bm{P}}(0)\right)$
$\displaystyle=\bm{P}_{1}$
$\displaystyle=\bm{P}_{0}^{\frac{1}{2}}\exp\left(\bm{P}_{0}^{\frac{1}{2}}\operatorname{Log}\left(\bm{P}_{0}^{-\frac{1}{2}}\bm{P}_{1}\bm{P}_{0}^{-\frac{1}{2}}\right)\bm{P}_{0}^{\frac{1}{2}}\right)\bm{P}_{0}^{\frac{1}{2}}$
$\displaystyle=\bm{P}_{0}^{\frac{1}{2}}\exp\left(\bm{P}_{0}^{-\frac{1}{2}}\dot{\bm{P}}(0)\bm{P}_{0}^{-\frac{1}{2}}\right)\bm{P}_{0}^{\frac{1}{2}},$

where $\dot{\bm{P}}(0)$ is the tangent vector at the point $\bm{P}_{0}$. Using
the Riemannian distance (3.3), we recall the RD mean and median of $m$ HPD
matrices $\\{\bm{R}_{i}\\}_{i=1}^{m}$ as follows

$\displaystyle\overline{\bm{R}}_{\mathrm{RDmean}}$
$\displaystyle:=\underset{\bm{R}\in\mathscr{P}(N,\mathbb{C})}{\operatorname{arg}\min}\dfrac{1}{m}\sum\limits_{i=1}^{m}d_{R}^{2}(\bm{R},\bm{R}_{i}),$
(3.4) $\displaystyle\overline{\bm{R}}_{\mathrm{RDmedian}}$
$\displaystyle:=\underset{\bm{R}\in\mathscr{P}(N,\mathbb{C})}{\operatorname{arg}\min}\dfrac{1}{m}\sum\limits_{i=1}^{m}d_{R}(\bm{R},\bm{R}_{i}).$
(3.5)

In other words, the RD mean and median are, respectively, solutions of the
optimization problems with objective functions

$\displaystyle
G_{\mathrm{RDmean}}\left(\bm{R}\right)=\frac{1}{m}\sum\limits_{i=1}^{m}d_{R}^{2}(\bm{R},\bm{R}_{i}),$
(3.6) $\displaystyle
G_{\mathrm{RDmedian}}\left(\bm{R}\right)=\frac{1}{m}\sum\limits_{i=1}^{m}d_{R}(\bm{R},\bm{R}_{i}).$
(3.7)

######  Proposition 1.

Since the RD mean and median are related to optimisation problems on HPD
spaces as Riemannian manifolds, it is natural to propose the Riemannian
gradient descent algorithms, i.e., gradient descent algorithms on Riemannian
manifolds [11, 19]:

$\displaystyle\bm{R}_{t+1}$
$\displaystyle=\operatorname{Exp}\left(\bm{R}_{t},-\eta_{t}\operatorname{grad}G_{\mathrm{RDmean}}\left(\bm{R}_{t}\right)\right)$
(3.8)
$\displaystyle=\bm{R}_{t}^{\frac{1}{2}}\exp\left(-\eta_{t}\bm{R}_{t}^{-\frac{1}{2}}\operatorname{grad}G_{\mathrm{RDmean}}\left(\bm{R}_{t}\right)\bm{R}_{t}^{-\frac{1}{2}}\right)\bm{R}_{t}^{\frac{1}{2}}$
$\displaystyle=\bm{R}_{t}^{\frac{1}{2}}\exp\left(-\frac{2\eta_{t}}{m}\sum_{i=1}^{m}\operatorname{Log}\left(\bm{R}_{t}^{\frac{1}{2}}\bm{R}_{i}^{-1}\bm{R}_{t}^{\frac{1}{2}}\right)\right)\bm{R}_{t}^{\frac{1}{2}},$

and

$\displaystyle\bm{R}_{t+1}$
$\displaystyle=\operatorname{Exp}\left(\bm{R}_{t},-\eta_{t}\operatorname{grad}G_{\mathrm{RDmedian}}\left(\bm{R}_{t}\right)\right)$
(3.9)
$\displaystyle=\bm{R}_{t}^{\frac{1}{2}}\exp\left(-\eta_{t}\bm{R}_{t}^{-\frac{1}{2}}\operatorname{grad}G_{\mathrm{RDmedian}}\left(\bm{R}_{t}\right)\bm{R}_{t}^{-\frac{1}{2}}\right)\bm{R}_{t}^{\frac{1}{2}}$
$\displaystyle=\bm{R}_{t}^{\frac{1}{2}}\exp\left(-\frac{\eta_{t}}{m}\sum_{i=1}^{m}\frac{\operatorname{Log}\left(\bm{R}_{t}^{\frac{1}{2}}\bm{R}_{i}^{-1}\bm{R}_{t}^{\frac{1}{2}}\right)}{d_{R}\left(\bm{R}_{t},\bm{R}_{i}\right)}\right)\bm{R}_{t}^{\frac{1}{2}}.$

Here $\eta_{t}$ is the step size, and $\operatorname{grad}$ denotes the
Riemannian gradient on Riemannian manifolds.

For more details on Riemannian gradient descent algorithms, the reader may
refer to [37, 38]. The RD mean and median can also be calculated numerically
through the fixed-point algorithms [22]. It was noted that the detection
performance of the RD mean in the matrix-CFAR is better than that of the RD
median [22].

###  3.2 The total Bregman divergence, TBD means and medians

The Bregman divergence for matrices was introduced [17] and this idea has been
extended to total Bregman divergences [25]. The total Bregman divergence of
two matrices $\bm{A},\bm{B}\in\operatorname{GL}\left(N,\mathbb{C}\right)$, the
general linear group of $N\times N$ invertible complex-valued matrices, is
defined by

$\delta_{F}\left(\bm{A},\bm{B}\right)=\frac{F(\bm{A})-F(\bm{B})-\left\langle\nabla
F(\bm{B}),\bm{A}-\bm{B}\right\rangle}{\sqrt{1+\left\|\nabla
F(\bm{B})\right\|^{2}}},$ (3.10)

where $F\left(\bm{A}\right)$ is a differentiable and strictly convex function
in $\operatorname{GL}\left(N,\mathbb{C}\right)$ and $\nabla$ denotes the
(Euclidean) gradient with respect to the Frobenius metric, defined by

$\langle\nabla
F(\bm{A}),\bm{B}\rangle=\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}F(\bm{A}+\varepsilon\bm{B}),\
\ \bm{A},\bm{B}\in\operatorname{GL}(N,\mathbb{C}).$ (3.11)

The well-known total square loss (TSL), and the total log-determinant (TLD)
and total von-Neumann (TVN) divergences defined in the HPD manifold
$\mathscr{P}(N,\mathbb{C})$ are shown in Table 1.

Table 1: The TBDs | The function $F(\bm{Y})$ | TBDs in $\mathscr{P}(N,\mathbb{C})$
---|---|---
TSL | $\frac{1}{2}\left\|\bm{Y}\right\|^{2}$ |  $\delta_{\mathrm{TSL}}\left(\bm{Y},\bm{Z}\right)=\cfrac{\left\|\bm{Y}-\bm{Z}\right\|^{2}}{2\sqrt{1+\left\|\bm{Z}\right\|^{2}}}$
TLD divergence | $-\ln\det\bm{Y}$ |  $\delta_{\mathrm{TLD}}\left(\bm{Y},\bm{Z}\right)=\cfrac{\ln\det\left(\bm{Z}\bm{Y}^{-1}\right)+\operatorname{tr}\left(\bm{Z}^{-1}\bm{Y}\right)-N}{2\sqrt{1+\left\|\bm{Z}^{\mathsf{H}}\right\|^{2}}}$
TVN divergence | $\operatorname{tr}\left(\bm{Y}\operatorname{Log}\bm{Y}-\bm{Y}\right)$ |  $\delta_{\mathrm{TVN}}\left(\bm{Y},\bm{Z}\right)=\cfrac{\operatorname{tr}\left(\bm{Y}\left(\operatorname{Log}\bm{Y}-\operatorname{Log}\bm{Z}\right)-\bm{Y}+\bm{Z}\right)}{2\sqrt{1+\left\|\left(\operatorname{Log}\bm{Z}\right)^{\mathsf{H}}\right\|^{2}}}$

Before introducing the TBD medians, we briefly recall their mean counterparts.
Details can be found in [25]. The TBD mean of $m$ HPD matrices
$\\{\bm{R}_{i}\\}_{i=1}^{m}$ is defined by

$\displaystyle\underset{\bm{R}\in\mathscr{P}(N,\mathbb{C})}{\operatorname{arg}\min}\dfrac{1}{m}\sum_{i=1}^{m}\delta_{F}\left(\bm{R},\bm{R}_{i}\right).$
(3.12)

As the TBD $\delta_{F}(\bm{R},\bm{R}_{i})$ is defined in the convex HPD space
and is strictly convex with respect to $\bm{R}$, the TBD mean exists uniquely.
It solves the next matrix equation,

$\displaystyle\nabla
F(\bm{R})=\left({\sum\limits_{j=1}^{m}\frac{1}{\sqrt{1+\left\lVert\nabla
F(\bm{R}_{j})\right\rVert^{2}}}}\right)^{-1}\left({\sum_{i=1}^{m}\frac{\nabla
F(\bm{R}_{i})}{\sqrt{1+\left\lVert\nabla
F(\bm{R}_{i})\right\rVert^{2}}}}\right).$ (3.13)

In particular, the TBD means corresponding to the convex functions $F$ given
by Table 1 can be calculated analytically [25]. We summarize them up as
follows. The TSL mean, TLD mean, and TVN mean of $m$ HPD matrices
$\\{\bm{R}_{i}\\}_{i=1}^{m}$ are respectively given by

$\displaystyle\overline{\bm{R}}_{\mathrm{TSL}}$ $\displaystyle=$
$\displaystyle\left(\sum\limits_{j=1}^{m}\frac{1}{\sqrt{1+\left\lVert\bm{R}_{j}\right\rVert^{2}}}\right)^{-1}\left(\sum_{i=1}^{m}\frac{\bm{R}_{i}}{\sqrt{1+\left\lVert\bm{R}_{i}\right\rVert^{2}}}\right),$
(3.14) $\displaystyle\overline{\bm{R}}_{\mathrm{TLD}}$ $\displaystyle=$
$\displaystyle\left(\sum_{j=1}^{m}\frac{1}{\sqrt{1+\left\lVert\bm{R}^{-1}_{j}\right\rVert^{2}}}\right)\left(\sum_{i=1}^{m}\frac{\bm{R}_{i}^{-1}}{\sqrt{1+\left\lVert\bm{R}_{i}^{-1}\right\rVert^{2}}}\right)^{-1},$
(3.15) $\displaystyle\overline{\bm{R}}_{\mathrm{TVN}}$ $\displaystyle=$
$\displaystyle\exp\left(\left(\sum_{j=1}^{m}\frac{1}{\sqrt{1+\left\lVert\operatorname{Log}\bm{R}_{j}\right\rVert^{2}}}\right)^{-1}\sum_{i=1}^{m}\frac{\operatorname{Log}\bm{R}_{i}}{\sqrt{1+\left\lVert\operatorname{Log}\bm{R}_{j}\right\rVert^{2}}}\right).$
(3.16)

######  Definition 2.

Assume $F$ is a differentiable and strictly convex function. Let $\delta_{F}$
be the associated TBD. The TBD median of $m$ HPD matrices
$\\{\bm{R}_{i}\\}_{i=1}^{m}$ is defined by

$\underset{\bm{R}\in\mathscr{P}(N,\mathbb{C})}{\operatorname{arg}\min}\dfrac{1}{m}\sum_{i=1}^{m}\left\\{\delta_{F}\left(\bm{R},\bm{R}_{i}\right)\right\\}^{\frac{1}{2}}.$
(3.17)

######  Theorem 3.

If the TBD median (3.17) exists, then it satisfies the matrix equation,

$\nabla
F\left(\bm{R}\right)=\left(\sum_{j=1}^{m}\frac{1}{\left\\{\delta_{F}\left(\bm{R},\bm{R}_{j}\right)\right\\}^{\frac{1}{2}}\sqrt{1+\left\|\nabla
F(\bm{R}_{j})\right\|^{2}}}\right)^{-1}\left(\sum_{i=1}^{m}\frac{\nabla
F\left(\bm{R}_{i}\right)}{\left\\{\delta_{F}\left(\bm{R},\bm{R}_{i}\right)\right\\}^{\frac{1}{2}}\sqrt{1+\left\|\nabla
F(\bm{R}_{i})\right\|^{2}}}\right).$ (3.18)

###### Proof.

Let $G(\bm{R})$ be the function to be optimized, namely

$G(\bm{R}):=\frac{1}{m}\sum\limits_{i=1}^{m}\left\\{\delta_{F}\left(\bm{R},\bm{R}_{i}\right)\right\\}^{\frac{1}{2}}.$
(3.19)

The proof is completed by setting $\nabla G(\bm{R})=0$ where the gradient of
$G(\bm{R})$ is

$\nabla G(\bm{R})=\frac{1}{2m}\sum_{i=1}^{m}\frac{\nabla F(\bm{R})-\nabla
F(\bm{R}_{i})}{\left\\{\delta_{F}\left(\bm{R},\bm{R}_{i}\right)\right\\}^{\frac{1}{2}}\sqrt{1+\left\lVert\nabla
F(\bm{R}_{i})\right\rVert^{2}}}.$ (3.20)

∎

Unfortunately, it is difficult to solve the matrix equation (3.18)
analytically, and alternatively we seek for its numerical solutions using the
fixed-point algorithm [32].

######  Proposition 4.

The TSL median, the TLD median, and the TVN median of $m$ HPD matrices
$\\{\bm{R}_{i}\\}_{i=1}^{m}$ can respectively be calculated by the following
fixed-point algorithms in Table 2.

Table 2: The fixed-point algorithms for computing TBD medians TBD | The fixed-point algorithm
---|---
TSL | $\overline{\bm{R}}_{t+1}=\left(\sum\limits_{j=1}^{m}\cfrac{1}{\left\\{\delta_{\mathrm{TSL}}\left(\overline{\bm{R}}_{t},\bm{R}_{j}\right)\right\\}^{\frac{1}{2}}\sqrt{1+\left\lVert\bm{R}_{j}\right\rVert^{2}}}\right)^{-1}\left(\sum\limits_{i=1}^{m}\cfrac{\bm{R}_{i}}{\left\\{\delta_{\mathrm{TSL}}\left(\overline{\bm{R}}_{t},\bm{R}_{i}\right)\right\\}^{\frac{1}{2}}\sqrt{1+\left\lVert\bm{R}_{i}\right\rVert^{2}}}\right)$
TLD | $\overline{\bm{R}}_{t+1}=\left(\sum\limits_{j=1}^{m}\cfrac{1}{\left\\{\delta_{\mathrm{TLD}}\left(\overline{\bm{R}}_{t},\bm{R}_{j}\right)\right\\}^{\frac{1}{2}}\sqrt{1+\left\lVert\bm{R}_{j}^{-1}\right\rVert^{2}}}\right)\left({\sum\limits_{i=1}^{m}\cfrac{\bm{R}_{i}^{-1}}{\left\\{\delta_{\mathrm{TLD}}\left(\overline{\bm{R}}_{t},\bm{R}_{i}\right)\right\\}^{\frac{1}{2}}\sqrt{1+\left\lVert\bm{R}_{i}^{-1}\right\rVert^{2}}}}\right)^{-1}$
TVN | $\overline{\bm{R}}_{t+1}=\mathrm{exp}\left(\left(\sum\limits_{j=1}^{m}\cfrac{1}{\left\\{\delta_{\mathrm{TVN}}\left(\overline{\bm{R}}_{t},\bm{R}_{j}\right)\right\\}^{\frac{1}{2}}\sqrt{1+\left\lVert\left(\operatorname{Log}\bm{R}_{j}\right)^{\mathsf{H}}\right\rVert^{2}}}\right)^{-1}\sum\limits_{i=1}^{m}\cfrac{\operatorname{Log}\bm{R}_{i}}{\left\\{\delta_{\mathrm{TVN}}\left(\overline{\bm{R}}_{t},\bm{R}_{i}\right)\right\\}^{\frac{1}{2}}\sqrt{1+\left\lVert\left(\operatorname{Log}\bm{R}_{i}\right)^{\mathsf{H}}\right\rVert^{2}}}\right)$

In the matrix-CFAR, the matrix $\bm{R}_{g}$ in (2.11) will be replaced by the
TBD means, the TBD medians, the RD mean, and so forth; see also Fig. 1.

Numerical simulations in the next sections show that the Riemannian gradient
descent algorithms introduced in Proposition 1 and the fixed-point algorithms
in Table 2 are all convergent with the initial value chosen as the arithmetic
mean of a set of HPD matrices; the tolerance of each iteration is $10^{-3}$.

##  4 Numerical performance analysis

To confirm and compare the performance of the matrix-CFAR via TBD means and
medians, the following numerical simulations are conducted. For the
comparison, we also show results of the RD mean induced from the AIRM, the
GLRT, and the ANMF. Two types of clutter are considered. One is the Gaussian
clutter modeled by an $N$ dimensional random complex vector possessing the
complex circular Gaussian distribution $CN\left(0,\bm{\Sigma}\right)$ of zero-
mean and a covariance matrix $\bm{\Sigma}$. Its probability density function
is given by

$p\left(\bm{z}\ |\
0,\bm{\Sigma}\right)=\frac{1}{\pi^{N}\det\left(\bm{\Sigma}\right)}\exp\left(-\bm{z}^{\mathsf{H}}\bm{\Sigma}^{-1}\bm{z}\right).$
(4.1)

The other is the compound-Gaussian clutter which is a model commonly used for
characterising heavy-tailed clutter distributions, for instance sea clutter
[40], defined by

$\bm{c}=\sqrt{\tau}\bm{z}.$ (4.2)

Here $\bm{z}$ is called the speckle component which is fast fluctuating and
modeled by a complex Gaussian distribution. The positive real random variable
$\tau$ is called the texture component which is comparatively slowly varying.
In the simulation, the texture component $\tau$ is assumed to subject to a
gamma distribution with a shape parameter $\alpha$ and a scale parameter
$\beta$ whose probability density function is

$q\left(\tau\ |\
\alpha,\beta\right)=\frac{1}{\beta^{\alpha}\Gamma\left(\alpha\right)}\tau^{\alpha-1}\exp\left(\frac{-\tau}{\beta}\right),$
(4.3)

where $\Gamma(\cdot)$ denotes the gamma function. Consequently, the
distribution of clutter is defined by the complex K-distribution which is a
special compound-Gaussian model [40]. In the simulations, we generate the
Gaussian clutter, as well as the speckle component in the compound-Gaussian
clutter, from $CN\left(0,\bm{\Sigma}\right)$ of which the known covariance
matrix $\bm{\Sigma}$ is given by

$\bm{\Sigma}=\bm{\Sigma}_{0}+\bm{I}_{N},$ (4.4)

where

$\bm{\Sigma}_{0}\left(i,j\right)=\sigma_{c}^{2}\rho^{\left|i-j\right|}\exp\left(\mathrm{i}2\pi
f_{c}(i-j)\right),\ \ i,j=1,2,\ldots,N.$

Here, $\bm{I}_{N}$ is the $N$ dimensional identity matrix. In the simulation,
the clutter-to-noise ratio $\sigma_{c}$ is set to
$\sigma_{c}^{2}=20\mathrm{dB}$, the one-lag coefficient is $\rho=0.9$, and the
normalized Doppler frequency of clutter is $f_{c}=0.2$. Furthermore, the shape
parameter is set to $\alpha=4$ and the scale parameter is set to $\beta=3$.
The signal to clutter ratio (SCR) is defined by

$\mathrm{SCR}=|\xi|^{2}\bm{s}^{\mathsf{H}}\bm{R}^{-1}\bm{s},$ (4.5)

where $\xi$ is the amplitude coefficient, $\bm{s}$ is the target vector
defined as (2.3) and $\bm{R}$ is the clutter covariance matrix. Furthermore,
the false alarm rate $P_{fa}$ is $10^{-3}$, and the signal normalized Doppler
frequency $f_{d}$ is $0.2$. Two interferences are inserted into clutter data
and its normalized Doppler frequency $f_{i}$ is also $0.2$. Dimension of the
signal is set to $N=8$, and we study the cases when the numbers of
observations are $m=8,12,16$, respectively.

The threshold $\gamma$ is derived by $100/P_{fa}$, and the probability of
detection $P_{d}$ is derived by 2000 independent trails. Fig. 2 shows that in
both types of clutter, the detection performance improves as the number of
observations increases. The GLRT and the ANMF are not valid when the number of
observations is small but the matrix-CFAR detectors work well. In the case of
Gaussian clutter (see Figs. 2, 2, 2), the TBD means and medians are
advantageous over the RD mean and median, the GLRT and the ANMF and in
particular, the TLD median has the best performance. In the compound-Gaussian
clutter (see Figs. 2, 2, 2), the TVN median has the best performance when the
number of observations is not large enough. In the case of $m=16$, the ANMF is
the best among all detectors. Among the TBD means and medians, the TLD and TVN
medians behave better than the corresponding means, but surprisingly the TSL
median is worse than the TSL mean.

######  Remark 5.

Generally speaking, the simulations show that the matrix-CFAR performs better
than the GLRT and ANMF particularly when limited data is available. This is
probably due to that the effect of clutter and noise in the matrix-CFAR is
less than that in the GLRT and ANMF, since the threshold of the GLRT and ANMF
is derived by using contaminated sample data directly, while the matrix-CFAR
leverages the covariance matrix of sample data.

(a) Gaussian clutter, $m=8$

(b) Compound-Gaussian clutter, $m=8$

(c) Gaussian clutter, $m=12$

(d) Compound-Gaussian clutter, $m=12$

(e) Gaussian clutter, $m=16$

(f) Compound-Gaussian clutter, $m=16$

Figure 2: $P_{d}$ versus signal-to-clutter ratio for the TBD means and
medians, the RD mean and median, the GLRT and the ANMF.

##  5 Robustness analysis

In real detection, the received signal may also include outliers, and hence
robustness about estimate of the covariance matrix is crucial. To evaluate the
robustness of the detectors to outliers, the influence function has usually
been employed [25]. It was noted that the influence functions with respect to
the TSL mean and the TLD mean can easily be calculated analytically but not
for the TVN mean [25]. The situation for the median counterparts is often more
complicated. In this section, we will propose an orthogonal basis for
Hermitian matrices; by doing so, we are allowed to rewrite the equations of
influence functions as a linear system, that can be solved analytically in a
routine manner.

###  5.1 An orthogonal basis for Hermitian matrices

As we will see later that influences functions are defined as a small
perturbation of the means or medians that depends on a Hermitian matrix. A
well-known basis for the $N^{2}$ dimensional space of Hermitian matrices
$\mathscr{H}(N,\mathbb{C})$ is given by $N^{2}$ number of matrices:

$\displaystyle\bm{E}_{ii},$ $\displaystyle 1\leq i\leq N;$ (5.1)
$\displaystyle\bm{E}_{ij}+\bm{E}_{ji},~{}$ $\displaystyle 1\leq i<j\leq N;$
$\displaystyle\mathrm{i}\left(\bm{E}_{ij}-\bm{E}_{ji}\right),~{}$
$\displaystyle 1\leq i<j\leq N,$

where $\bm{E}_{ij}$ represents the $N\times N$ matrix with 1 at row $i$ and
column $j$ and 0 elsewhere. To ease the application in the following
robustness analysis, we normalize them into the following orthogonal basis

$\displaystyle\bm{E}_{ii},$ $\displaystyle 1\leq i\leq N;$ (5.2)
$\displaystyle\frac{1}{\sqrt{2}}\left(\bm{E}_{ij}+\bm{E}_{ji}\right),~{}$
$\displaystyle 1\leq i<j\leq N;$
$\displaystyle\frac{\mathrm{i}}{\sqrt{2}}\left(\bm{E}_{ij}-\bm{E}_{ji}\right),~{}$
$\displaystyle 1\leq i<j\leq N,$

and denote it by $\\{\mathsf{E}_{k}\\}_{k=1}^{N^{2}}$. Each Hermitian matrix
$\bm{H}\in\mathscr{H}(N,\mathbb{C})$ can be represented as a linear
combination using the orthogonal basis, namely

$\bm{H}=h^{k}\mathsf{E}_{k},\quad h^{k}\in\mathbb{R},$ (5.3)

where the summation is over $\\{1,\ldots,N^{2}\\}$ and the Einstein summation
convention is applied here and all through the paper. It is immediate to check
that the basis $\\{\mathsf{E}_{k}\\}_{k=1}^{N^{2}}$ is orthogonal by showing
that

$\langle\mathsf{E}_{k},\mathsf{E}_{l}\rangle=\delta_{kl},\quad\forall
k,l=1,\ldots,N^{2},$ (5.4)

where $\delta_{kl}$ is the Kronecker delta that gives 1 when $k=l$ and 0
otherwise. Given a Hermitian matrix $\bm{H}$, it is immediate to write down
its coordinate expression as

$\bm{H}=h^{k}\mathsf{E}_{k},\text{ with
}h^{k}=\langle\bm{H},\mathsf{E}_{k}\rangle.$ (5.5)

###  5.2 Robustness via influence functions

Let $\overline{\bm{R}}$ be the TBD medians (or TBD means, RD mean, RD median)
of a given set of $m$ HPD matrices $\\{\bm{R}_{i}\\}_{i=1}^{m}$, and let
$\widehat{\bm{R}}$ be the TBD medians (or TBD means, RD mean, RD median) of
the mixed HPD data including $n$ number of outliers
$\\{\bm{P}_{j}\\}_{j=1}^{n}$. By thinking of the outliers as a perturbation
$\varepsilon\left(\varepsilon<<1\right)$ of each mean or median,
$\widehat{\bm{R}}$ and $\overline{\bm{R}}$ can be related by

$\widehat{\bm{R}}=\overline{\bm{R}}+\varepsilon\bm{H}\left(\\{\bm{R}_{i}\\}_{i=1}^{m},\\{\bm{P}_{j}\\}_{j=1}^{n}\right)+\mathcal{O}\left(\varepsilon^{2}\right),$
(5.6)

where
$\bm{H}\left(\\{\bm{R}_{i}\\}_{i=1}^{m},\\{\bm{P}_{j}\\}_{j=1}^{n}\right)$ is
a Hermitian matrix depending on the outliers, and the $m$ HPD matrices
$\\{\bm{R}_{i}\\}_{i=1}^{m}$, the latter of which are assumed to be given and
fixed. The function

$f\left(\\{\bm{R}_{i}\\}_{i=1}^{m},\\{\bm{P}_{j}\\}_{j=1}^{n}\right)\coloneqq\dfrac{\left\lVert\bm{H}\left(\\{\bm{R}_{i}\\}_{i=1}^{m},\\{\bm{P}_{j}\\}_{j=1}^{n}\right)\right\rVert}{\left\lVert\overline{\bm{R}}\right\rVert}$
(5.7)

is defined as the (normalized) influence function.

To calculate the corresponding influence functions, define the objective
function $G(\bm{R})$ for the contaminated $m+n$ number of HPD matrices as
follows

$G(\bm{R}):=(1-\varepsilon)\frac{1}{m}\sum_{i=1}^{m}d\left(\bm{R},\bm{R}_{i}\right)+\varepsilon\frac{1}{n}\sum_{j=1}^{n}d\left(\bm{R},\bm{P}_{j}\right),$
(5.8)

where $d$ is either the distance or divergence functions for medians, and
square of them for means. Since $\bm{\widehat{R}}$ is any mean or median of
the mixed $m+n$ HPD matrices, we have $\nabla
G(\bm{\widehat{R}})=\bm{0},$111Note that in the Riemannian case, this should
be $\operatorname{grad}G(\widehat{\bm{R}})=0$, which is nevertheless
equivalent to $\nabla G(\widehat{\bm{R}})=0$ for HPD manifolds equipped with
the AIRM. that is,

$(1-\varepsilon)\frac{1}{m}\sum_{i=1}^{m}\nabla
d\left(\bm{\widehat{R}},\bm{R}_{i}\right)+\varepsilon\frac{1}{n}\sum_{j=1}^{n}\nabla
d\left(\bm{\widehat{R}},\bm{P}_{j}\right)=\bm{0}.$ (5.9)

Differentiating $\nabla G(\widehat{\bm{R}})=\bm{0}$ about$\varepsilon$ at
$\varepsilon=0$ , we obtain a matrix equation for $\bm{H}$:

$\frac{1}{m}\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\left(\sum_{i=1}^{m}\nabla
d\left(\bm{\widehat{R}},\bm{R}_{i}\right)\right)+\frac{1}{n}\sum_{j=1}^{n}\nabla
d\left(\overline{\bm{R}},\bm{P}_{j}\right)=\bm{0}.$ (5.10)

Next, we are going to show that the above equation can be solved analytically
using the orthogonal basis (5.2) defined above.

Recall that the gradient of a function defined in HPD spaces is a Hermitian
matrix. Using the orthogonal basis $\\{\mathsf{E}_{k}\\}_{k=1}^{N^{2}}$, we
can write

$\sum_{i=1}^{m}\nabla
d\left(\bm{\widehat{R}},\bm{R}_{i}\right)=\xi^{l}(\widehat{\bm{R}})\mathsf{E}_{l}$
(5.11)

and

$\sum_{j=1}^{m}\nabla
d\left(\overline{\bm{R}},\bm{P}_{j}\right)=\phi^{l}\mathsf{E}_{l}.$ (5.12)

Letting $\bm{H}=h^{s}\mathsf{E}_{s}$, the first term in (5.10) then becomes

$\displaystyle\dfrac{1}{m}\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\xi^{l}(\widehat{\bm{R}})\mathsf{E}_{l}$
$\displaystyle=$
$\displaystyle\dfrac{1}{m}\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\xi^{l}(\overline{\bm{R}}+\varepsilon\bm{H})\mathsf{E}_{l}$
(5.13) $\displaystyle=$
$\displaystyle\dfrac{1}{m}\langle\nabla\xi^{l}(\overline{\bm{R}}),\bm{H}\rangle\mathsf{E}_{l}$
$\displaystyle=$
$\displaystyle\dfrac{h^{k}}{m}\langle\nabla\xi^{l}(\overline{\bm{R}}),\mathsf{E}_{k}\rangle\mathsf{E}_{l}.$

Definition of the gradient (3.11) (induced to HPD spaces) is applied here.
Although a linear system about $\\{h^{k}\\}_{k=1}^{N^{2}}$ can be obtained,
unfortunately, as $\\{\xi^{l}(\overline{\bm{R}})\\}_{l=1}^{N^{2}}$ is the
coordinates whose gradient can be practically difficult to be calculated, in
particular, in numerical simulations. Alternatively, noticing the first term
in (5.10) is Hermitian and linear about $\\{h^{k}\\}_{k=1}^{N^{2}}$ (from
(5.13)), we can simply write

$\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\left(\sum_{i=1}^{m}\nabla
d\left(\bm{\widehat{R}},\bm{R}_{i}\right)\right)=h^{k}\theta_{k}^{l}\mathsf{E}_{l}.$
(5.14)

Consequently, we have successfully written the matrix equation (5.10) as a
linear system about $\\{h^{k}\\}_{k=1}^{N^{2}}$ with matrix coefficients:

$\frac{1}{m}h^{k}\theta_{k}^{l}\mathsf{E}_{l}+\frac{1}{n}\phi^{l}\mathsf{E}_{l}=\bm{0}.$
(5.15)

Conducting Frobenius metric on both sides with each $\mathsf{E}_{s}$ and
making use of the orthogonality, we obtain a linear system for
$\\{h^{k}\\}_{k=1}^{N^{2}}$, reading

$\frac{1}{m}\theta_{k}^{s}h^{k}+\frac{1}{n}\phi^{s}=0,\quad s=1,\ldots,N^{2},$
(5.16)

which can be immediately solved if the coefficient matrix $(\theta_{k}^{s})$
is invertible. We summarize the above deduction into the following theorem.

######  Theorem 6.

Under the orthogonal basis $\\{\mathsf{E}_{k}\\}_{k=1}^{N^{2}}$ defined by
(5.2), the Hermitian matrix $\bm{H}=h^{k}\mathsf{E}_{k}$ in the influence
function (5.7) is given by linear system (5.16).

Influence function of the RD mean. Taking the RD mean as an illustrative
example, we show how the above deduction helps us to derive the matrix
$\bm{H}$ and hence the corresponding influence function. Regarding the
contaminated HPD matrices, define $G(\bm{R})$ as the objective function, i.e.,

$G(\bm{R}):=(1-\varepsilon)\frac{1}{m}\sum_{i=1}^{m}d^{2}_{R}\left(\bm{R},\bm{R}_{i}\right)+\varepsilon\frac{1}{n}\sum_{j=1}^{n}d^{2}_{R}\left(\bm{R},\bm{P}_{j}\right),$
(5.17)

where

$d^{2}_{R}\left(\bm{R},\bm{R}_{i}\right)=\left\lVert\operatorname{Log}\left(\bm{R}_{i}^{-1}\bm{R}\right)\right\rVert^{2}=\operatorname{tr}\left(\operatorname{Log}^{2}\left(\bm{R}_{i}^{-1}\bm{R}\right)\right).\\\
$ (5.18)

Differentiate $\operatorname{grad}G(\widehat{\bm{R}})=\bm{0}$ with respect to
$\varepsilon$ at $\varepsilon=0$ and obtain

$\frac{1}{m}\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\sum_{i=1}^{m}\bm{\widehat{R}}\operatorname{Log}\left(\bm{R}_{i}^{-1}\bm{\widehat{R}}\right)+\frac{1}{n}\sum_{j=1}^{n}\overline{\bm{R}}\operatorname{Log}\left(\bm{P}_{j}^{-1}\overline{\bm{R}}\right)=\bm{0}.$
(5.19)

The second term can be simply written as

$\frac{1}{n}\sum_{j=1}^{n}\overline{\bm{R}}\operatorname{Log}\left(\bm{P}_{j}^{-1}\overline{\bm{R}}\right)=\frac{1}{n}\phi^{k}\mathsf{E}_{k},$
(5.20)

while the first term can be expanded as [25, 31]

$\frac{1}{m}\sum_{i=1}^{m}\left(\bm{{H}}\operatorname{Log}\left(\bm{R}_{i}^{-1}\overline{\bm{R}}\right)+\overline{\bm{R}}\int_{0}^{1}\left[(\bm{R}_{i}^{-1}\overline{\bm{R}}-\bm{I})\tau+\bm{I}\right]^{-1}\bm{R}_{i}^{-1}\bm{H}\left[(\bm{R}_{i}^{-1}\overline{\bm{R}}-\bm{I})\tau+\bm{I}\right]^{-1}d\tau\right).$

Substituting $\bm{H}=h^{k}\mathsf{E}_{k}$ to the first term, it can be written
in the form $\frac{1}{m}h^{k}\theta_{k}^{s}\mathsf{E}_{s}$, where
$\theta_{k}^{s}$ is given by

$\theta_{k}^{s}=\left\langle\sum_{i=1}^{m}\left(\mathsf{E}_{k}\operatorname{Log}\left(\bm{R}_{i}^{-1}\overline{\bm{R}}\right)+\overline{\bm{R}}\int_{0}^{1}\left[(\bm{R}_{i}^{-1}\overline{\bm{R}}-\bm{I})\tau+\bm{I}\right]^{-1}\bm{R}_{i}^{-1}\mathsf{E}_{k}\left[(\bm{R}_{i}^{-1}\overline{\bm{R}}-\bm{I})\tau+\bm{I}\right]^{-1}d\tau\right),\mathsf{E}_{s}\right\rangle.$

Finally, we obtain the linear system (5.16) about $\\{h^{k}\\}_{k=1}^{N^{2}}$
by conducting the Frobenius metric with respect to each element of the
orthogonal basis, which will give us the influence function for the RD mean.
Those of the RD median, as well as the TBD means and medians can also be
analytically computed in a similar way; details are omitted here.

In the following simulations, we firstly generate $m$ data of sample, each of
which is $N$ dimensional and complex valued, generated from
$CN\left(0,\bm{\Sigma}\right)$ where the covariance matrix $\bm{\Sigma}$ is
given by (4.4). The means and medians $\overline{\bm{R}}$ are calculated from
the $m$ HPD covariance matrices $\\{\bm{R}_{i}\\}_{i=1}^{m}$ of the $m$ number
of sample data. We then mix $n$ outliers $\\{\bm{P}_{j}\\}_{j=1}^{n}$ with the
sample data; the outliers are modeled as the covariance matrix of
$\xi\bm{s}+\bm{c}$ where $\bm{s}$ denotes the steering vector, $\bm{c}$ is the
Gaussian clutter with the covariance matrix $\bm{\Sigma}$, and $\xi$ is the
amplitude coefficient derived by the signal-to-clutter ratio which is set to
40 dB as (4.5). In the numerical simulations, the number of sample is $m=50$
while the number of outliers $n$ varies from 1 to 40. Regarding the influence
function of the SCR, the $\overline{\bm{R}}$ is calculated from the $m$ number
of sample data and the $\widehat{\bm{R}}$ is calculated from the $m+n$ number
contaminated data by using (2.4). Therefore the influence function of the SCR
is derived similarly to that in (5.7). Fig. 3 plots the influence functions
which are calculated by averaging 1000 trails and are shown in the logarithmic
scale. It is noticed that the influence functions of the TBD medians are
smaller than the TBD means and the geometric means and medians except for the
TSL mean are more robust than the SCM. However, regardless of its worse
detection performance compared with the TBD means and medians, the RD mean and
median are more robust, although the RD mean and median take much longer time
than the TBD means and medians. Therefore, from a comprehensive perspective
the TBD medians can serve as effective estimators.

Figure 3: Robustness analysis of the SCM, the RD mean and median, and the TBD
means and medians via the accurate influence functions.

##  6 Conclusions

The TBD medians were defined and applied to matrix-CFAR detectors
corresponding to the well-known TBDs, i.e., the TSL, and the TLD and TVN
divergences. To illustrate the improvement of detection performance, they were
compared with their mean counterparts, the RD mean and median associated to
the AIRM, as well as the conventional GLRT detector and the ANMF. Simulations
showed that the TLD and TVN medians performed better than all means regardless
of the number of observations. To analyze the robustness of those detectors to
outliers, we transformed the intractable governing matrix equation into a
linear system by introducing an orthogonal basis for Hermitian matrices. This
allowed us to calculate the corresponding influence functions analytically and
made the robustness analysis accurate and precise. We believe this methodology
works equally well for other geometric detectors in matrix-CFAR. By
simulations, the RD mean is most robust among all of them, while the TBD
medians are more robust compared with the TBD means. From an overall point of
view, the TBD medians, in particular the TVN median, have best performance and
satisfactory robustness.

As the optimization problem is convex and defined in a convex space, the TBD
means exist uniquely. Unfortunately, an analytical verification of the
existence and uniqueness of TBD medians is yet to be constructed, although
numerical simulations suggest so. Future research in applications includes
detection problems for real-world data using the TBD medians. Finally, it
would also be interesting to extend the current study to other covariance
estimators, such as the diagonal loading [15, 1], and structured covariance
interference (e.g., [12]).

## Acknowledgements.

The authors thank Prof. Sergey Yurish for the kind invitation and for the
hospitality during the ASPAI’2021 Conference, in which part of the work was
presented [34]. The authors would also like to thank Xiaoqiang Hua and Hiroki
Sato for helpful discussions and to the referees for the constructive
comments. This work was partially supported by Japan Society for the Promotion
of Science KAKENHI (No. JP20K14365), Japan Science and Technology Agency CREST
(No. JPMJCR1914), and the Fukuzawa Fund of Keio University.

## References

  * [1] Y. I. Abramovich, N. K. Spencer, and A. Y. Gorokhov.  Modified GLRT and AMF framework for adaptive detectors.  IEEE Transactions on Aerospace and Electronic Systems, 43(3):1017–1051, 2007.
  * [2] M. Arnaudon, F. Barbaresco, and L. Yang.  Medians and means in Riemannian geometry: Existence, uniqueness and computation.  In R. Bhatia and F. Nielsen, editors, Matrix Information Geometry, pages 169–198, 2012.
  * [3] M. Arnaudon, F. Barbaresco, and L. Yang.  Riemannian medians and means with applications to radar signal processing.  IEEE Journal of Selected Topics in Signal Processing, 7(4):595–604, 2013.
  * [4] F. Bandiera, O. Besson, and G. Ricci.  Knowledge-aided covariance matrix estimation and adaptive detection in compound-Gaussian noise.  IEEE Transactions on Signal Processing, 58(10):5391–5396, Oct. 2010\.
  * [5] F. Barbaresco.  Innovative tools for radar signal processing based on Cartan’s geometry of SPD matrices $\&$ information geometry.  In 2008 IEEE International Radar Conference, pages 1–6, 2008.
  * [6] F. Barbaresco.  Interactions between symmetric cone and information geometries: Bruhat-tits and siegel spaces models for high resolution autoregressive doppler imagery.  In F. Nielsen, editor, Emerging Trends In Visual Computing, pages 124–163, Palaiseau, France, 2009.
  * [7] F. Barbaresco and U. Meier.  Radar monitoring of a wake vortex: Electromagnetic reflection of wake turbulence in clear air.  Comptes Rendus Physique, 11(1):54–67, 2010.
  * [8] O. Besson, J. Tourneret, and S. Bidon.  Knowledge-aided Bayesian detection in heterogeneous environments.  IEEE Signal Processing Letters, 14(5):355–358, May 2007.
  * [9] Y. Cabanes, F. Barbaresco, M. Arnaudon, and J. Bigot.  Toeplitz Hermitian positive definite matrix machine learning based on Fisher metric.  In Geometric Science of Information: 4th International Conference, GSI 2019, Toulouse, France, August 27-29, 2019, Proceedings, Geometric Science of Information, pages 261–270. Springer International Publishing, 2019.
  * [10] H. Chahrour, R. M. Dansereau, S. Rajan, and B. Balaji.  Target detection through Riemannian geometric approach with application to drone detection.  IEEE Access, 9:123950–123963, 2021.
  * [11] M. Charfi, Z. Chebbi, M. Moakher, and B. C. Vemuri.  Bhattacharyya median of symmetric positive-definite matrices and application to the denoising of diffusion-tensor fields.  In 2013 IEEE 10th International Symposium on Biomedical Imaging, pages 1227–1230, 2013.
  * [12] D. Ciuonzo, D. Orlando, and L. Pallotta.  On the maximal invariant statistic for adaptive radar detection in partially homogeneous disturbance with persymmetric covariance.  IEEE Signal Processing Letters, 23(12):1830–1834, Dec 2016.
  * [13] E. Conte, M. Lops, and G. Ricci.  Asymptotically optimum radar detection in compound-Gaussian clutter.  IEEE Transactions on Aerospace and Electronic Systems, 31(2):617–625, 1995.
  * [14] A. De Maio, A. Farina, and G. Foglia.  Knowledge-aided Bayesian radar detectors & their application to live data.  IEEE Transactions on Aerospace and Electronic Systems, 46(1):170–183, Jan. 2010.
  * [15] A. De Maio, L. Pallotta, J. Li, and P. Stoica.  Loading factor estimation under affine constraints on the covariance eigenvalues with application to radar target detection.  IEEE Transactions on Aerospace and Electronic Systems, 55(3):1269–1283, 2019.
  * [16] A. Decurninge and F. Barbaresco.  Robust Burg estimation of radar scatter matrix for autoregressive structured SIRV based on Fréchet medians.  IET Radar, Sonar & Navigation, 11(1):78–89, 2017.
  * [17] I. S. Dhillon and J. A. Tropp.  Matrix nearness problems with Bregman divergences.  SIAM Journal on Matrix Analysis and Applications, 29(4):1120–1146, 2008.
  * [18] X. Du, A. Aubry, A. De Maio, and G. Cui.  Toeplitz structured covariance matrix estimation for radar applications.  IEEE Signal Processing Letters, 27:595–599, 2020.
  * [19] P. T. Fletcher, S. Venkatasubramanian, and S. Joshi.  Robust statistics on Riemannian manifolds via the geometric median.  In 2008 IEEE Conference on Computer Vision and Pattern Recognition, pages 1–8. IEEE Computer Society, 2008.
  * [20] N. R. Goodman.  Statistical analysis based on a certain multivariate complex Gaussian distribution.  The Annals of Mathematical Statistics, 34(1):152–177, 1963.
  * [21] N. J. Higham.  Functions of Matrices: Theory and Computation.  SIAM, Philadelphia, 2008.
  * [22] X. Hua, Y. Cheng, H. Wang, Y. Qin, and Y. Li.  Geometric means and medians with applications to target detection.  IET Signal Processing, 11(6):711–720, 2017.
  * [23] X. Hua, Y. Cheng, H. Wang, Y. Qin, Y. Li, and W. Zhang.  Matrix CFAR detectors based on symmetrized Kullback–Leibler and total Kullback–Leibler divergences.  Digital Signal Processing, 69:106–116, 2017.
  * [24] X. Hua, H. Fan, Y. Cheng, H. Wang, and Y. Qin.  Information geometry for radar target detection with total Jensen–Bregman divergence.  Entropy, 20(4):256, 2018.
  * [25] X. Hua, Y. Ono, L. Peng, Y. Cheng, and H. Wang.  Target detection within nonhomogeneous clutter via total Bregman divergence-based matrix information geometry detectors.  IEEE Transactions on Signal Processing, 69:4326–4340, 2021.
  * [26] X. Hua, Y. Ono, L. Peng, and Y. Xu.  Unsupervised learning discriminative MIG detectors in nonhomogeneous clutter.  IEEE Transactions on Communications, 70:4107–4120, 2022.
  * [27] E. J. Kelly.  Adaptive detection in non-stationary interference, part I and part II.  Technical report, Lincoln Laboratory, MIT, June 25, 1985.
  * [28] J. Lapuyade-Lahorgue and F. Barbaresco.  Radar detection using Siegel distance between autoregressive processes, application to HF and X-band radar.  In 2008 IEEE Radar Conference, pages 1–6, May 2008.
  * [29] N. Li, H. Yang, G. Cui, L. Kong, and Q. H. Liu.  Adaptive two-step Bayesian MIMO detectors in compound-Gaussian clutter.  Signal Processing, 161:1–13, 2019.
  * [30] M. Martorella, S. Gelli, and A. Bacci.  Ground moving target imaging via SDAP-ISAR processing: Review and new trends.  Sensors, 21(7):2391, 2021.
  * [31] M. Moakher.  A differential geometric approach to the geometric mean of symmetric positive-definite matrices.  SIAM Journal on Matrix Analysis and Applications, 26(3):735–747, 2005.
  * [32] M. Moakher.  On the averaging of symmetric positive-definite tensors.  Journal of Elasticity, 82(3):273–296, 2006.
  * [33] M. Moakher and M. Zéraï.  The Riemannian geometry of the space of positive-definite matrices and its application to the regularization of positive-definite matrix-valued data.  Journal of Mathematical Imaging and Vision, 40(2):171–187, 2011\.
  * [34] Y. Ono, L. Peng, and H. Sato.  The matrix-CFAR via total Bregman divergence medians in signal detection.  In S. Yurish, editor, 3rd International Conference on Advances in Signal Processing and Artificial Intelligence, pages 18–22. International Frequency Sensor Association Publishing, 2021.
  * [35] I. S. Reed, J. D. Mallett, and L. E. Brennan.  Rapid convergence rate in adaptive arrays.  IEEE Transactions on Aerospace and Electronic Systems, AES-10(6):853–863, 1974.
  * [36] M. Richards.  Fundamentals of Radar Signal Processing.  McGraw-Hill Professional Professional, 2nd edition, 2014.
  * [37] S. T. Smith.  Geometric Optimization Methods for Adaptive Filtering.  PhD Thesis, Harvard University, Cambridge, Massachusetts, 1993.
  * [38] C. Udrişte.  Convex Functions and Optimization Methods on Riemannian Manifolds.  Springer Science+Business Media, B.V., Dordrecht, 1994.
  * [39] P. Wang, Z. Wang, H. Li, and B. Himed.  Knowledge-aided parametric adaptive matched filter with automatic combining for covariance estimation.  IEEE Transactions on Signal Processing, 62(18):4713–4722, 2014\.
  * [40] K. Ward, R. Tough, and S. Watts.  Sea Clutter: Scattering, the K Distribution and Radar Performance.  The institution of Engineering and Technology, London, 2013.
  * [41] Z. Yang, X. Li, H. Wang, and R. Fa.  Knowledge-aided STAP with sparse-recovery by exploiting spatio-temporal sparsity.  IET Signal Processing, 10(2):150–161, 2016.
  * [42] L. Ye, Q. Yang, and W. Deng.  Matrix constant false alarm rate (MCFAR) detector based on information geometry.  In 2016 IEEE Information Technology, Networking, Electronic and Automation Control Conference, pages 211–215, 2016.
  * [43] M. A. Zatman.  ABF limitations when using either Toeplitz covariance matrix estimators or the parametric vector AR technique.  In Conference Record of Thirty-Fifth Asilomar Conference on Signals, Systems and Computers, volume 2, pages 1111–1115, 2001.