A Geometric Approach to Covariance Matrix Estimation and its Applications to Radar Problems

TL;DR

This paper introduces a geometric class of covariance matrix estimators for radar signal processing that improve SINR performance, especially in data-scarce scenarios, by projecting sample covariance matrices into structured sets using efficient algorithms.

Contribution

It proposes a novel geometric framework for covariance estimation using unitary invariant norms, with efficient algorithms for Frobenius and spectral cases, enhancing radar processing performance.

Findings

Significant SINR improvements in training-starved regimes.

Efficient algorithms for near closed-form covariance estimates.

Performance validated in spatial and Doppler radar scenarios.

Abstract

A new class of disturbance covariance matrix estimators for radar signal processing applications is introduced following a geometric paradigm. Each estimator is associated with a given unitary invariant norm and performs the sample covariance matrix projection into a specific set of structured covariance matrices. Regardless of the considered norm, an efficient solution technique to handle the resulting constrained optimization problem is developed. Specifically, it is shown that the new family of distribution-free estimators shares a shrinkagetype form; besides, the eigenvalues estimate just requires the solution of a one-dimensional convex problem whose objective function depends on the considered unitary norm. For the two most common norm instances, i.e., Frobenius and spectral, very efficient algorithms are developed to solve the aforementioned one-dimensional optimization leading to almost closed form covariance estimates. At the analysis stage, the performance of the new estimators is assessed in terms of achievable Signal to Interference plus Noise Ratio (SINR) both for a spatial and a Doppler processing assuming different data statistical characterizations. The results show that interesting SINR improvements with respect to some counterparts available in the open literature can be achieved especially in training starved regimes.

Equations

$${\mathbb{E}}[\mbox{\boldmath{$r$}}_{i}\mbox{\boldmath{$r$}}_{i}^{\dagger}]=\mbox{\boldmath{$M$}},\quad\text{for }i=1,\ldots,K,$$

$$\widehat{\mbox{\boldmath{$S$}}}=\frac{1}{K}\sum_{i=1}^{K}\mbox{\boldmath{$r$}}_{i}\mbox{\boldmath{$r$}}_{i}^{\dagger},$$

