{'fraction_non_alphanumeric': 0.03725868725868726, 'fraction_numerical': 0.011583011583011582, 'mean_word_length': 4.480533220482438, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Towards a median signal detector through the total Bregman divergence and its robustness analysis

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.


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

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 x is only the clutter c and the alternative hypothesis H 0 represents that x contains both the target signal s and the clutter c, namely where ξ denotes the unknown amplitude coefficient which is complex scalar-valued. The observation data x and the target signal s denoting the steering vector are respectively given by where C N represents the N dimensional complex space, i denotes the imaginary unit, f d is the signal normalized Doppler frequency, and (·) T denotes the transpose of matrices or vectors. Conventional SCM estimator, derived by secondary data set {x i } m i=1 , is given by [20] where (·) H represents the conjugate transpose of matrices or vectors. This estimator has been widely utilized; for instance, the GLRT detector reads [27] x and the ANMF reads [13] where γ denotes the detection threshold.
In the matrix-CFAR, to determine whether the received signal x includes a target s or not, we model the observation data by the covariance matrix as a Toeplitz HPD matrix defined as [3,9] where P(N, C) denotes the space of N × N HPD matrices. Its component is the k-th correlation coefficient. Here, r k is the conjugate of r k , and E[·] 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 The covariance matrix of clutter is estimated as R g by the observations {R i } m i=1 . Consequently, the detection problem can be modeled as where the matrix R CUT is the covariance matrix of the observation in the cell under test (CUT). By defining a threshold γ, the signal detection is modeled as discriminating R CUT from R g by where d (R g , R CUT ) is the difference between R CUT and 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 is firstly calculated by (2.7) and (2.9). Mean or median R g of those covariance matrices will be obtained by using either the distance or divergences. The matrix R g can be considered as an estimator of the covariance matrix of the clutter since we consider almost all observation data {x i } m i=1 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 R CUT and R g and the threshold γ.

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.

HPD manifolds and the RD mean and median
We firstly define the AIRM on the HPD matrix manifold P(N, C), and introduce the corresponding geodesic distance, also called the Riemannian distance. The tangent space at a point P ∈ P(N, C) is T P P(N, C) ∼ = H (N, C), which denotes the space of N × N Hermitian matrices. The tangent bundle is denoted by T P(N, C) = ∪ P T P P(N, C). The well-known AIRM is defined by for A, B ∈ T P P(N, C), where tr denotes the trace of a matrix. The geodesic curve between two points P 0 , P 1 ∈ P(N, C) is given by where exp denotes matrix exponential and Log denotes the principal logarithm [21]. Obviously, P (0) = P 0 , and P (1) = P 1 . Under the AIRM metric, the geodesic distance between two points P 0 , P 1 ∈ P(N, C) is given by where · is the Frobenius norm associated to the Frobenius metric A, B := tr A H B . For more details, the reader may refer to [33].
In the Riemannian manifold P(N, C) equipped with the AIRM, the exponential map Exp : T P(N, C) → P(N, C) is defined by making use of the geodesics. Locally, assuming P (t) is the geodesic (3.2) connecting two arbitrary HPD matrices P 0 = P (0), P 1 = P (1), the exponential map gives endpoint of the geodesic, namely whereṖ (0) is the tangent vector at the point P 0 . Using the Riemannian distance (3.3), we recall the RD mean and median of m HPD matrices {R i } m i=1 as follows In other words, the RD mean and median are, respectively, solutions of the optimization problems with objective functions 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]: and Here η t is the step size, and 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].

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 A, B ∈ GL (N, C), the general linear group of N × N invertible complex-valued matrices, is defined by where F (A) is a differentiable and strictly convex function in GL (N, C) and ∇ denotes the (Euclidean) gradient with respect to the Frobenius metric, defined by The well-known total square loss (TSL), and the total log-determinant (TLD) and total von-Neumann (TVN) divergences defined in the HPD manifold P(N, C) are shown in Table 1.
Before introducing the TBD medians, we briefly recall their mean counterparts. Details can be found in [25]. The TBD mean of m HPD matrices {R i } m i=1 is defined by arg min  The function F (Y ) TBDs in P(N, C) As the TBD δ F (R, R i ) is defined in the convex HPD space and is strictly convex with respect to R, the TBD mean exists uniquely. It solves the next matrix equation, 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 {R i } m i=1 are respectively given by 14) Proof. Let G(R) be the function to be optimized, namely The proof is completed by setting ∇G(R) = 0 where the gradient of G(R) is (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 {R i } m i=1 can respectively be calculated by the following fixed-point algorithms in Table 2.
In the matrix-CFAR, the matrix 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 .

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 (0, Σ) of zero-mean and a covariance matrix Σ. Its probability density function is given by 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 Here z is called the speckle component which is fast fluctuating and modeled by a complex Gaussian distribution. The positive real random variable τ is called the texture component which is comparatively slowly varying. In the simulation, the texture component τ is assumed to subject to a gamma distribution with a shape parameter α and a scale parameter β whose probability density function is where Γ(·) 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 (0, Σ) of which the known covariance matrix Σ is given by Here, I N is the N dimensional identity matrix. In the simulation, the clutter-to-noise ratio σ c is set to σ 2 c = 20dB, the one-lag coefficient is ρ = 0.9, and the normalized Doppler frequency of clutter is f c = 0.2. Furthermore, the shape parameter is set to α = 4 and the scale parameter is set to β = 3. The signal to clutter ratio (SCR) is defined by SCR = |ξ| 2 s H R −1 s, (4.5) where ξ is the amplitude coefficient, s is the target vector defined as (2.3) and R is the clutter covariance matrix. Furthermore, the false alarm rate P f a 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 γ is derived by 100/P f a , 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. 2a, 2c, 2e), 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. 2b, 2d, 2f), 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.

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.

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 H (N, C) is given by N 2 number of matrices: where E ij represents the N × 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  H (N, C) can be represented as a linear combination using the orthogonal basis, namely where the summation is over {1, . . . , N 2 } and the Einstein summation convention is applied here and all through the paper. It is immediate to check that the basis {E k } N 2 k=1 is orthogonal by showing that E k , E l = δ kl , ∀k, l = 1, . . . , N 2 , where δ kl is the Kronecker delta that gives 1 when k = l and 0 otherwise. Given a Hermitian matrix H, it is immediate to write down its coordinate expression as

Robustness via influence functions
Let R be the TBD medians (or TBD means, RD mean, RD median) of a given set of m HPD matrices {R i } m i=1 , and let R be the TBD medians (or TBD means, RD mean, RD median) of the mixed HPD data including n number of outliers {P j } n j=1 . By thinking of the outliers as a perturbation ε (ε << 1) of each mean or median, R and R can be related by is a Hermitian matrix depending on the outliers, and the m HPD matrices {R i } m i=1 , the latter of which are assumed to be given and fixed. The function is defined as the (normalized) influence function.
To calculate the corresponding influence functions, define the objective function G(R) for the contaminated m + n number of HPD matrices as follows where d is either the distance or divergence functions for medians, and square of them for means. Since R is any mean or median of the mixed m + n HPD matrices, we have ∇G( R) = 0, 1 that is, Differentiating ∇G( R) = 0 aboutε at ε = 0 , we obtain a matrix equation for H: 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 {E k } N 2 k=1 , we can write and m j=1 ∇d R, P j = φ l E l . (5.12) Letting H = h s E s , the first term in (5.10) then becomes Definition of the gradient (3.11) (induced to HPD spaces) is applied here. Although a linear system about {h k } N 2 k=1 can be obtained, unfortunately, as {ξ l (R)} N 2 l=1 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 } N 2 k=1 (from (5.13)), we can simply Consequently, we have successfully written the matrix equation (5.10) as a linear system about {h k } N 2 k=1 with matrix coefficients: Conducting Frobenius metric on both sides with each E s and making use of the orthogonality, we obtain a linear system for which can be immediately solved if the coefficient matrix (θ s k ) is invertible. We summarize the above deduction into the following theorem. Theorem 6. Under the orthogonal basis {E k } N 2 k=1 defined by (5.2), the Hermitian matrix H = h k 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 H and hence the corresponding influence function. Regarding the contaminated HPD matrices, define G(R) as the objective function, i.e., Differentiate grad G( R) = 0 with respect to ε at ε = 0 and obtain The second term can be simply written as while the first term can be expanded as [25,31] 1 m Substituting H = h k E k to the first term, it can be written in the form 1 m h k θ s k E s , where θ s k is given by Finally, we obtain the linear system (5.16) about {h k } N 2 k=1 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 (0, Σ) where the covariance matrix Σ is given by (4.4). The means and medians R are calculated from the m HPD covariance matrices of the m number of sample data. We then mix n outliers {P j } n j=1 with the sample data; the outliers are modeled as the covariance matrix of ξs + c where s denotes the steering vector, c is the Gaussian clutter with the covariance matrix Σ, and ξ 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 R is calculated from the m number of sample data and the 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.

Conclusions
The TBD medians were defined and applied to matrix-CFAR detectors corresponding to the wellknown 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 [1,15], and structured covariance interference (e.g., [12]).