$$\left\{\begin{array}[c]{ll}&\mbox{\boldmath{$M$}}=\mbox{\boldmath{$R$}}+\sigma_{n}^{2}\mbox{\boldmath{$I$}},\\ &\mbox{\boldmath{$R$}}\succeq{\bf 0},\\ &\sigma_{n}^{2}\geq\sigma^{2},\\ &\frac{\lambda_{max}\left(\mbox{\boldmath{$M$}}\right)}{\lambda_{min}\left(\mbox{\boldmath{$M$}}\right)}\leq\kappa_{M},\end{array}\right.$$

$${\cal{P}}\left\{\begin{array}[c]{ll}\underset{\mbox{\boldmath{$M$}}}{\min}&\begin{array}[c]{c}\|\mbox{\boldmath{$M$}}-\widehat{\mbox{\boldmath{$S$}}}||\end{array}\\ \mbox{s.t.}&\begin{array}[t]{l}\frac{\lambda_{max}\left(\mbox{\boldmath{$M$}}\right)}{\lambda_{min}\left(\mbox{\boldmath{$M$}}\right)}\leq\kappa_{M}\\ \sigma_{n}^{2}\mbox{\boldmath{$I$}}+\mbox{\boldmath{$R$}}=\mbox{\boldmath{$M$}}\\ \mbox{\boldmath{$R$}}\succeq{\bf 0}\\ \sigma_{n}^{2}\geq\sigma^{2}\end{array}\end{array}\right.,$$

$$\displaystyle\|\mbox{\boldmath{$M$}}-\widehat{\mbox{\boldmath{$M$}}}\|$$

$$\le$$

$$\displaystyle{\mathbb{E}}[\|\mbox{\boldmath{$M$}}-\widehat{\mbox{\boldmath{$M$}}}\|^{2}]$$

$$\le$$

$${\cal{P}}_{1}\left\{\begin{array}[c]{ll}\underset{\mbox{\boldmath{$X$}}}{\min}&\begin{array}[c]{c}\|\mbox{\boldmath{$X$}}-\mbox{\boldmath{$S$}}\|\end{array}\\ \mbox{s.t.}&\begin{array}[t]{l}\mbox{\boldmath{$X$}}\succeq\mbox{\boldmath{$I$}}\\ \frac{\lambda_{max}\left(\mbox{\boldmath{$X$}}\right)}{\lambda_{min}\left(\mbox{\boldmath{$X$}}\right)}\leq\kappa_{M}\\ \end{array}\end{array}\right.,$$

$${\cal{P}}_{2}\left\{\begin{array}[c]{ll}\underset{\mbox{\boldmath{$\Lambda$}}}{\min}&\|\mbox{\boldmath{$\Lambda$}}-\mbox{\boldmath{$\Lambda$}}_{S}\|\\ \mbox{s.t.}&\mbox{\boldmath{$\Lambda$}}\succeq\mbox{\boldmath{$I$}}\\ &\frac{\lambda_{max}(\mbox{\boldmath{$\Lambda$}})}{\lambda_{min}(\mbox{\boldmath{$\Lambda$}})}\leq{\kappa_{M}}\end{array}\right.,$$

$${\cal{P}}_{2}^{\prime}(u)\left\{\begin{array}[c]{ll}\underset{\mbox{\boldmath{$\Lambda$}},u}{\min}&\|\mbox{\boldmath{$\Lambda$}}-\mbox{\boldmath{$\Lambda$}}_{S}\|\\ \mbox{s.t.}&\lambda_{i}\geq 1,\\ &u\leq{\lambda_{i}}\leq{\kappa_{M}u},\\ &u\geq{\frac{1}{\kappa_{M}}},\end{array}\right.\quad i=1,\ldots,N.$$

$$\mbox{\boldmath{$\Lambda$}}^{\star}(\bar{u})=\mbox{\boldmath\bf diag}\,(\mbox{\boldmath{$\lambda$}}^{\star}(\bar{u})),$$

$$\begin{split}\mbox{\boldmath{$\lambda$}}^{\star}(u)=&\left[\lambda_{1}(u),\lambda_{2}(u),\ldots,\lambda_{N}(u)\right]^{T}\in\mathbb{R}^{N},\end{split}$$

$${\lambda }_i(u)=\min ({\kappa }_{M}u,\max (d_i,\max (1,u))),i=1,\dots ,N.$$

$${\cal{P}}_{3}\left\{\begin{array}[c]{ll}\underset{u}{\min}&g(h_{1}(u),h_{2}(u),\ldots,h_{N}(u))\\ \mbox{s.t.}&u\geq{\frac{1}{\kappa_{M}}}\\ \end{array}\right.,$$

$$\mbox{\boldmath{$X$}}^{\star}={\mbox{\boldmath{$U$}}}_{S}\,\mbox{\boldmath\bf diag}\,(\mbox{\boldmath{$\lambda$}}^{\star}(u^{\star})){\mbox{\boldmath{$U$}}}_{S}^{\dagger}.$$

$$h_{i}(u)=\left\{\begin{array}[c]{ll}d_{i}-\kappa_{M}u&\mbox{if}\quad\frac{1}{\kappa_{M}}\leq u<\frac{d_{i}}{\kappa_{M}}\\ 0&\mbox{if}\quad\frac{d_{i}}{\kappa_{M}}\leq u<d_{i}\\ u-d_{i}&\mbox{if}\quad u\geq d_{i}\\ \end{array}\right.,$$

$$h_{i}(u)=\left\{\begin{array}[c]{ll}1-d_{i}&\mbox{if}\quad\frac{1}{\kappa_{M}}\leq u<1\\ u-d_{i}&\mbox{if}\quad u\geq 1\\ \end{array}\right..$$

$$g(h_1(u),h_2(u),\dots ,h_N(u))=\sqrt{\sum _{i=1}^N(h_i(u){)}^2}$$

$${\cal{P}}_{4}\left\{\begin{array}[c]{ll}\underset{u}{\min}&G_{1}(u)\\ \mbox{s.t.}&u\geq\frac{1}{\kappa_{M}}\\ \end{array}\right.,$$

$$\mbox{\boldmath{$v$}}=\left[d_{1},d_{2},\ldots,d_{\bar{N}},1\right]^{T}\in\mathbb{R}^{\bar{N}+1},$$

$$u^{\star }=\frac{\sum _{i=1}^{\beta }{\kappa }_M{d}_i+\sum _{i=\alpha }^N{d}_i}{N-\alpha +1+\beta {\kappa }_M^2},$$

$$u_{\alpha ,\beta }=\frac{\sum _{i=1}^{\beta }{\kappa }_M{d}_i+\sum _{i=\alpha }^N{d}_i}{N-\alpha +1+\beta {\kappa }_M^2},$$

$$v_{\alpha }<u_{\alpha ,\beta }\le v_{\alpha -1}\text{mbox}and\frac{v_{\beta +1}}{{\kappa }_M}\le u_{\alpha ,\beta }<\frac{v_{\beta }}{{\kappa }_M},$$

$$g(h_1(u),h_2(u),\dots ,h_N(u))=\underset{i=1,\dots ,N}{\max }\left\{h_i(u)\right\}$$

$$\bar{{\cal{P}}}_{3}\left\{\begin{array}[c]{ll}\underset{u}{\min}&G_{2}(u)\\ \mbox{s.t.}&u\geq\frac{1}{\kappa_{M}}\\ \end{array}\right.,$$

$$\mbox{SINR}_{av}=\frac{1}{MC}\sum_{i=1}^{MC}\frac{|\widehat{\mbox{\boldmath{$w$}}}_{i}^{{\dagger}}\mbox{\boldmath{$s$}}(x)|^{2}}{\left(\widehat{\mbox{\boldmath{$w$}}}_{i}^{{\dagger}}\mbox{\boldmath{$M$}}\widehat{\mbox{\boldmath{$w$}}}_{i}\right)},$$

$$\begin{split}\mbox{\boldmath{$M$}}_{s}\left(n,m\right)=\displaystyle{\sum_{i=1}^{J}}&\sigma_{i}^{2}\mbox{sinc}\left[0.5B_{f}(n-m)\phi_{i}\right]e^{j(n-m)\phi_{i}},\\ \end{split}$$

$$\mbox{\boldmath{$s$}}(\theta)=\left[1,\exp\left(j\pi\sin(\theta)\right),\ldots,\exp\left(j\pi\sin(\theta)(N-1)\right)\right]^{T}.$$

$$\mbox{\boldmath{$r$}}_{i}=\mbox{\boldmath{$n$}}_{i}+\sqrt{\tau_{i}}\mbox{\boldmath{$x$}}_{i}\,\quad i=1,\ldots,K,$$

$$\widehat{\mbox{\boldmath{$M$}}}_{\text{NSCM}}=\frac{N}{K}\sum_{i=1}^{K}\frac{\mbox{\boldmath{$r$}}_{i}\mbox{\boldmath{$r$}}_{i}^{{\dagger}}}{\mbox{\boldmath{$r$}}_{i}^{{\dagger}}\mbox{\boldmath{$r$}}_{i}};$$

Full text

I Introduction

Interference covariance matrix estimation is a longstanding and basic problem in adaptive radar signal processing and naturally arises in several areas such as target detection, direction of arrival estimation, sidelobe cancelling, and secondary data selection [1, 2, 3, 4] (just to list a few). Conventional adaptive architectures (such as Sample Matrix Inversion (SMI) Doppler filter [1], Kelly’s receiver [4], and spatial beamformers [5]) resort to the sample covariance matrix of a secondary data set collected from range gates spatially close to the one under test to estimate the interference covariance. These algorithms are often very prohibitive because they lean on the assumption that the environment remains stationary and homogeneous during the adaptation process. Precisely, they provide satisfactory performance when the secondary vectors share the same spectral properties of the interference in the test cell, are statistical independent, and their number is higher than twice the useful signal dimension [1]. These requisites however may represent important limitations since in real environments the number of data where the disturbance is homogeneous (often referred to as sample support) is very limited. Besides, poor training data selection, in such adaptive algorithms, can result in severe radar performance degradation [6] and [7].

A viable means to thwart the lack of a sufficient number of homogeneous secondary data is to capitalize some a-priori information about the radar surrounding environment, namely to realize a knowledge-based/cognitive processing [8] so as to restrict the uncertainty region of the unknown parameters. According to this processing paradigm, several approaches have been pursued in the open literature assuming different structural models as well as statistical distributions of the data. In particular, both homogeneous and heterogeneous interference environments are dealt with. As to the former scenario, the training data are modeled as independent and identically distributed (i.i.d.), zero-mean, circularly symmetric Gaussian vectors. Whereas, the latter context mainly considers clutter power variations within the sample support. Thus, assuming homogeneity, in [9] the Maximum Likelihood (ML) covariance matrix estimator is derived modeling the disturbance as the sum of a coloured interference plus white disturbance; in [10], the ML estimation of an unstructured covariance matrix with a condition number upper bound requirement is considered; in [11], using the same covariance structure as in [9], the ML estimator is derived when a constraint on the condition number is imposed too; in [12], a rank-constrained ML estimator is developed; furthermore, relying on a Mean Square Error (MSE) design criterion, in [13] and [14], some shrinkage estimators are proposed. With reference to heterogeneous scenarios, compound Gaussian statistical models (such as K-distributed or Gamma amplitudes) are usually exploited to account for clutter returns spikiness. In [15], assuming a clutter dominated environment, a unified framework to regularize the ML estimate in scaled Gaussian models (e.g., elliptical distributions, compound-Gaussian processes and spherically invariant random vectors) is developed in order to enhance Tyler’s estimator in the presence of a small sample support exploiting a-priori information on the covariance structure. In [16], Tyler’s robust covariance M-estimator under group symmetry constraints, such as circulant, persymmetric, and proper quaternion matrices, is considered. Precisely, it is provided an iterative fixed point algorithm to compute the constrained estimate. Finally, in [17], an iterative algorithm to estimate, according to the ML approach, both the clutter subspace and the covariance is proposed, assuming the disturbance composed of a low rank compound Gaussian clutter plus a white Gaussian noise contribution. Further technically sound and effective covariance estimators can be found in [18, 19, 20, 21, 22, 23].

These mentioned algorithms lean on specific assumptions about the data statistical characterization and may suffer performance loss in the presence of model mismatches due to for instance quantization effects, constant modulus jamming signals, phase noise [24, 25, 26]. To overcome this shortcoming and endow robustness to the estimation process, some strategies derived from geometric considerations on the metric space of the covariance matrices have been also developed. In this respect, in [27] and [28] techniques based on Riemannian $p$-mean (e.g., Fréchet median and Karcher barycenter) evaluation are developed to estimate suitably structured interference covariance matrices showing that substantial ameliorations over classical algorithms can be attained. In [27] and [29] also an extension of the conventional Ordered Statistic (OS) framework [30] is proposed relying on the Riemannian $p$-mean computation of Toepliz or Toeplitz-Block-Toeplitz space-time covariance matrices. Besides, in [31] and [32], covariance estimates defined through geometric barycenters/medians (associated with specific distances in the space of Hermitian matrices) of structured covariance estimates are exploited both for training data selection and adaptive radar detection highlighting significant gains with respect to the classic sample covariance matrix. Finally, in [33, 34, 35, 36, 37, 38] other interesting geometric-inspired procedures are devised.

In this paper, leveraging on a geometric criterion, a novel class of covariance estimators that do not consider any assumption on the statistical characterization of the secondary data is proposed and analyzed. Each estimator is associated with a given unitary invariant norm and performs the sample covariance matrix projection into a specific set of structured covariance matrices. Precisely, this set encompasses the matrices modeled as the sum of an unknown positive semi-definite matrix (describing coloured interference and clutter) plus a term proportional to the identity matrix (related to white disturbance). Besides, a constraint on the condition number is accounted for so as to control the numerical stability of the resulting adaptive algorithms [5]. Regardless of the considered norm, an efficient solution technique to handle the formulated constrained optimization problem is developed. Precisely, each estimator exhibits a shrinkage-type form and its evaluation requires the sample covariance matrix spectral decomposition as well as the solution of a one-dimensional convex problem whose objective function depends on the considered unitary norm.

At the analysis stage, the performance of the new class of distribution-free estimators is evaluated in terms of achievable Signal-to-Interference-plus-Noise Ratio (SINR) for different sample support sizes, assuming both a spatial and a Doppler processing scenario. The results show that interesting SINR improvements can be achieved with respect to some counterparts available in the open literature also with a reduced computational complexity.

The remainder of this paper is organized as follows. Section II is devoted to the description of the system model as well as the formulation of the covariance matrix estimation problem. In Section III, an efficient procedure to solve the resulting constrained optimization is developed. In Section IV, the performance of the proposed distribution-free estimators is assessed. Finally, Section V concludes the paper and provides some possible future research tracks.

Notation

We adopt the notation of using boldface for vectors $a$ (lower case), and matrices $A$ (upper case). The transpose and the conjugate transpose operators are denoted by the symbols $(\cdot)^{T}$ and $(\cdot)^{\dagger}$ respectively. $\mbox{\rm tr}\,\{\cdot\}$ is the trace of the square matrix argument. $I$ and ${\bf 0}$ denote respectively the identity matrix and the matrix with zero entries (their size is determined from the context). \mbox{\boldmath\bf diag}\,(\mbox{\boldmath{a}}) indicates the diagonal matrix whose $i$-th diagonal element is the $i$-th entry of $a$. ${\mathbb{R}}^{N}$, ${\mathbb{C}}^{N}$, ${\mathbb{C}}^{N,K}$, and ${\mathbb{H}}^{N}$ are respectively the sets of $N$-dimensional vectors of real numbers, of $N$-dimensional vectors of complex numbers, of $N\times K$ matrices of complex numbers, and of $N\times N$ Hermitian matrices. The curled inequality symbol $\succeq$ (and its strict form $\succ$) is used to denote generalized matrix inequality: for any \mbox{\boldmath{A}}\in{\mathbb{H}}^{N}, \mbox{\boldmath{A}}\succeq\mbox{\boldmath{0}} means that $A$ is a positive semi-definite matrix (\mbox{\boldmath{A}}\succ\mbox{\boldmath{0}} for positive definiteness). $\|\cdot\|$ denotes an arbitrary unitary invariant matrix norm operator, while the specific spectral and Frobenius instances are indicated by $\|\cdot\|_{2}$ and $\|\cdot\|_{F}$, respectively. The letter $j$ represents the imaginary unit (i.e. $j=\sqrt{-1}$). For any complex number $x$, $|x|$ represents the modulus of $x$. Finally, ${\mathbb{E}}\left[\cdot\right]$ denotes statistical expectation and for any optimization problem ${\cal{P}}$ $v({\cal{P}})$ represents its optimal value.

II Problem Formulation

In this section, the problem of estimating the covariance matrix \mbox{\boldmath{M}}\in\mathbb{H}^{N} of $K$ secondary data \mbox{\boldmath{r}}_{1},\ldots,\mbox{\boldmath{r}}_{K} modeled as $N$-dimensional, circularly symmetric, zero-mean random vectors, is addressed. It is assumed that these vectors share the same second order statistical characterization, i.e.,

[TABLE]

but are drawn from an arbitrary and unknown joint probability distribution.

According to the previous assumptions, classic estimation approaches such as the ML or the minimum MSE strategies can be no longer pursued since they need the data statistical distribution knowledge. Hence, in this work a new family of covariance estimators based on geometric considerations is introduced. Specifically, the idea is to estimate the data covariance matrix performing the projection111Notice that, the projection operator is usually defined assuming the reference set being convex and closed as well as the norm being induced by an inner product. Nevertheless, with a slight abuse of notation, in this paper we continue to define as projection the point \mbox{\boldmath{x}}_{p} minimizing the distance (based on a specific metric) between a given point $x$ and the reference set as long as the point \mbox{\boldmath{x}}_{p} can be uniquely identified., as induced by a specific metric, of the sample covariance matrix

[TABLE]

into the set

[TABLE]

where $\sigma^{2}>0$ is a lower bound to the white interference power, $R$ accounts for the colored interfering contribution, and $\kappa_{M}\geq 1$ is an upper bound to the covariance condition number. This uncertainty set accounts for an interference covariance structure that is commonly met in adaptive radar signal processing applications and ensures a well conditioned estimate necessary to compute the adaptive radar weight vector. The effectiveness of this model has already been proved in [11], where the training data are assumed i.i.d., zero-mean, circularly symmetric Gaussian vectors and the ML estimate is derived. Hence, the main goal of this work is the development of a distribution-free approach which provides a robust alternative to [11] when inference on the interference statistics is not possible. A pictorial representation of the geometric-based estimation process is reported in Fig. 1.

Following the above guidelines, the covariance matrix estimate \widehat{\mbox{\boldmath{M}}} is obtained as solution to the following optimization problem

[TABLE]

where $\|\cdot\|$ refers to an arbitrary but given unitary invariant matrix norm [39] that induces a specific metric in the space of the positive semi-definite matrices over which performing the projection. Examples of norms that can be considered in the estimation process are the Frobenius, the spectral, and the Ky Fan norms. In this respect, observe that the first two instances are the most common and widely used norms in the space of Hermitian matrices corroborating the interest toward the family of unitary invariant norms. Remarkably, different estimators can be jointly exploited within a bank/battery of adaptive receivers. More in detail, each detector may resort to a specific norm and the presence of a prospective target is established by means of a suitable fusion logic, such as binary integration or $m$-of-$n$ detection [40].

Before concluding this section, it is also worth pointing out that the proposed estimators possess the consistency property as long as the secondary data vectors are statistically independent. Indeed,

[TABLE]

where the first equation stems from the triangular inequality, while the second inequality follows from the definition of \widehat{\mbox{\boldmath{M}}}, i.e., it is a minimizer. Now, provided that \widehat{\mbox{\boldmath{M}}} is a measurable function of the secondary data222A sufficient condition for estimator measurability is reported in the supplementary material., it follows that

[TABLE]

where the second inequality results from the finite dimension of the space $\mathcal{H}^{N}$, $\alpha$ is a specific constant linking the considered norm with the Frobenius one [39, Corollary 5.4.5], and finally the convergence to zero comes from the consistency of the sample covariance estimator. Based on Chebyshev’s inequality [41], the estimation error converges to zero also in probability which is often invoked as the classic definition of consistency.

III Derivation of the Structured Covariance Matrix Estimators Based on Unitary Invariant Norm Projection

This section deals with the development of an efficient procedure to solve Problem ${\cal{P}}$ almost in closed form regardless of the considered norm. Specifically, it is proved that any sample covariance projector belongs to the class of the shrinkage estimators333A shrinkage covariance estimator \widehat{\mbox{\boldmath{M}}} is a matrix sharing the same eigenvectors as the sample covariance matrix \mbox{\boldmath{S}}\succeq{\bf 0}, but transforming the eigenvalues, i.e. \widehat{\mbox{\boldmath{M}}}=\mbox{\boldmath{U}}_{S}\,\mbox{\boldmath\bf diag}\,([g_{1}(d_{1},d_{2},\ldots,d_{N}),g_{2}(d_{1},d_{2},\ldots,d_{N}),\ldots,\\ g_{N}(d_{1},d_{2},\ldots,d_{N})])\,\mbox{\boldmath{U}}_{S}^{\dagger}\succeq{\bf 0}. and the eigenvalues estimate is obtained solving a one-dimensional convex problem whose objective function is tied up to the considered unitary norm. In general, this optimization problem can be solved in polynomial-time using convex solvers such as CVX [42]. In addition, for the two most relevant norm instances, i.e., Frobenius and spectral, very efficient algorithms are provided to tackle the associated one-dimensional optimizations.

As first step toward the solution of ${\cal{P}}$, let us observe that it is equivalent to the solvable convex problem

[TABLE]

where444The interested reader may refer to the supplementary material for the proof. \mbox{\boldmath{S}}=\widehat{\mbox{\boldmath{S}}}/\sigma^{2} and \mbox{\boldmath{X}}=\mbox{\boldmath{M}}/\sigma^{2}. Next, let us indicate with {\mbox{\boldmath{S}}}={\mbox{\boldmath{U}}}_{S}{\mbox{\boldmath{\Lambda}}}_{S}{\mbox{\boldmath{U}}}^{{\dagger}}_{S} the spectral decomposition of the normalized sample covariance matrix, where {\mbox{\boldmath{\Lambda}}}_{S}=\mbox{\boldmath\bf diag}\,(\left[d_{1},d_{2},\cdots,d_{N}\right]^{T}), with $d_{1}\geq d_{2}\geq\cdots\geq d_{N}$ the eigenvalues of $S$ arranged in decreasing order, and {\mbox{\boldmath{U}}}_{S} is the unitary matrix whose columns contain the corresponding eigenvectors. Hence, the following lemma holds true.

Lemma III.1

An optimal solution to ${\cal{P}}_{1}$ is a shrinkage estimator \mbox{\boldmath{X}}^{\star}={\mbox{\boldmath{U}}}_{S}{\mbox{\boldmath{\Lambda}}}^{\star}{\mbox{\boldmath{U}}}_{S}^{{\dagger}} where \mbox{\boldmath{\Lambda}}^{\star}=\mbox{\boldmath\bf diag}\,\left(\left[\lambda_{1}^{\star},\lambda_{2}^{\star},\cdots,\lambda_{N}^{\star}\right]^{T}\right) is a solution to the following optimization problem

[TABLE]

with \mbox{\boldmath{\Lambda}}=\mbox{\boldmath\bf diag}\,\left(\left[\lambda_{1},\lambda_{2},\cdots,\lambda_{N}\right]^{T}\right).

Proof:

See Appendix -A. ∎

To proceed further, let us introduce the auxiliary variable $u>0$ and cast ${\cal{P}}_{2}$ as555The interested reader may refer to the supplementary material for the proof.

[TABLE]

This formulation paves the way for an efficient solution of ${\cal{P}}_{2}$. Indeed, for any fixed $u$ a closed form optimal matrix can be derived as shown in the following lemma.

Lemma III.2

For any $\bar{u}\geq\frac{1}{\kappa_{M}}$, an optimal solution \mbox{\boldmath{\Lambda}}^{\star}(\bar{u}) to ${\cal{P}}_{2}^{\prime}(\bar{u})$, is

[TABLE]

where

[TABLE]

with

[TABLE]

Proof:

See Appendix -B. ∎

Leveraging on Lemmas III.1 and III.2, and denoting by $g(\cdot)$ the gauge function associated with the considered unitary invariant norm [39], the following fundamental result can be shown.

Theorem III.3

Let $u^{\star}$ be the lowest optimal solution to the convex optimization problem

[TABLE]

where $h_{i}(u)=\left|\lambda_{i}^{\star}(u)-d_{i}\right|$, $i=1,\ldots,N$. Then, an optimal solution to ${\cal{P}}_{1}$ is

[TABLE]

Proof:

See Appendix -C. ∎

According to Theorem III.3, a unique solution to ${\cal{P}}_{1}$ can be constructed in almost closed form. Consequently, this new class of estimators effectively performs specific unitary invariant norm-based projections. It is also worth pointing out that the functional dependence over the selected norm is concentrated in the optimal value of the auxiliary $u$ that accounts for the corresponding gauge function. In Fig. 2, a schematic illustration of the steps involved in the proposed procedure for the computation of \widehat{\mbox{\boldmath{M}}} is reported. As already highlighted, \widehat{\mbox{\boldmath{M}}} is a shrinkage estimator which regularizes the sample covariance matrix according to the specific unitary invariant norm and explicitly accounting for a condition number constraint so as to provide a well-conditioned structured estimate.

Before continue further, it is worth observing that the functions $h_{i}(u)$, $i=1,\ldots,N$, involved in ${\cal{P}}_{3}$ can be expressed in closed form666The interested reader may refer to the supplementary material for the proof.. This is a useful result in general and more important it lays the ground for the derivation of specialized procedures to handle ${\cal{P}}_{3}$ assuming either the Frobenius or the spectral norm (see Subsections III-A and III-B). Precisely, if $d_{i}>1$

[TABLE]

otherwise $d_{i}\leq 1$ and

[TABLE]

III-A Frobenius Norm

The gauge function associated with the Frobenius norm is given by

[TABLE]

and the problem to solve boils down to

[TABLE]

where $G_{1}(u)=\displaystyle{\sum_{i=1}^{N}}{(h_{i}(u))^{2}}$. To proceed further, let us observe that the optimal solution to problem ${\cal{P}}_{4}$ is $u^{\star}=\frac{\max\{1,d_{1}\}}{\kappa_{M}}$ if $d_{1}\leq\kappa_{M}$, otherwise, the optimal solution to ${\cal{P}}_{4}$ lies within the interval777The interested reader may refer to the supplementary material for the proof. $[1,d_{1}]$. In this last case, as shown in supplementary material, $G_{1}(u)$ is a convex function with a continuous derivative function within the interval of interest $\left[1,d_{1}\right]$. Hence, denoting by $\bar{N}$ the number of sample eigenvalues $d_{i}$ greater than $1$, i.e. $d_{i}>1$, $i=1,\ldots,\bar{N}$, and defining the vector

[TABLE]

the following theorem holds true.

Theorem III.4

Assuming $d_{1}>\kappa_{M}>1$, the optimal solution $u^{\star}$ to ${\cal{P}}_{4}$ is

$u^{\star}=1$, if \frac{dG_{1}(u)}{du}\Big{|}_{u=1}\geq 0; 2. 2.

$u^{\star}=d_{1}$, if \frac{dG_{1}(u)}{du}\Big{|}_{u=d_{1}}\leq 0; 3. 3.

$u^{\star}=\frac{d_{1}}{\kappa_{M}}$, if $\frac{d_{1}}{\kappa_{M}}\leq d_{N}$; 4. 4.

if $1)$, $2)$, and $3)$ are not satisfied, $u^{\star}$ is the optimal solution if and only if

[TABLE]

with $\alpha\in\{1,2,\ldots,\bar{N},\bar{N}+1\}$ the smallest index such that $v_{\alpha}<u^{\star}$, and $\beta\in\{1,2,\ldots,\bar{N},\bar{N}+1\}$ the largest index such that $\frac{v_{\beta}}{\kappa_{M}}>u^{\star}$.

Proof:

See Appendix -D. ∎

Theorem III.4 provides the guidelines to find $u^{\star}$. Indeed, the selection of the integers $\alpha$ and $\beta$ such that

[TABLE]

with

[TABLE]

is required. To this end, an efficient procedure is now described. Precisely, the strategy consists in iteratively verifying the conditions (20) and (21) once the values of $\alpha$ and $\beta$ have been efficiently fixed. In this respect, notice that, if the intersection of the intervals (21) is empty, then $u_{\alpha,\beta}$ cannot be the optimal solution. On the contrary, the intersection is given by one of the following sub-intervals \Big{]}v_{\alpha},v_{\alpha-1}\Big{]},\left]v_{\alpha},\frac{v_{\beta}}{\kappa_{M}}\right[,\left[\frac{v_{\beta+1}}{\kappa_{M}},v_{\alpha-1}\right],\left[\frac{v_{\beta+1}}{\kappa_{M}},\frac{v_{\beta}}{\kappa_{M}}\right[, and the optimal point must belong to one of them. According to this line of reasoning, the procedure follows these simple steps:

Set $\beta=1$, $\alpha=2$, and increase $\alpha$ until $v_{\alpha}\geq\frac{v_{\beta}}{\kappa_{M}}$. 2. 2.

Compute $u_{\alpha,\beta}$. If $u_{\alpha,\beta}$ belongs to the current intersection, let $u^{\star}=u_{\alpha,\beta}$ and exit; otherwise go to step $3)$. 3. 3.

If $\frac{v_{\beta+1}}{\kappa_{M}}<v_{\alpha}$, increase $\alpha$ and go to step $2)$, otherwise increase $\beta$ and go to step $2)$.

It is worth pointing out that this algorithm provides the optimal solution to problem ${\cal{P}}_{4}$ with a linear computational complexity with the number of the sample covariance matrix eigenvalues greater than $1$.

III-B Spectral Norm

In this case, the gauge function is

[TABLE]

and ${\cal{P}}_{3}$ can be specialized as

[TABLE]

with $G_{2}(u)=\underset{i}{\max}\left\{h_{i}(u)\right\}$.

The following proposition provides an efficient procedure to solve $\bar{{\cal{P}}}_{3}$.

Theorem III.5

Let $u^{\star}$ the lowest optimal solution to $\bar{{\cal{P}}}_{3}$. Then

$u^{\star}=\frac{1}{\kappa_{M}}$, if $d_{1}\leq 1$. 2. 2.

$u^{\star}=\max\left\{\frac{d_{1}+d_{N}-1}{\kappa_{M}},\frac{1}{\kappa_{M}}\right\}$, if $1<d_{1}\leq\kappa_{M}$ and $d_{N}\leq 1$. 3. 3.

$u^{\star}=\frac{d_{1}}{\kappa_{M}}$, if $1<d_{1}\leq\kappa_{M}$ and $d_{N}>1$. In this case, the covariance estimate reduces to $S$. 4. 4.

If $d_{1}>\kappa_{M}$ and $d_{N}\leq 1$, then

• •

$u^{\star}=\max\left\{\eta_{1},\frac{1}{\kappa_{M}}\right\}$, if $\eta_{1}=\frac{d_{1}+d_{N}-1}{\kappa_{M}}\leq 1$;

• •

otherwise, $u^{\star}=\frac{d_{1}+d_{N}}{1+\kappa_{M}}>1$. 5. 5.

If $d_{1}>\kappa_{M}$ and $d_{N}>1$, then

• •

$u^{\star}=\frac{d_{1}+d_{N}}{1+\kappa_{M}}$. if $d_{N}\leq\frac{d_{1}}{\kappa_{M}}$;

• •

otherwise, $u^{\star}=\frac{d_{1}}{\kappa_{M}}$. Besides, in this case, the estimate coincides with $S$.

Proof:

See Appendix -E. ∎

Based on Theorem III.5, the optimal solution to $\bar{{\cal{P}}}_{3}$ is substantially available in closed form. Indeed, just the comparison between some linear functions of the highest and the lowest sample covariance eigenvalues with some fixed thresholds is required. It is also worth pointing out that, $u^{\star}$ is a continuous function of $d_{1},\ldots,d_{N}$ implying that \mbox{\boldmath{X}}^{\star} is a continuous function of $S$.

IV Performance Analysis

This section is devoted to the analysis of the proposed covariance estimators. The average SINR is adopted as performance metric and some counterparts available in the open literature are considered for comparison purposes. Two typical radar signal processing scenarios are studied: the former focuses on spatial processing with wideband jammers impairing the received data, the latter considers Doppler processing with the interfering returns originated by clutter. More formally, the average SINR (over $MC$ i.i.d. realizations888In the numerical results the number of Monte Carlo trials is set to $MC=500$. of $K$ secondary data) is given by

[TABLE]

where \mbox{\boldmath{s}}(x) is the $N$-dimensional target steering vector whose expression depends on:

• •

the considered processing scenario (spatial/temporal);

• •

the radar configuration (array type, Pulse Repetition Time (PRT), etc.);

• •

the target state $x$ (angle of arrival $\theta$/normalized Doppler frequency $\nu$).

Moreover, \widehat{\mbox{\boldmath{w}}}_{i}=\widehat{\mbox{\boldmath{M}}}_{i}^{-1}\mbox{\boldmath{s}}(x) is the adaptive estimate of the optimal weight vector, where \widehat{\mbox{\boldmath{M}}}_{i} is the data-dependent estimate of $M$ at the $i$-th run.

In the following analysis, it is assumed $\sigma^{2}=0$ dB and \kappa_{M}=\lambda_{max}(\mbox{\boldmath{M}})/\lambda_{min}(\mbox{\boldmath{M}}). Furthermore, for each case study two different values of the actual white noise power level are considered, i.e., $\sigma_{a}^{2}=0$ and $\sigma_{a}^{2}=10$ dB, so as to account for both a matched and a mismatched scenario.

IV-A Spatial Processing

A radar system equipped with a uniform linear array of $N=8$ elements (with a spacing between the antennas of $d=\lambda_{0}/2$ where $\lambda_{0}$ is the radar operating wavelength) pointing in the boresight direction is considered. The interference covariance matrix is given by \mbox{\boldmath{M}}=\mbox{\boldmath{M}}_{s}+\sigma_{a}^{2}\mbox{\boldmath{I}} [9] where $\sigma_{a}^{2}$ is the actual power level of the white disturbance term, whereas \mbox{\boldmath{M}}_{s} is the covariance matrix associated to $J$ (possibly wideband) jammers. Specifically, $\forall\left(n,m\right)\in\left\{1,\ldots,N\right\}^{2}$,

[TABLE]

with $B_{f}=B/f_{0}$ the fractional bandwidth, $B$ the instantaneous bandwidth of the desired signal (coinciding with the jammer’s bandwidth), $f_{0}=c/\lambda_{0}$, $\sigma_{i}^{2}$ the power associated with the $i$-th jammer, and $\phi_{i}$ the jammer phase angle with respect to the antenna phase center. Precisely, $\phi_{i}=2\pi d(\sin\theta_{i})/\lambda_{0}$, with $\theta_{i}$ the angle off-boresight of the jammer. Finally, according to the specified system model, the steering vector (23) reduces to

[TABLE]

As case study, it is considered a wideband jammer with a fractional bandwidth $B_{f_{1}}=0.3$, a power $\sigma_{1}^{2}=30$ dB, and a direction of arrival $\theta_{1}=20$ deg that impinges of the radar receive array. In Fig. 3, the average SINR is plotted versus the Direction Of Arrival (DOA) $\theta$ for both the Frobenius Norm based Estimator (FNE) and the Spectral Norm based Estimator (SNE). Therein, the secondary data are modeled as i.i.d., zero-mean, circularly symmetric Gaussian random vectors. For comparison purposes, the SINR behavior associated with the Constrained ML estimator (CML) [11], the Fast Maximum Likelihood estimator (FML) [9], and the classic Sample Covariance Matrix (SCM)999Notice that, when $K<N$ the pseudo inverse of the sample matrix is utilized in place of its inverse [1]. is displayed too. Besides, the SINR upper bound \mbox{\boldmath{s}}^{{\dagger}}\mbox{\boldmath{M}}^{-1}\mbox{\boldmath{s}} is reported as benchmark. Figs. 3(a), 3(c), 3(e) assume a matched condition, i.e. $\sigma_{a}^{2}=\sigma^{2}=0$ dB, with $K=4$, $K=8$, and $K=16$ secondary data, respectively. Figs. 3(b), 3(d), 3(e) account for a mismatched situation, i.e. $\sigma_{a}^{2}=10$ dB and $\sigma^{2}=0$ dB.

The results show that FNE and SNE substantially exhibit the same performance in terms of average SINR regardless of the considered DOA and scenario. In the matched cases the curves of the new devised estimators almost overlap with those of FML and CML (the maximum gain of FNE and SNE over CML and FML is 0.08 dB) and significantly dominate the SCM performance. Besides, FNE and SNE outperform all the counterparts in the mismatched scenario with gains (with respect to CML, FML, and SCM) of: $0.93$ dB, $1.26$ dB, and $5.65$ dB for $K=4$; $0.80$ dB, $2.40$ dB, and $5.45$ dB for $K=8$; $0.45$ dB, $1.60$ dB, and $1.60$ dB for $K=16$. As expected, increasing the sample support size the gain reduces since all the curves tend to approach the SINR upper bound due to the consistency of the estimators.

In Fig. 4, the same spatial processing scenario as in Fig. 3 is analyzed, but for a different secondary data statistical distribution that is no longer Gaussian. Specifically,

[TABLE]

where \mbox{\boldmath{n}}_{i}\sim\mathcal{CN}(0,\sigma_{a}^{2}\mbox{\boldmath{I}}), \mbox{\boldmath{x}}_{i}\sim\mathcal{CN}(0,\mbox{\boldmath{M}}_{s}), and $\tau_{i}\sim\Gamma(1/\mu_{\tau},\mu_{\tau})$ ($\mu_{\tau}=2$), $i=1,\ldots,K$, are statistically independent random variables/vectors. Otherwise stated, a compound Gaussian jamming is now accounted for. As in Fig. 3, the average SINR versus $\theta$ is reported for three different sample support sizes, i.e., $K=4$, $K=8$, and $K=16$. Moreover, the subplots on the left refer to the matched scenario ($\sigma_{a}^{2}=\sigma^{2}=0$ dB) whereas the plot on the right address the mismatched case ($\sigma_{a}^{2}=10$ dB and $\sigma^{2}=0$ dB).

For comparison purposes, other than the FML, CML, and SCM estimators, three additional strategies designed to operate in compound Gaussian clutter are considered. Specifically,

• •

the Normalized Sample Covariance Matrix estimator (NSCM), [43]

[TABLE]

• •

the Fixed-Point Estimator (FPE) [44, 45, 46] \widehat{\mbox{\boldmath{M}}}_{\text{FPE}}, that is obtained iteratively solving a fixed point equation;

• •

the Low Rank clutter Estimator (LRE), [17], addressing a mixed Gaussian/compound Gaussian disturbance model. The covariance estimate can be computed as

[TABLE]

where $\widehat{\tau}_{k}$ is the estimated texture of the $k$-th clutter datum and \widehat{\mbox{\boldmath{\Sigma}}} is the covariance estimate of the speckle obtained through the iterative algorithm proposed in [17].

Inspection of Fig. 4 reveals that FNE and SNE are basically equivalent and outperform in terms of average SINR all the counterparts, included those specific for compound Gaussian disturbance101010As to the LRE, the clutter covariance matrix rank is evaluated as the number of the eigenvalues greater than \mbox{\rm tr}\,(\mbox{\boldmath{M}}_{s})/10^{5}\geq 10^{-4}. Hence, the estimators exploiting the true rank, i.e., 6 (LRE-6) and rank 7 (LRE-7) are displayed.. Precisely, in the matched condition, as already observed in Fig. 3, FNE, SNE, FML, and CML are almost coincident. Besides, compared to NSCM, FPE, LRE-6, LRE-7, and SCM, they provide SINR gains up to: $11$ dB, $11$ dB, $1.89$ dB, $1.90$ dB, and $6.93$ dB for $K=4$; $6.27$ dB, $6.27$ dB, $4.74$ dB, $7.50$ dB, and $5.38$ dB for $K=8$; $4.77$ dB, $2.01$ dB, $6.65$ dB, $9.20$ dB, $1.63$ dB for $K=16$. Without surprise, LRE-6 outperforms LRE-7 reflecting the presence of a mismatch loss. Finally, in the mismatched scenario, FNE and SNE also grant better performance than FML and CML with gains of $1.35$ dB and $0.94$ dB for $K=4$, $2.35$ dB and $0.84$ dB, for $K=8$, $1.51$ dB and $0.42$ dB, for $K=16$, clearly highlighting the effectiveness of the new devised strategies. As to the comparisons with the other estimators, considerations similar to those for the matched scenario holds true.

IV-B Doppler Processing

A radar system transmitting a coherent burst of $N=16$ pulses is considered. In this case $x$ refers to the normalized Doppler frequency of the target, i.e. $\nu\in[-1/2,1/2[$ and the steering vector in (23) reduces to

[TABLE]

As to the interference environment, it is assumed that the radar operates in the presence of both ground and sea clutter in addition to white noise [47]. Therefore the overall disturbance covariance matrix is

[TABLE]

where, $\forall\left(n,m\right)\in\left\{1,\ldots,N\right\}^{2}$,

[TABLE]

with

• •

$\text{CNR}_{S}$ and $\text{CNR}_{G}$ the power of the sea and ground clutter, respectively;

• •

$\rho_{S}$ and $\rho_{G}$ the one-lag correlation coefficients of the sea and the ground clutter, respectively;

• •

$f_{S}$ the normalized Doppler frequency of the sea clutter.

In Fig. 5, the average SINR versus $\nu$ is reported for FNE, SNE, FML, CML, and SCM estimators. The training data are drawn from a complex circular Gaussian distribution and the parameters in (26) are $\text{CNR}_{S}=10$ dB, $\text{CNR}_{G}=25$ dB, $\rho_{S}=0.8$, $\rho_{G}=0.95$, and $f_{S}=0.2$. Figs. 5(a), 5(c), and 5(e) refer to $K=8$, $K=16$, and $K=32$, respectively and assume $\sigma_{a}^{2}=0$ dB. The mismatched analysis, i.e., $\sigma_{a}^{2}=10$ dB, is instead reported in Figs. 5(b), 5(d), and 5(f).

The plots clearly illustrate the effectiveness of the new devised estimators. Indeed, both in the matched and in mismatched scenario, FNE and SNE achieve higher SINR values than the counterparts at each target Doppler frequency. Interestingly, unlike the spatial-processing, in this case FNE and SNE also outperform CML and FML in the matched conditions. Precisely, the SINR gains with respect to the CML, that is the major competitor, are, for $K=8,16,32$, respectively: $1.29$ dB, $1.38$ dB, and $0.52$ dB in the matched case and $1.43$ dB, $1.42$ dB, and $0.61$ dB in mismatched situation. As expected, the gains are lower and lower as $K$ increases, due to the consistency of all the involved estimators.

V Conclusions

The disturbance covariance matrix estimation problem for radar signal processing applications has been addressed according to a geometric approach. Specifically, a new family of distribution-free covariance estimators has been introduced performing the sample covariance matrix projection, according to a specific unitary invariant norm, into a structured covariance set of practical relevance.

To tackle the resulting constrained optimization problem an efficient solution technique has been designed which represents the main technical contribution of this paper. In particular, it has been proved that each estimator exhibits a shrinkage-type form with the eigenvalues estimate obtained via the solution of a one-dimensional convex problem tuned to the considered unitary norm. Furthermore, almost closed form covariance estimates have been provided assuming either Frobenius or spectral norms at the design stage. Remarkably, the proposed estimators possess the consistency property as long as the training vectors are statistically independent.

Some interesting case studies have been considered to illustrate the effectiveness of the new proposed framework. The results have shown that the new estimators may provide better SINR values than some structured estimators available in the open literature. Precisely, the lower the sample support the higher the gain. Additionally, accounting for the computational efforts as well as the SINR performance, the best estimator appears that based on spectral norm.

As possible future research tracks, it might be worth analyzing the performance of the new family of estimators on real radar data as well as to account for other unitary invariant constraint at the design stage such as an upper bound on the clutter rank.

-A Proof of Lemma III.1

Proof:

Let {\mbox{\boldmath{X}}}_{1}^{\star}={\mbox{\boldmath{U}}}_{1}{\mbox{\boldmath{\Lambda}}}_{1}^{\star}{\mbox{\boldmath{U}}}^{{\dagger}}_{1} be the spectral decomposition of an optimal solution to ${\cal{P}}_{1}$, where {\mbox{\boldmath{\Lambda}}}_{1}^{\star}=\mbox{\boldmath\bf diag}\,(\lambda_{1}^{\star},\lambda_{2}^{\star},\cdots,\lambda_{N}^{\star}), with $\lambda_{1}^{\star}\geq\lambda_{2}^{\star}\geq\cdots\geq\lambda_{N}^{\star}$. Based on [39, Theorem 7.4.51]

[TABLE]

Besides,

[TABLE]

due to the unitary invariance of the norm. Hence, \mbox{\boldmath{X}}^{\star}={\mbox{\boldmath{U}}}_{S}{\mbox{\boldmath{\Lambda}}}_{1}^{\star}{\mbox{\boldmath{U}}}_{S}^{{\dagger}} is an optimal solution to ${\cal{P}}_{1}$ since it is a feasible point achieving the minimum value. Finally, ${\cal{P}}_{2}$ is obtained replacing $X$ with {\mbox{\boldmath{U}}}_{S}{\mbox{\boldmath{\Lambda}}}{\mbox{\boldmath{U}}}_{S}^{{\dagger}} in ${\cal{P}}_{1}$. ∎

-B Proof of Lemma III.2

Proof:

Let $g(\cdot)$ be the gauge function associated with the considered unitary invariant norm $\|\cdot\|$ (see [39, Theorem 7.4.24]). According to [39, Property 7.4.21],

[TABLE]

implying that ${\cal{P}}_{2}^{\prime}({u})$ can be cast as

[TABLE]

To proceed further, let us observe that $g(\cdot)$ is a monotone norm on $\mathbb{C}^{N}$ [39, Theorem 5.5.10]. As a consequence, given two vectors \mbox{\boldmath{x}}\in\mathbb{C}^{N} and \mbox{\boldmath{y}}\in\mathbb{C}^{N} such that $|x_{i}|\leq|y_{i}|,i=1,\ldots,N$, {g(\mbox{\boldmath{x}})}\leq{g(\mbox{\boldmath{y}})}. Additionally, for any fixed $u=\bar{u}$ the constraints on the variables $\lambda_{i}$, $i=1,\ldots,N$ are not coupled. Thus, an optimal solution to

[TABLE]

can be found solving the following $N$ scalar optimization problems

[TABLE]

The closed form solution to ${\cal{P}}_{3}^{i}(\bar{u})$ can be obtained analyzing

[TABLE]

where the variables $x,y,a,b$ are given by

[TABLE]

Since $|x-y|$ is a monotonically decreasing function if $x\leq y$ and a monotonically increasing function if $x\geq y$, it follows that the optimal solution is $x_{min}=y$ if $a\leq y\leq b$. If $y\leq a\leq b$ $x_{min}=a$. Thus, $x_{min}=\max(y,a)$ as long as $y\leq b$. Moreover, $x_{min}=b$ if $y\geq b$ implying that the optimal solution to $\tilde{{\cal{P}}}_{3}$ can be written in closed form as

[TABLE]

Replacing (31) in (32), the optimal solution to ${\cal{P}}_{3}^{i}(\bar{u})$ is

[TABLE]

∎

-C Proof of Theorem III.3

Proof:

Let \mbox{\boldmath{\Lambda}}^{\star}(\bar{u}) be an optimal solution to problem ${\cal{P}}_{2}^{\prime}(\bar{u})$ in (11) as provided by Lemma III.2. Hence, a minimizer of ${\cal{P}}_{2}$ is \mbox{\boldmath{\Lambda}}^{\star}=\mbox{\boldmath{\Lambda}}^{\star}(u_{1}^{\star}) with $u_{1}^{\star}$ an optimal solution to the one-dimensional problem

[TABLE]

Based on [39, Property 7.4.21] ${\cal{P}}_{3}^{\prime}$ can be reformulated as

[TABLE]

where $g(\cdot)$ is the gauge function associated with the considered unitary invariant norm and $\lambda_{i}^{\star}(u)$, $i=1,\ldots,N$, are the entries of the vectorial function \mbox{\boldmath{\lambda}}^{\star}(u) given in (13). Now, notice that, $g\left(|\lambda_{1}^{\star}(u)-d_{1}|,\ldots,|\lambda_{N}^{\star}(u)-d_{N}|\right)$ is a continuous function since obtained as the composition of continuous functions. Besides, due to [48, p. 88], the objective function in ${\cal{P}}_{3}$ is convex. Hence, an optimal solution to ${\cal{P}}_{1}$ just requires the solution of the convex optimization problem ${\cal{P}}_{3}$. Additionally, the set of the optimal solutions to ${\cal{P}}_{3}$ defines a closed convex interval bounded below ensuring the existence of the lowest optimal solution $u^{\star}$. ∎

-D Proof of Theorem III.4

Proof:

The proof is organized in four different parts accordingly to the provided claims.

As proved in the supplementary material, $G_{1}(u)$ is a convex and differentiable function within the interval $u\in[1,d_{1}]$. Recalling that a differentiable function is convex if and only if its derivative is an increasing function, it can be claimed that $\frac{dG_{1}(u)}{du}$ is increasing within $u\in[1,d_{1}]$. Since \frac{dG_{1}(u)}{du}\Big{|}_{u=1}\geq 0, then

[TABLE]

namely $G_{1}(u)$ is an increasing function in $u\in[1,d_{1}]$. Hence, the minimum to ${\cal{P}}_{4}$ is attained in $u^{\star}=1$. 2. 2.

Assume that

[TABLE]

Since, $G_{1}(u)$ is a convex and differentiable function within the interval $u\in[1,d_{1}]$,

[TABLE]

therefore $G_{1}(u)$ is decreasing in $[1,d_{1}]$ and $u_{1}^{\star}=d_{1}$ is the optimal solution as long as111111Without loss of generality, it is assumed that $v_{i}\neq v_{j}$ and $v_{i}\neq\kappa_{M}v_{j}$ for all $i\neq j$ with $1\leq i,j\leq\bar{N}+1$. $d_{1}>d_{N}$. Indeed, if there exists $u^{\star}<d_{1}$, $G_{1}(u)$ would be constant over $[u^{\star},u_{1}^{\star}]$ but this is not possible since $d_{1}>d_{N}$. 3. 3.

If $\frac{d_{1}}{\kappa_{M}}\leq d_{N}$, $G_{1}(u)$ is strictly decreasing up to $\frac{d_{1}}{\kappa_{M}}$ and increasing after this point implying that $u^{\star}=\frac{d_{1}}{\kappa_{M}}$. 4. 4.

Finally, if $\frac{d_{1}}{\kappa_{M}}>d_{N}$, the optimal solution is unique. Indeed, if $G_{1}(u)$ is constant over an interval its second order derivative is null in this set. Now, observe that $G_{1}(u)$ is composed of convex functions whose second order derivative exists except for a finite number of points. Additionally, in any regular interval there is at least one of such functions that is strictly convex due to the assumption $\frac{d_{1}}{\kappa_{M}}>d_{N}$, and hence the initial claim is contradicted. Now, owing to

[TABLE]

the optimal point belongs to $]1,d_{1}[$. Since $G_{1}(u)$ is a convex and differentiable function, $u^{\star}$ is the optimal point if and only if

[TABLE]

Let us now fully characterize condition (36). To this end, let $\alpha\in\{1,2,\ldots,\bar{N},\bar{N}+1\}$ be the smallest index such that $v_{\alpha}<u^{\star}$ ($\alpha\geq 2$ since $u^{\star}<d_{1}=v_{1}$; moreover $\alpha\leq\bar{N}+1$ because $u^{\star}>1=v_{\bar{N}+1}$). Besides, let $\beta\in\{1,2,\ldots,\bar{N},\bar{N}+1\}$ be the largest index such that $\frac{v_{\beta}}{\kappa_{M}}>u^{\star}$ ($\beta\geq 1$ because121212\frac{dG_{1}(u)}{du}\Big{|}_{u=\frac{d_{1}}{\kappa_{M}}}>0 since each term composing $G_{1}(u)$ is an increasing function over $u\geq\frac{d_{1}}{\kappa_{M}}$ and at least the term associated with $d_{N}$ has a strictly positive derivative. $u^{\star}<\frac{d_{1}}{\kappa_{M}}$ and $\beta\leq\bar{N}$ since $u^{\star}>1>\frac{1}{\kappa_{M}}=\frac{v_{\bar{N}+1}}{\kappa_{M}}$). Notice that, $\alpha>\beta$, otherwise $v_{\alpha}\geq v_{\beta}>\frac{v_{\beta}}{\kappa_{M}}>u^{\star}$, which contradicts $u^{\star}>v_{\alpha}$. Therefore, it will exist a neighborhood $B_{u{{}^{\star}}}$ of $u^{\star}$, contained within the interval $[v_{\alpha},\frac{v_{\beta}}{\kappa_{M}}]$, such that $G_{1}(u)$ can be expressed as

[TABLE]

Computing the derivative of (37) with respect to $u$, and imposing the optimality condition (36), the minimizer $u^{\star}$ is given by

[TABLE]

On the other hand, let $\bar{u}$ a value such that

[TABLE]

with $\bar{\alpha}\in\{1,2,\ldots,\bar{N},\bar{N}+1\}$ the smallest indexes such that $v_{\bar{\alpha}}<\bar{u}$ and $\bar{\beta}\in\{1,2,\ldots,\bar{N},\bar{N}+1\}$ the largest indexes such that $\frac{v_{\bar{\beta}}}{\kappa_{M}}>\bar{u}$. This means that it exists a neighborhood $B_{\bar{u}}$ of $\bar{u}$ contained within the interval $\left[v_{\bar{\alpha}},\frac{v_{\bar{\beta}}}{\kappa_{M}}\right]$, such that $\forall u\in B_{\bar{u}}$, $G_{1}(u)$ is given by

[TABLE]

Now, computing the derivative of (40), condition (39) implies that \frac{dG_{1}(u)}{du}\Big{|}_{u=\bar{u}}=0, i.e., $\bar{u}=u^{\star}$.

∎

-E Proof of Theorem III.5

Proof:

The proof is organized in five different parts accordingly to the provided claims.

If $d_{i}\leq 1$, $i=1,\ldots,N$, the functions $h_{i}(u)$, $i=1,\ldots,N$, in (17) are monotonically increasing over $u\geq\frac{1}{\kappa_{M}}$. Hence, $G_{2}(u)$ is monotonically increasing implying that $u^{\star}=\frac{1}{\kappa_{M}}$. 2. 2.

Let $\mathcal{I}_{q}=\{i\,:\,\,d_{i}\leq 1\}$ and $\mathcal{I}_{p}=\{i\,:\,\,d_{i}>1\}$; by assumption, $1<d_{1}\leq\kappa_{M}$ and $\mathcal{I}_{q}\neq\emptyset$. Now, observe that

[TABLE]

where

[TABLE]

To proceed further, notice that

[TABLE]

where $p_{max}=\max\{j\in\mathcal{I}_{p}\}$. In fact,

[TABLE]

implying that

[TABLE]

moreover, $\forall j\in\mathcal{I}_{p}$, $\frac{d_{j}}{\kappa_{M}}<1$ and

[TABLE]

thus $h^{p}(u)=h_{p_{max}}(u)$ over $u\geq 1$. Finally, since $h^{q}(u)=h_{N}(u)$ if $u\geq\frac{1}{\kappa_{M}}$,

[TABLE]

To proceed further, observe that

[TABLE]

since $u-d_{N}\geq\max\left\{u-d_{p_{max}},0\right\}$, $u\geq 1$. Now, let $\eta=\frac{d_{1}+d_{N}-1}{\kappa_{M}}$ be the point such that $1-d_{N}=d_{1}-\eta\kappa_{M}$, where $\eta\leq\frac{d_{1}}{\kappa_{M}}\leq 1$ since $d_{N}\leq 1$. If $\eta\leq\frac{1}{\kappa_{M}}$, $G_{2}(u)=h_{N}(u)$, $u\geq\frac{1}{\kappa_{M}}$, implying that $u^{\star}=\frac{1}{\kappa_{M}}$. Otherwise,

[TABLE]

and $u^{\star}=\eta$ since $G_{2}(u)$ is a strictly decreasing function up to $\eta$ and it monotonically increases if $u\geq\eta$. Thus, $u^{\star}=\max\left\{\eta,\frac{1}{\kappa_{M}}\right\}$. 3. 3.

Assume $1<d_{1}\leq\kappa_{M}$ and $d_{N}>1$. Since, $d_{i}>1$, $i=1,\ldots,N$, $G_{2}(u)=h_{1}(u)$ as $\frac{1}{\kappa_{M}}\leq u\leq d_{N}$; hence, $G_{2}(u)$, is a strictly decreasing function over $\frac{1}{\kappa_{M}}\leq u\leq\frac{d_{1}}{\kappa_{M}}$ and $G_{2}(\frac{d_{1}}{\kappa_{M}})=0$ since $\frac{d_{1}}{\kappa_{M}}\leq 1<d_{N}$. As a consequence $u^{\star}=\frac{d_{1}}{\kappa_{M}}$ and \mbox{\boldmath{X}}^{\star}=\mbox{\boldmath{S}}. 4. 4.

Consider $d_{1}>\kappa_{M}$ and $d_{N}\leq 1$. In this case,

[TABLE]

where, as in item 2), $p_{max}$ is the index of the lowest eigenvalue higher than 1. To proceed further, let $\eta=\frac{d_{1}+d_{N}-1}{\kappa_{M}}$; if $\eta\leq 1$,

[TABLE]

Furthermore

[TABLE]

Hence, based on (42), (43), and (44)

[TABLE]

implying that $u^{\star}=\max\left\{\frac{1}{\kappa_{M}},\eta_{1}\right\}$.

Now, assume $\eta>1$ and let $\eta_{1}=\frac{d_{1}+d_{N}}{1+\kappa_{M}}$, where $1\leq\eta_{1}\leq\frac{d_{1}}{\kappa_{M}}$. According to (41), $G_{2}(u)=h_{1}(u)$, if $\frac{1}{\kappa_{M}}\leq u\leq 1$. Additionally,

[TABLE]

Besides,

• •

if $\frac{d_{p_{max}}}{\kappa_{M}}\leq 1$

[TABLE]

• •

otherwise, $\frac{d_{p_{max}}}{\kappa_{M}}>1$ and

[TABLE]

Summarizing,

[TABLE]

As a consequence, $G_{2}(u)$ is a strictly decreasing function over $\frac{1}{\kappa_{M}}\leq u\leq\eta_{1}$ and monotonically increases over $\eta_{1}\leq u\leq d_{1}$. Thus, $u^{\star}=\eta_{1}$. 5. 5.

Let $d_{1}>\kappa_{M}$ and $d_{N}>1$.

• •

If $d_{N}\leq{\frac{d_{1}}{\kappa_{M}}}$, let $\eta_{1}=\frac{d_{1}+d_{N}}{1+\kappa_{M}}\leq{\frac{d_{1}}{\kappa_{M}}}$. Hence, $G_{2}(u)=h_{1}(u)$, if $\frac{1}{\kappa_{M}}\leq u\leq\eta_{1}$. In fact, within this interval, $\forall i\in\{1,\ldots,N\}$

[TABLE]

Moreover, $\eta_{1}\geq d_{N}$, since $\frac{d_{1}+d_{N}}{1+\kappa_{M}}-d_{N}=\frac{d_{1}/\kappa_{M}-d_{N}}{\kappa_{M}(1+\kappa_{M})}\geq 0$ and consequently $G_{2}(u)=h_{N}(u)$ if $\eta_{1}\leq u\leq d_{1}$. As a result, $u^{\star}=\frac{d_{1}+d_{N}}{1+\kappa_{M}}$;

• •

otherwise, $d_{N}>\frac{d_{1}}{\kappa_{M}}$ and following the same line of reasoning as in item 3) it follows that $u^{\star}=\frac{d_{1}}{\kappa_{M}}$ as well as \mbox{\boldmath{X}}^{\star}=\mbox{\boldmath{S}}.

∎

SUPPLEMENTARY MATERIAL ORGANIZATION

The following additional appendices contain supplementary material for this paper. More specifically, the following sections provide detailed proofs of some claims within the manuscript.

APPENDIX S1. Suffiecient condition for estimator measurability

Proof:

Let us show that the estimator \widehat{\mbox{\boldmath{M}}} is a continuous (thus measurable) function of \widehat{\mbox{\boldmath{S}}} if $\|\cdot\|$ is strictly convex. To this end, let \widehat{\mbox{\boldmath{S}}}_{i}\in\mathbb{H}^{N}, $i=1,2\ldots$ be a sequence of positive semi-definite matrices converging to \widehat{\mbox{\boldmath{S}}}^{\star} as $i\rightarrow\infty$ (e.g., \widehat{\mbox{\boldmath{S}}}_{i}\rightarrow\widehat{\mbox{\boldmath{S}}}^{\star}) and \widehat{\mbox{\boldmath{M}}}_{i} the resulting sequence of covariance estimates. Now, the goal is to show that \widehat{\mbox{\boldmath{M}}}_{i}\rightarrow\widehat{\mbox{\boldmath{M}}}^{\star} with \widehat{\mbox{\boldmath{M}}}^{\star} the estimate associated with \widehat{\mbox{\boldmath{S}}}^{\star}. Based on [49, Lemma IV.1], \displaystyle{\lim_{i\rightarrow\infty}}\displaystyle{\min_{\widehat{\mbox{\boldmath{M}}}\in\mathcal{M}}}\|\widehat{\mbox{\boldmath{S}}}_{i}-\widehat{\mbox{\boldmath{M}}}\|=\displaystyle{\min_{\widehat{\mbox{\boldmath{M}}}\in\mathcal{M}}}\|\widehat{\mbox{\boldmath{S}}}^{\star}-\widehat{\mbox{\boldmath{M}}}\|, where $\mathcal{M}$ is a compact set contained in (3) that encompasses the optimal solutions associated with \widehat{\mbox{\boldmath{S}}}^{\star} and all the \widehat{\mbox{\boldmath{S}}}_{i}. Now, let \left(\widehat{\mbox{\boldmath{S}}}_{i^{\prime}},\widehat{\mbox{\boldmath{M}}}_{i^{\prime}}\right) be a sequence extracted from \left(\widehat{\mbox{\boldmath{S}}}_{i},\widehat{\mbox{\boldmath{M}}}_{i}\right) that converges to a point \left(\widehat{\mbox{\boldmath{S}}}^{\star},\widehat{\mbox{\boldmath{M}}}_{1}^{\star}\right). Due to the continuity of $\|\cdot\|$, \|\widehat{\mbox{\boldmath{S}}}^{\star}-\widehat{\mbox{\boldmath{M}}}_{1}^{\star}\|=\displaystyle{\lim_{i^{\prime}\rightarrow\infty}}\displaystyle{\min_{\widehat{\mbox{\boldmath{M}}}\in\mathcal{M}}}\|\widehat{\mbox{\boldmath{S}}}_{i^{\prime}}-\widehat{\mbox{\boldmath{M}}}\|=\displaystyle{\min_{\widehat{\mbox{\boldmath{M}}}\in\mathcal{M}}}\|\widehat{\mbox{\boldmath{S}}}^{\star}-\widehat{\mbox{\boldmath{M}}}\|, namely \widehat{\mbox{\boldmath{M}}}_{1}^{\star}=\widehat{\mbox{\boldmath{M}}}^{\star} is the optimal solution. Finally, if \widehat{\mbox{\boldmath{M}}}_{i} does not converge to \widehat{\mbox{\boldmath{M}}}^{\star} there exists at least one extract sequence \widehat{\mbox{\boldmath{M}}}_{i_{1}^{\prime}} such that \widehat{\mbox{\boldmath{M}}}_{i_{1}^{\prime}}\rightarrow\widehat{\mbox{\boldmath{M}}}_{2}^{\star}\neq\widehat{\mbox{\boldmath{M}}}^{\star}, which is an absurd since the optimal solution is unique. ∎

APPENDIX S2. Proof of the equivalence between problems ${\cal{P}}$ and ${\cal{P}}_{1}$

Proof:

Let \bar{\mbox{\boldmath{M}}} be a feasible solution to ${\cal{P}}$ and \left(\bar{\mbox{\boldmath{R}}},\bar{\sigma}_{n}^{2}\right) be such that \left(\bar{\mbox{\boldmath{M}}},\bar{\mbox{\boldmath{R}}},\bar{\sigma}_{n}^{2}\right) satisfies the constraints in (3). Since \left(\bar{\mbox{\boldmath{M}}},\bar{\mbox{\boldmath{R}}}+(\bar{\sigma}_{n}^{2}-{\sigma^{2}})\mbox{\boldmath{I}},{\sigma^{2}}\right) is also feasible to (3), ${\cal{P}}$ is equivalent to

[TABLE]

Next, observe that the set

[TABLE]

can be recast as

[TABLE]

Hence, defining \mbox{\boldmath{X}}=\frac{\mbox{\boldmath{M}}}{\sigma^{2}}, ${\cal{P}^{\prime}}$ boils down to

[TABLE]

where \mbox{\boldmath{S}}=\frac{\widehat{\mbox{\boldmath{S}}}}{\sigma^{2}}. Since \frac{\lambda_{max}(\mbox{\boldmath{X}})}{\lambda_{min}(\mbox{\boldmath{X}})} is quasi-convex function over the set \mbox{\boldmath{X}}\succeq\mbox{\boldmath{I}}, ${\cal{P}}_{1}$ is a convex problem. Finally, notice that Problem ${\cal{P}}_{1}$ is equivalent to

[TABLE]

Now, since the objective in ${\cal{P}}_{1}^{\prime}$ is a continuous function and the feasible set is a compact set, Weierstrass theorem ensures the existence of a feasible point \mbox{\boldmath{X}}^{\star} to ${\cal{P}}_{1}^{\prime}$ such that v({\cal{P}}_{1})=v({\cal{P}}_{1}^{\prime})=\|\mbox{\boldmath{X}}^{\star}-\mbox{\boldmath{S}}\|, which concludes the proof. ∎

APPENDIX S3. Proof of the equivalence between ${\cal{P}}_{2}$ and ${\cal{P}}_{2}^{\prime}(u)$

Proof:

Let us observe that the constraint set

[TABLE]

is equivalent to

[TABLE]

Now, introducing an auxiliary variable $u>0$, the set (52) can be cast as

[TABLE]

Indeed, if \bar{\mbox{\boldmath{\lambda}}} is a feasible point to (52), \left(\bar{\mbox{\boldmath{\lambda}}},\min\left\{\bar{\mbox{\boldmath{\lambda}}}\right\}\right) is a feasible point to (53). On the other hand, if \left(\bar{\mbox{\boldmath{\lambda}}}^{1},\bar{u}^{1}\right) is a feasible point to (53),

[TABLE]

As a result, \bar{\mbox{\boldmath{\lambda}}}_{1} is a feasible point to (52). Finally, since (53) is empty if $u<\frac{1}{\kappa_{M}}$, ${\cal{P}}_{2}^{\prime}(u)$ follows. ∎

APPENDIX S4. Proof of the results in equations (16) and (17)

Proof:

To prove equations (16) and (17), let us introduce the following notation

[TABLE]

Hence, $\lambda_{i}(u)=\min(\kappa_{M}u,\beta_{1})$. Now, assuming $d_{i}>1$:

• •

if $\frac{1}{\kappa_{M}}\leq u<\frac{d_{i}}{\kappa_{M}}$, then $u<d_{i}$ and $\beta_{1}=d_{i}$. Moreover, $\quad u\kappa_{M}<d_{i}$, implying that $\lambda_{i}(u)=u\kappa_{M}$ as well as $h_{i}(u)=|u\kappa_{M}-d_{i}|=d_{i}-u\kappa_{M}$;

• •

if $\frac{d_{i}}{\kappa_{M}}\leq u<d_{i}$, as in the previous item $\beta_{1}=d_{i}$. Now, ${u\kappa_{M}}\geq{d_{i}}$, and thus $\lambda_{i}(u)=d_{i}$, i.e., $h_{i}(u)=0$;

• •

if $u\geq d_{i}$, then $\alpha_{1}=u$ and $\beta_{1}=u$; thus, $\lambda_{i}(u)=u$ and $h_{i}(u)=u-d_{i}$.

Instead, when $d_{i}\leq 1$:

• •

if $\frac{1}{\kappa_{M}}\leq u<1$, then $\alpha_{1}=1$ and $\beta_{1}=1$; hence, $\lambda_{i}(u)=1$ and $h_{i}(u)=1-d_{i}$;

• •

if $u\geq 1$, then $\beta_{1}=u$; thus, $\lambda_{i}(u)=u$ and $h_{i}(u)=u-d_{i}$.

∎

APPENDIX S5. Proof of the claims concerning $G_{1}(u)$

Proof:

Let us start assuming $d_{1}\leq 1$. This implies that $d_{i}\leq 1$, $i=1,\ldots,N$ and

[TABLE]

Consequently, for $i=1,\ldots,N$, each $G_{1}(u)$ is a constant function over $u\in\left[\frac{1}{\kappa_{M}},1\right[$ and it is monotonically increasing in $u\in\left[1,+\infty\right]$. Thus, $u^{\star}=\frac{1}{\kappa_{M}}$.

Consider now $1<d_{1}\leq\kappa_{M}$. Let $\mathcal{I}_{p}=\{i:d_{i}>1\}$ be the set of indexes corresponding to the eigenvalues greater than $1$. Since $\frac{d_{1}}{\kappa_{M}}\leq 1<d_{p_{max}}$, where $p_{max}\in\mathcal{I}_{p}$ is the index associated with the lowest eigenvalue greater than 1,

[TABLE]

is monotonically decreasing over $\frac{1}{\kappa_{M}}\leq u\leq\frac{d_{1}}{\kappa_{M}}$. Moreover, (55) is constant over $\frac{d_{1}}{\kappa_{M}}\leq u\leq d_{p_{max}}$ and monotonically increases for $u\geq d_{p_{max}}$. Furthermore, if $d_{i}\leq 1$, $(h_{i}(u))^{2}$ is a constant function over $\frac{1}{\kappa_{M}}\leq u<1$ and monotonically increases over $u\geq 1$. Hence, $u^{\star}=\frac{d_{1}}{\kappa_{M}}$.

Finally, if $d_{1}>\kappa_{M}$ $G_{1}(u)$ is a strictly decreasing function over $\frac{1}{\kappa_{M}}\leq u<1$, whereas it is a strictly increasing function within $d_{1}\leq u<+\infty$. Hence, any minimum belongs to the interval $\left[1,d_{1}\right]$.

Let us now focus on the differentiability of $G_{1}(u)$. If $d_{i}\leq 1$, the derivative of $(h_{i}(u))^{2}$ is

[TABLE]

Hence, $\frac{d(h_{i}(u))^{2}}{du}$ is continuous over $u\geq 1$ where in $u=1$ the right derivative is considered. If $d_{i}>1$,

[TABLE]

Hence, in each of the sub-intervals $\left[\frac{1}{\kappa_{M}},\frac{d_{i}}{\kappa_{M}}\right[$, $\left]\frac{d_{i}}{\kappa_{M}},d_{i}\right[$, and $\left]d_{i},+\infty\right[$, $\frac{d(h_{i}(u))^{2}}{du}$ is continuous; moreover, in correspondence of the points $u=\frac{d_{i}}{\kappa_{M}}$ and $u=d_{i}$, the right and left derivatives coincide implying that the overall derivative is continuous over $u\geq\frac{1}{\kappa_{M}}$. Since $\frac{1}{\kappa_{M}}\leq 1$, the function $\frac{dG_{1}(u)}{du}=\displaystyle{\sum_{i=1}^{N}}\frac{d(h_{i}(u))^{2}}{du}$, is continuous over $u\in\left[1,d_{1}\right]$.

As to the convexity of $G_{1}(u)$, it easily follows from the convexity of each term. ∎