# New Methods for MLE of Toeplitz Structured Covariance Matrices with
Applications to RADAR Problems

Augusto Aubry, Prabhu Babu, Antonio De Maio, and Rikhabchand Jyothi Augusto
Aubry and Antonio De Maio are with the Department of Electrical Engineering
and Information Technology, Universita degli Studi di Napoli “Federico II”,
DIETI, Via Claudio 21, I-80125 Napoli, Italy (E-mail<EMAIL_ADDRESS>ademaio@unina.it). Prabhu Babu and Rikhabchand Jyothi are with CARE, IIT
Delhi, New Delhi, 110016, India (E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This work considers Maximum Likelihood Estimation (MLE) of a Toeplitz
structured covariance matrix. In this regard, an equivalent reformulation of
the MLE problem is introduced and two iterative algorithms are proposed for
the optimization of the equivalent problem. Both the strategies are based on
the Majorization Minimization (MM) framework and hence enjoy nice properties
such as monotonicity and ensured convergence to a stationary point of the
equivalent MLE problem. The proposed algorithms are also extended to deal with
MLE of other related covariance structures, namely, the banded Toeplitz,
Toeplitz-block-Toeplitz, low rank Toeplitz structure plus a scalar matrix
(accounting for white noise), and finally Toeplitz matrices satisfying a
condition number constraint. Through numerical simulations, it is shown that
new methods provide satisfactory performance levels in terms of both mean
square estimation error and signal-to-interference-plus-noise ratio.

###### Index Terms:

Toeplitz covariance matrix, Maximum likelihood estimation, Banded Toeplitz,
Toeplitz-block-Toeplitz, Adaptive radar signal processing, Array processing,
Spectral estimation

##  I Introduction

Estimation of the interference covariance matrix has diverse applications in
radar signal processing such as direction of arrival estimation, target
detection, adaptive beamforming, and sidelobe canceller design ([1, 2, 3, 4]).
In these situations, the interference covariance matrix is estimated from the
secondary/training data, which are assumed target-free and collected from
spatial and/or temporal returns corresponding to range cells close to the one
of interest. The Sample Covariance Matrix (SCM), which is the unstructured
Maximum Likelihood Estimate (MLE) when the data follows a complex circular
Gaussian distribution, does not always represent a good choice for the
covariance estimate especially in the presence of a small number of training
data and/or when mismatches in training data spectral properties occur [5, 6].
A well-known strategy often discussed in the open literature to improve the
performance of a covariance estimator relies on the incorporation of some _a
priori_ knowledge about its structure. For instance, in some radar
applications, it is customary to suppose that data come from a stationary
Gaussian random process, leading to a Toeplitz Structured Covariance (TSC)
matrix. Leveraging this information, one can obtain (under the design
conditions) a more reliable estimator as compared with the SCM [7]. Other than
radar applications, the estimation of a TSC matrix is encountered in speech
recognition [8], spectral estimation [2], and in hyperspectral imaging [9].
So far several algorithms have been proposed for estimating a TSC matrix. Let
us first discuss those for MLE. According to the Caratheodory parametrization
([2, 10, 11]) a Toeplitz covariance matrix ${\mathbf{T}}\in\mathbb{H}^{m\times
m}$ can always be decomposed as

$\begin{array}[]{ll}{\mathbf{T}}={{\mathbf{A}}}\tilde{{\mathbf{P}}}{{\mathbf{A}}}^{H};\>\>[\tilde{{\mathbf{P}}}]_{k,k}\geq
0\end{array}$ (1)

where

${\mathbf{A}}=\begin{bmatrix}1&1&\cdots&1\\\
e^{j\omega_{1}}&e^{j\omega_{2}}&\cdots&e^{j\omega_{r}}\\\
\vdots&\vdots&\ddots&\vdots\\\
e^{j(m-1)\omega_{1}}&e^{j(m-1)\omega_{2}}&\cdots&e^{j(m-1)\omega_{r}}\end{bmatrix}$
(2) $\tilde{{\mathbf{P}}}=\begin{bmatrix}\tilde{p}_{1}&0&\dots&0\\\
0&\tilde{p}_{2}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\dots&\tilde{p}_{r}\end{bmatrix}$ (3)

where $\omega_{i}$ and $\tilde{p}_{i}$, $i=1,2,\cdots,r\leq m$, denote some
angular frequencies and their corresponding powers and $r$ indicates the rank
of ${\mathbf{T}}$. Capitalizing on this parametrization, Circulant Embedding
(CE) of Toeplitz matrix ([12, 13, 14]) can be used to compute approximately
the MLE of ${\mathbf{T}}$. According to CE, a positive semidefinite $m\times
m$ Toeplitz matrix is modeled as:

$\begin{array}[]{ll}{\mathbf{T}}=\tilde{{\mathbf{F}}}{\mathbf{P}}\tilde{{\mathbf{F}}}^{H};\>\>{\mathbf{P}}={\textrm{diag}}([p_{1},p_{2},\cdots,p_{L}]),p_{k}\geq
0\end{array}$ (4)

where $\tilde{{\mathbf{F}}}=[{\mathbf{I}}_{m\times
m}\>{\boldsymbol{0}}_{m\times L-m}]{\mathbf{F}}$, ${\mathbf{I}}_{m\times m}$
is the identity matrix of size $m\times m$, ${\boldsymbol{0}}_{m\times L-m}$
is the zero matrix of size $m\times L-m$, ${\mathbf{F}}$ is the normalized DFT
matrix of size $L\geq 2m-1$ and ${\mathbf{P}}$ is a diagonal matrix of size
$L\times L$ with diagonal elements $p_{k}\geq 0$. Therefore, the matrix
${\mathbf{T}}$ is completely parametrized by the diagonal matrix
${\mathbf{P}}$. Although estimating the Toeplitz covariance matrix using CE
seems attractive, the representation in (4) is valid only for a subset of
Toeplitz covariance matrices. This can be intuitively justified because the
Caratheodory parametrization in (1) does not give restrictions on the
frequencies spacing, while the CE in (4) strictly requires the frequencies to
lie on the Fourier grid. Hence, for some Toeplitz matrices, the
parametrization in (4) is only approximate. Based on CE, [15] and [16]
proposed an iterative algorithm based on Expectation-Maximization (EM) for MLE
of ${\mathbf{T}}$. By modifying the M step in the EM procedure, in [17] the
technique has been extended to deal with the banded Toeplitz covariance case.
In [18], still leveraging CE framework, a Majorization Minimization (MM) based
optimization, with faster convergence than the EM in [15] and [16], has been
introduced. In [19] a closed-form estimator is designed by invoking the
extended invariance principle to deal with the Toeplitz constraint. Finally,
in [20], an efficient approximation of a Toeplitz covariance matrix under a
rank constraint has been proposed forcing the eigenvectors to be the same as
those of the SCM whereas the Toeplitz constraint has been explicitly imposed
while estimating the eigenvalues. Other than the MLE, several other
alternative paradigms have been considered for the problem at hand. Recently,
in [21] the Toeplitz structure is forced together with a condition number
constraint via SCM projection onto a suitable constraint set. Other geometric
based approaches for the TSC estimate have also been proposed in [22, 23].
In this manuscript, two iterative algorithms denoted as Alternating Projection
Based TOeplitz Covariance Matrix Estimation 1 (ATOM1) and ATOM2 are devised
leveraging a suitable reformulation of the MLE problem and the MM framework.
Both ATOM1 and ATOM2 involve the construction of a surrogate function together
with its optimization. Specifically, the two procedures construct distinct
surrogate functions and therefore solve different surrogate minimization
problems. While ATOM1 addresses the surrogate minimization problem using the
Alternating Direction Method of Multipliers (ADMM), ATOM2 handles it via
alternating projection or Dykstra’s algorithm. However, both the procedures
directly estimate the Toeplitz covariance matrix without reparametrization via
the CE. ATOM2 is also extended to include other constraints, such as an upper
bound to the condition number and a low rank plus white noise covariance
structure. Finally, ATOM2 is modified to deal with banded Toeplitz and
Toeplitz-block-Toeplitz structures. The major contributions of this paper can
be summarized as follows:

  1. 1.

Two iterative algorithms ATOM1 and ATOM2 are proposed based on the MM
framework to address MLE of a Toeplitz covariance matrix. Their computational
complexities are thoroughly discussed. Also, the convergence of the procedures
to a stationary point of the MLE problem is established.

  2. 2.

The extensions of ATOM2 to handle additional covariance structures, such as
banded Toeplitz, Toeplitz-block-Toeplitz, and low rank Toeplitz term plus a
scalar matrix are discussed, together with the possibility to account for a
condition number constraint.

  3. 3.

The Cramer-Rao lower bounds for the MLE estimation of Toeplitz, banded
Toeplitz, and Toeplitz-block-Toeplitz covariance matrices are derived.

  4. 4.

The proposed algorithms and their extensions are compared with some state-of-
the-art procedures via numerical simulations using the Mean Square Error (MSE)
and Signal-to-Interference-plus-Noise Ratio (SINR) (for case studies related
to radar application) as performance metrics.

The organization of the paper is as follows. The MLE problem of Toeplitz
covariance matrix for complex circular Gaussian observations is formulated in
Section II. In Section III, ATOM1 and ATOM2 algorithms are proposed, along
with a discussion on their computational complexity and implementation
aspects. Also, their convergence properties are studied. At the end of the
section, the extension of ATOM2 to handle additional constraints along with
the Toeplitz requirement is discussed too. In Section IV, the Cramer-Rao lower
bounds for the MLE of Toeplitz, banded Toeplitz, and Toeplitz-block-Toeplitz
covariance matrices are computed. In Section V, the proposed algorithms are
compared with some state-of-the-art algorithms, and finally, concluding
remarks are given in Section VI.
A. Notation
Throughout the paper, bold capital and bold small letter denote matrix and
vector, respectively. A scalar is represented by a small letter. The value
taken by an optimization vector ${{\boldsymbol{x}}}$ at the ${t^{th}}$
iteration is denoted by ${{\boldsymbol{x}}_{t}}$. The notations $\mathbb{R}$,
$\mathbb{R}^{m\times 1}$, $\mathbb{C}^{m\times 1}$, $\mathbb{R}^{m\times m}$,
$\mathbb{C}^{m\times m}$, and $\mathbb{H}^{m\times m}$ are used to represent
the sets of real numbers, $m$ dimensional vectors of real numbers, $m$
dimensional vectors of complex numbers, $m\times m$ matrices of real numbers,
$m\times m$ matrices of complex numbers, and $m\times m$ Hermitian matrices,
respectively. Superscripts $(\cdot)^{T}$, $(\cdot)^{*}$, $(\cdot)^{H}$, and
$(\cdot)^{-1}$ indicate the transpose, complex conjugate, complex conjugate
transpose, and inverse, respectively. The trace and the determinant of a
matrix ${\boldsymbol{X}}$ are denoted by ${\rm{Tr}}({\boldsymbol{X}})$ and
$|{\boldsymbol{X}}|$, respectively. The notation $[{{\boldsymbol{X}}}]_{i}$ is
used to represent the $i^{th}$ column of the matrix ${\boldsymbol{X}}$. The
$\otimes$ symbol indicates the Kronecker product while the gradient of a
function $f$ is denoted by $\nabla f$. The $\succeq$ symbol (and its strict
form $\succ$) is used to denote the generalized matrix inequality: for any
${\boldsymbol{X}}\in\mathbb{H}^{m\times m}$, ${\boldsymbol{X}}\succeq 0$ means
that ${\boldsymbol{X}}$ is a positive semi-definite matrix
(${\boldsymbol{X}}\succ 0$ for positive definiteness). Besides, for any
${\boldsymbol{X}}\in\mathbb{H}^{m\times m}$, ${\rm{eig}}({\boldsymbol{X}})$ is
the vector collecting the eigenvalues of ${\boldsymbol{X}}$ with the maximum
and the minimum eigenvalue indicated as $\lambda_{max}({\boldsymbol{X}})$ and
$\lambda_{min}({\boldsymbol{X}})$, respectively. The Euclidean norm of the
vector ${\boldsymbol{x}}$ is denoted by ${\|{\boldsymbol{x}}\|}_{2}$,
$|{\boldsymbol{x}}|$ indicates the element wise modulus of the vector
${\boldsymbol{x}}$. The notation $\textbf{E}[\cdot]$ stand for statistical
expectation. Finally, for any
${\boldsymbol{X}},{\boldsymbol{Y}}\in\mathbb{R}^{m\times m}$,
$\rm{max}({\boldsymbol{X}},{\boldsymbol{Y}})$ refers to the matrix containing
the element wise maximum between ${\boldsymbol{X}}$ and ${\boldsymbol{Y}}$.

##  II Problem Formulation

Let us assume the availability of $n$ independent and identically distributed
vectors $\\{{\mathbf{y}}_{1},{\mathbf{y}}_{2},\cdots,{\mathbf{y}}_{n}\\}$,
where each ${\mathbf{y}}_{i}$ is of size $m$ and follows a $m$-variate complex
circular Gaussian distribution with zero mean and covariance matrix
${\mathbf{R}}$. The maximum likelihood covariance estimation problem can be
formulated as

$\begin{array}[]{ll}\underset{{\mathbf{R}}\succ 0}{\rm
minimize}\>f({\mathbf{R}}){=}\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{\mathbf{R}}^{-1}{\mathbf{y}}_{i}+\log|{\mathbf{R}}|\end{array}$
(5)

If $n\geq m$, problem (5) has a unique minimizer with probability one which is
given by the SCM, i.e.,
${\mathbf{R}}_{\textrm{SCM}}=\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}{\mathbf{y}}_{i}^{H}$.
However, if the random process, from which each observation is drawn, is
stationary (at least in wide sense) then the covariance matrix also exhibits a
Toeplitz structure which can be capitalized in the estimation process. By
doing so, problem (5) becomes

$\begin{array}[]{ll}\textrm{MLE:}\quad\underset{{\mathbf{R}}\in
Toep,{\mathbf{R}}\succ 0}{\rm minimize}\>f({\mathbf{R}})\end{array}$ (6)

where $Toep$ is used to denote the set of Toeplitz matrices of size $m\times
m$. The above problem has two constraints: a structural constraint and a
positive definite constraint. Even though the structural constraint is convex,
the non-convexity of the objective function makes Problem (6) challenging to
solve and no analytical solution seems to be available. In the following two
iterative solution procedures for (6) are designed exploiting the MM
principle. Briefly, the MM technique mainly consists of two steps - 1)
Constructing a surrogate function $g({\mathbf{R}}|{\mathbf{R}}_{t})$ (where
${\mathbf{R}}_{t}$ is the estimate of ${\mathbf{R}}$ at the $t^{th}$
iteration) for the objective function in (6) and 2) Minimizing the resulting
surrogate problem at each iteration. For more details, [24, 25, 26, 27]
provide an in-depth discussion on MM based algorithms.

##  III Algorithms for Toeplitz covariance matrix estimation

In this section, ATOM1 and ATOM2 are proposed to tackle the MLE problem of TSC
matrix. Both leverage the MM principle (applied to an equivalent reformulation
of the MLE problem) and differ in the way they construct and handle the
surrogate minimization problem. ATOM1 solves the surrogate optimization using
ADMM while ATOM2 tackles it using alternating projection or Dykstra’s
algorithm. Hence the computational complexity and proof of convergence of the
procedures are established. Finally, the extension of ATOM2 to deal with
additional covariance constraints along with the Toeplitz structure is
provided.
To begin with, let us present the following two lemmas which pave the way for
the development of ATOM1 and ATOM2.

###### Lemma III.1

Consider the following two optimization problems

$\begin{array}[]{ll}{\mathbf{R}}^{*}=\underset{{\mathbf{R}}\in
Toep,{\mathbf{R}}\succ 0}{\rm
arg\>min}\>\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{\mathbf{R}}^{-1}{\mathbf{y}}_{i}+\log|{\mathbf{R}}|\end{array}$
($P_{1}$)
$\begin{array}[]{ll}{\boldsymbol{X}}^{*}=\underset{{{\boldsymbol{X}}}\in
Toep,{\boldsymbol{X}}\succ 0}{\rm arg\>min}\>\log|{\boldsymbol{X}}|\>\>\>{\rm
s.t.}\quad\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{{\boldsymbol{X}}}^{-1}{\mathbf{y}}_{i}\leq
1\end{array}$ ($P_{2}$)

Problems $(P_{1})$ and $(P_{2})$ are equivalent and ${\mathbf{R}}^{*}$ can be
obtained via ${\boldsymbol{X}}^{*}$ using the following relationship

${\mathbf{R}}^{*}=\dfrac{{\boldsymbol{X}}^{*}}{m}$ (7)

###### Proof:

The proof can be found in Appendix A. ∎

Before proceeding with the next important lemma, it is worth pointing out that
Lemma III.1 holds true even if the Toeplitz structural constraint in $(P_{1})$
and $(P_{2})$ is replaced by an arbitrary closed-cone set.

###### Lemma III.2

Given a concave differentiable111For a non-differentiable function, the
inequality in (8) can be relaxed to $f({\boldsymbol{X}})\leq
f({\boldsymbol{X}}_{t})+{\rm{Tr}}\left({\boldsymbol{G}}({\boldsymbol{X}}_{t})^{H}({\boldsymbol{X}}-{\boldsymbol{X}}_{t})\right)$,
where ${\boldsymbol{G}}({\boldsymbol{X}}_{t})$ is the subgradient of the
concave function $f({\boldsymbol{X}})$ at ${\boldsymbol{X}}_{t}$.  function
$f({\boldsymbol{X}}):\mathbb{H}^{m\times m}\rightarrow\mathbb{R}$, it can be
majorized as

$\begin{array}[]{ll}f({\boldsymbol{X}})\leq
f({\boldsymbol{X}}_{t})+{\rm{Tr}}\left(\nabla
f({\boldsymbol{X}}_{t})^{H}({\boldsymbol{X}}-{\boldsymbol{X}}_{t})\right)\end{array}$
(8)

where ${\boldsymbol{X}}={{\boldsymbol{X}}}_{t}$ is the value of
${\boldsymbol{X}}$ at the $t^{th}$ iteration. The upper bound to
$f({\boldsymbol{X}})$ is linear and differentiable with respect to (w.r.t.)
${\boldsymbol{X}}$.

###### Proof:

Since $f({\boldsymbol{X}})$ is a concave function w.r.t. ${\boldsymbol{X}}$,
(8) stems from linearizing $f({\boldsymbol{X}})$ via its first order Taylor
expansion. ∎

Leveraging Lemma III.1, a minimizer of Problem ($P_{1}$), i.e., the MLE
problem in (6), can be obtained by solving Problem ($P_{2}$) and using the
relationship in (7). Therefore, in the following, let us focus on solving
Problem ($P_{2}$). Now, since $\log|{\boldsymbol{X}}|$ in ($P_{2}$) is a
concave function w.r.t. ${\boldsymbol{X}}$ [28], it can be majorized using
Lemma III.2 to get the following surrogate function

$\begin{array}[]{ll}g({\boldsymbol{X}}|{\boldsymbol{X}}_{t})=\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{\boldsymbol{X}}\right)\>{\rm
s.t.}\quad\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{{\boldsymbol{X}}}^{-1}{\mathbf{y}}_{i}\leq
1\end{array}$ (9)

Therefore, at any iteration, given ${\boldsymbol{X}}_{t}$, the MM method
demands for the solution of the following surrogate minimization problem

$\begin{array}[]{ll}{\boldsymbol{X}}_{t+1}=\underset{{\boldsymbol{X}}\in
Toep,{\boldsymbol{X}}\succ 0}{\rm
arg\>min}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{\boldsymbol{X}}\right)\>{\rm
s.t.}\quad\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{{\boldsymbol{X}}}^{-1}{\mathbf{y}}_{i}\leq
1\end{array}$ (10)

This a convex optimization problem and can be cast as a SemiDefinite Program
(SDP) (the interested reader may refer to Appendix B for details)

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}\in Toep}}{\rm
minimize}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)\\\
{\rm subject\
to}\>\>\boldsymbol{\hat{X}}={\mathbf{I}}\otimes{\boldsymbol{X}}\succeq\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\end{array}$
(11)

where $\boldsymbol{\hat{X}}$ is a block diagonal matrix with the diagonal
blocks equal to ${\boldsymbol{X}}$ and $\bar{{\mathbf{y}}}$ is a vector
obtained by stacking the vectors
${\mathbf{y}}_{1},{\mathbf{y}}_{2},\cdots,{\mathbf{y}}_{n}$. However, the
computational complexity necessary to handle SDP using interior point methods
is $\mathcal{O}(m^{3}n^{3})$ [29]. In the following, an iterative algorithm
based on ADMM is proposed to solve the surrogate minimization problem in (11).

###  III-A ATOM1

The surrogate minimization problem in (11) is solved using ADMM [30, 31]. To
do so, an auxiliary variable ${\boldsymbol{U}}$ is introduced in (11) and the
problem is re-written in the equivalent form

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}}\in
Toep,\>{\boldsymbol{U}}}{\rm
minimize}&\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)\\\
{\rm subject\
to}&\quad{{\mathbf{I}}\otimes{\boldsymbol{X}}}-{\boldsymbol{U}}=\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\\\
&\quad\quad{\boldsymbol{U}}\succeq\boldsymbol{0}\end{array}$ (12)

The augmented Lagrangian for (12) is

$\begin{array}[]{ll}L_{\rho}({\boldsymbol{X}},{\boldsymbol{U}},\hat{\boldsymbol{\lambda}})=\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)+{\rm{Tr}}\left({(\hat{{\boldsymbol{\lambda}}})}^{T}\left(({\mathbf{I}}\otimes{\boldsymbol{X}})-{\boldsymbol{U}}-\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\right)\right)+\dfrac{\rho}{2}\|({\mathbf{I}}\otimes{\boldsymbol{X}})-{\boldsymbol{U}}-\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\|_{F}^{2}\end{array}$
(13)

where $\rho>0$ is the penalty parameter and
$\hat{{\boldsymbol{\lambda}}}\succeq 0$ is the Lagrange multiplier of size
$mn\times mn$. The iterative steps of ADMM algorithm are

$\begin{array}[]{ll}{\boldsymbol{U}}^{t}_{k+1}=\underset{{{\boldsymbol{U}}\succeq
0}}{\rm
arg\>min}\>{\rm{Tr}}\left({(\hat{{\boldsymbol{\lambda}}}^{t}_{k})}^{T}\left(({\mathbf{I}}\otimes{\boldsymbol{X}}_{k}^{t})-{\boldsymbol{U}}\right)\right)+\dfrac{\rho}{2}\|({\mathbf{I}}\otimes{\boldsymbol{X}}_{k}^{t})-{\boldsymbol{U}}-\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\|_{F}^{2}\\\
\end{array}$ (14)
$\begin{array}[]{ll}{\boldsymbol{X}}^{t}_{k+1}=\underset{{{\boldsymbol{X}}\in
Toep}}{\rm
arg\>min}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)+{\rm{Tr}}\left({(\hat{{\boldsymbol{\lambda}}}^{t}_{k})}^{T}\left(({\mathbf{I}}\otimes{\boldsymbol{X}})\right)\right)+\dfrac{\rho}{2}\|({\mathbf{I}}\otimes{\boldsymbol{X}})-{\boldsymbol{U}}^{t}_{k+1}-\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\|_{F}^{2}\\\
\end{array}$ (15)
$\begin{array}[]{ll}\hat{{\boldsymbol{\lambda}}}^{t}_{k+1}=\hat{{\boldsymbol{\lambda}}}^{t}_{k}+\rho\left(({\mathbf{I}}\otimes{\boldsymbol{X}}^{t}_{k+1})-{\boldsymbol{U}}^{t}_{k+1}-\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\right)\end{array}$
(16)

where $(\cdot)^{t}_{k}$ is used to denote the $k$-th iteration of the ADMM
algorithm in correspondence of the $t$-th MM outer loop. Problems (14) and
(15) have closed-form solutions which can be computed via the projection of
appropriate matrices onto the respective feasible sets. Indeed, Problem (14)
can be equivalently cast as

$\begin{array}[]{ll}{\boldsymbol{U}}^{t}_{k+1}=\underset{{{\boldsymbol{U}}}}{\rm
arg\>min}\>\|{\boldsymbol{U}}-{\boldsymbol{\Psi}}\|_{F}^{2}\\\
\quad\quad\quad\quad{\rm subject\ to}\quad{{\boldsymbol{U}}}\succeq 0\\\
\end{array}$ (17)

${\boldsymbol{\Psi}}=\left(\left({\mathbf{I}}\otimes{\boldsymbol{X}}^{t}_{k}\right)+\dfrac{(\hat{{\boldsymbol{\lambda}}}^{t}_{k})}{\rho}-\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\right)$.
Hence, solving (14) is tantamount to performing the orthogonal projection of
the matrix ${\boldsymbol{\Psi}}$ onto the set of the Positive SemiDefinite
(PSD) matrices which can be computed as
${\boldsymbol{U}}^{t}_{k+1}=\tilde{{\boldsymbol{V}}}\rm{max}({\rm{diag}}(\tilde{{\boldsymbol{U}}}),{\boldsymbol{0}})\tilde{{\boldsymbol{V}}}^{H}$,
where ${\rm{diag}}(\tilde{{\boldsymbol{U}}})$ and $\tilde{{\boldsymbol{V}}}$
are the matrices containing the eigenvalues and the corresponding eigenvectors
of the matrix ${\boldsymbol{\Psi}}$, respectively. Similarly, the update step
of ${\boldsymbol{X}}$ in (15) can be re-written as:

$\begin{array}[]{ll}{\boldsymbol{X}}^{t}_{k+1}=\underset{{{\boldsymbol{X}}\in
Toep}}{\rm arg\>min}\>\|{\boldsymbol{X}}-{\boldsymbol{\Lambda}}\|_{F}^{2}\\\
\end{array}$ (18)

where
${\boldsymbol{\Lambda}}=\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}\left(({{{\boldsymbol{W}}_{i}}})^{t}_{k+1}-\dfrac{1}{\rho}{({{\boldsymbol{\lambda}}_{i}})^{t}_{k}}-\dfrac{1}{\rho
n}({{\boldsymbol{X}}}_{t})^{-1}\right)$, $({\boldsymbol{W}}_{i})^{t}_{k+1}$
and $({{\boldsymbol{\lambda}}_{i}})^{t}_{k}$ are the $i^{th}$ diagonal block
(of size $m\times m$) of
$\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}+{\boldsymbol{U}}^{t}_{k+1}$
and $\hat{{\boldsymbol{\lambda}}}^{t}_{k}$, respectively. The solution to (18)
requires the projection of the matrix ${\boldsymbol{\Lambda}}$ onto the set of
Toeplitz matrices which can be calculated by replacing the elements along the
diagonals of the matrix ${\boldsymbol{\Lambda}}$ with the average along the
respective diagonal elements. Before concluding, it is worth pointing out that
since the surrogate minimization problem in (11) is convex, it is guaranteed
that ADMM converges to a supposed existing optimal unique solution to (11)
(see Section $3.2$ in [32], [33]). The pseudocode of the proposed algorithm is
shown below.

Algorithm 1 : Pseudocode of ATOM 1 algorithm
---
Input: Samples
$\\{{\mathbf{y}}_{1},{\mathbf{y}}_{2},\cdots,{\mathbf{y}}_{n}\\}$ and $\rho$
Initialize: Set ${t},{k}=0$. Initialize ${{{\boldsymbol{X}}}_{0}}$,
${{{\boldsymbol{X}}}^{t}_{0}}$ and $\hat{{\boldsymbol{\lambda}}}^{t}_{0}$.
Repeat: Given ${\boldsymbol{X}}_{t}$, perform the $t+1$-th step.
$k\leftarrow 0$
Repeat: Given ${{{\boldsymbol{X}}}}^{t}_{k}$ and
$\hat{{\boldsymbol{\lambda}}}^{t}_{k}$, perform the $k+1$-th step.
1) Obtain ${\boldsymbol{U}}^{t}_{k+1}$ by projecting the matrix
${\boldsymbol{\Psi}}=\left(\left({\mathbf{I}}\otimes{\boldsymbol{X}}^{t}_{k}\right)+\dfrac{{\hat{{\boldsymbol{\lambda}}}}^{t}_{k}}{\rho}-\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\right)$
onto the
set of PSD matrices.
2) Obtain ${\boldsymbol{X}}^{t}_{k+1}$ by projecting the matrix
${\boldsymbol{\Lambda}}=\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}\left(({\boldsymbol{W}}_{i})^{t}_{k+1}-\dfrac{1}{\rho}{({\boldsymbol{\lambda}}_{i})^{t}_{k}}-\dfrac{1}{\rho
n}({{\boldsymbol{X}}}_{t})^{-1}\right)$
onto the set of Toeplitz matrices.
3)
$\hat{{\boldsymbol{\lambda}}}^{t}_{k+1}=\hat{{\boldsymbol{\lambda}}}^{t}_{k}+\rho\left(({\mathbf{I}}\otimes{\boldsymbol{X}}^{t}_{k+1})-{\boldsymbol{U}}^{t}_{k+1}-\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\right)$
until convergence.
Set ${{\boldsymbol{X}}}_{t+1}={{\boldsymbol{X}}}^{t}_{k+1}$
$t\leftarrow t+1$
until convergence
Output: ${\mathbf{R}}_{{\rm{\bf{ATOM1}}}}=\dfrac{{\boldsymbol{X}}^{*}}{m}$,
where ${\boldsymbol{X}}^{*}$ denotes the value of ${\boldsymbol{X}}_{t}$ at
the completion of the outer loop.

From Algorithm 1 it can be seen that ATOM1 requires initialization of the
matrices ${\boldsymbol{X}}_{0}$, ${\boldsymbol{X}}^{t}_{0}$ and
$\hat{{\boldsymbol{\lambda}}}^{t}_{0}$. ${\boldsymbol{X}}_{0}$ can be set
using the initialization scheme discussed in [18] and, as $t=0$,
${\boldsymbol{X}}^{t}_{0}$ can be set equal to ${\boldsymbol{X}}_{0}$ while
$\hat{{\boldsymbol{\lambda}}}^{t}_{0}$ can be constructed as
$\hat{{\boldsymbol{\lambda}}}^{t}_{0}={\boldsymbol{V}}{\boldsymbol{V}}^{T}$,
where the elements of ${\boldsymbol{V}}$ are drawn randomly from a uniform
distribution over $[0,1]^{mn}$. For $t\geq 1$, the matrices
${\boldsymbol{X}}^{t}_{0}$ and $\hat{{\boldsymbol{\lambda}}}^{t}_{0}$ can be
initialized with their last value after convergence at the previous ADMM
iteration, respectively. Another input parameter required by ATOM1 is the
penalty parameter $\rho$, introduced during the formation of the Augmented
Lagrangian of the ADMM algorithm. It is shown in [32], that the ADMM algorithm
converges for any value of $\rho>0$. However, the numerical stability and the
convergence rate of the algorithm depends on the choice of $\rho$. Simulation
results have highlighted that for $\rho=m$, the ADMM algorithm is stable for
different values of $n$ and $m$. Hence, unless otherwise stated, in all the
numerical analysis the penalty parameter $\rho=m$ is used.

###  III-B Computational complexity of ATOM1

ATOM1 is iterative in nature with two loops - the outer loop updates the
Toeplitz matrix ${\boldsymbol{X}}_{t}$ while the inner loop solves the
surrogate minimization problem using ADMM. Note that in the inner loop, one
requires to compute the matrix
$\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}$ \- which is iteration
independent and hence can be pre-computed. Moreover, this matrix can be
evaluated efficiently leveraging an equivalent formulation of Problems (10)
and (11). In this respect, first obtain the Cholesky decomposition of the
matrix ${\mathbf{R}}_{\textrm{SCM}}$, i.e., compute
${\mathbf{R}}_{\textrm{SCM}}=\bar{{\mathbf{R}}}\bar{{\mathbf{R}}}^{H}$ (where
$\bar{{\mathbf{R}}}$ is a lower triangular matrix) with a complexity of
$\mathcal{O}(m^{3})$. Then, Problem (11) becomes

$\begin{array}[]{ll}{\boldsymbol{X}}_{t+1}=\underset{{{\boldsymbol{X}}}\in
Toep,{\boldsymbol{X}}\succeq 0}{\rm
arg\>min}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)\>\>{\rm
s.t.}\quad\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}[\bar{{\mathbf{R}}}]_{i}^{H}{{\boldsymbol{X}}}^{-1}[\bar{{\mathbf{R}}}]_{i}\leq
1\end{array}$ (19)

Casting the above problem as an SDP it follows that

$\begin{array}[]{ll}{\boldsymbol{X}}_{t+1}=\underset{{{\boldsymbol{X}}}\in
Toep,\boldsymbol{\hat{X}}}{\rm
arg\>min}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)\>\>{\rm
s.t.}\quad\boldsymbol{\hat{X}}\geq\bar{{\mathbf{R}}}_{\textrm{SCM}}\end{array}$
(20)

where
$\bar{{\mathbf{R}}}_{\textrm{SCM}}=\bar{{\mathbf{r}}}\bar{{\mathbf{r}}}^{H}$
and the vector $\bar{{\mathbf{r}}}\in\mathbb{C}^{m^{2}\times 1}$ is obtained
by vectorizing the matrix $\bar{{\mathbf{R}}}$. Comparing the constraint of
Problem (20) and (11), it can be observed that instead of computing the matrix
$\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}$, one needs to evaluate
$\bar{{\mathbf{R}}}_{\textrm{SCM}}$. Exploiting the sparseness of the vector
$\bar{{\mathbf{r}}}$, the matrix $\bar{{\mathbf{R}}}_{\textrm{SCM}}$ can be
determined with a cost of
$\mathcal{O}\left(\left(\dfrac{m(m+1)}{2}\right)^{2}\right)$ [34]. Let us now
discuss the complexity related to the outer and inner loops of ATOM1. In this
regard, note that the inner loop of ATOM1 requires the computation of the
matrix $({{\boldsymbol{X}}_{t}})^{-1}$ \- which is outer loop iteration
dependent. Therefore, this matrix can be computed once in the outer loop.
Consequently, apart from the computations of the inner loop, an outer loop
cycle requires only the estimation of the matrix
$({{\boldsymbol{X}}_{t}})^{-1}$. This in-turn can be efficiently computed
using Gohberg–Semencul formula with a complexity $\mathcal{O}(m\>{\rm{log}}m)$
[2]. The computational complexity of an inner loop cycle is related to the
projection of ${\boldsymbol{\Psi}}$ onto the set of PSD matrices and
projection of ${\boldsymbol{\Lambda}}$ onto the set of Toeplitz matrices. The
cost of projection of ${\boldsymbol{\Lambda}}$ onto the set of Toeplitz
matrices is mainly dictated by the computation of average of the elements
along the respective diagonals of ${\boldsymbol{\Lambda}}$. Hence, the cost of
matrix projection onto the Toepltiz matrices set is $\mathcal{O}(m^{2})$.
Next, the projection of ${\boldsymbol{\Psi}}$ onto the set of PSD matrices
mainly involves the computation of the eigenvalues and eigenvectors of the
matrix ${\boldsymbol{\Psi}}$ \- whose corresponding complexity is
$\mathcal{O}({m^{6}})$ [34]. Therefore, the computational complexity of ATOM1
is $\mathcal{O}(\eta{m^{6}})$ where $\eta$ is the total number of inner loop
iterations required by the algorithm to converge. A drawback of ATOM1 is that
it has a restriction in handling additional constraints on the covariance
matrix. This is because ATOM1 implements ADMM algorithm at every inner
iteration which requires the optimization problem to exhibit the standard form
[32, 35]

$\begin{array}[]{ll}\underset{{{\boldsymbol{Z}}},{\boldsymbol{E}}}{\rm
minimize}&f({\boldsymbol{Z}})+g({\boldsymbol{E}})\\\ {\rm subject\
to}&\>\>{\mathbf{A}}_{1}{\boldsymbol{Z}}+{\mathbf{A}}_{2}{\boldsymbol{E}}={\boldsymbol{C}}\end{array}$
(21)

where $f({\boldsymbol{Z}})$, $g({\boldsymbol{E}})$ are convex functions and
${\mathbf{A}}_{1}$, ${\mathbf{A}}_{2}$, ${\boldsymbol{C}}$ are matrices of
appropriate dimensions, respectively. Therefore, to incorporate additional
inequality constraints (such as an upper bound on the condition number of the
matrix ${\boldsymbol{Z}}$), one needs to replace the inequality constraint
with an equality constraint. This can be done by introducing a slack variable
to the existing optimization variables ${\boldsymbol{Z}}$ and
${\boldsymbol{E}}$. However, there is no convergence guarantee of ADMM when
there are more than two optimization variables [36]. This issue is addressed
in the following by proposing another low complexity algorithm, referred to as
ATOM2, to solve Problem (11).

###  III-C ATOM 2

To develop a computationally efficient and flexible estimation procedure
capable of including additional constraints, a different surrogate function is
constructed for the objective function in Problem ($P_{2}$). To this end, let
us once again refer to Problem (11) cast in the equivalent (but
computationally convenient) form

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}\in
Toep,{\boldsymbol{\hat{X}}}}}{\rm
minimize}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)\\\
{\rm subject\
to}\quad\boldsymbol{\hat{X}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\end{array}$
(22)

and suitably constructing its surrogate function. In this respect, adding and
subtracting $\textrm{Tr}({{\boldsymbol{X}}}^{2})$ in the objective function
yields

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}}\in
Toep,\boldsymbol{\hat{X}}}{\rm
minimize}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)+\textrm{Tr}({{\boldsymbol{X}}}^{2})-\textrm{Tr}({{\boldsymbol{X}}}^{2})\\\
{\rm subject\
to}\quad\boldsymbol{\hat{X}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\end{array}$
(23)

Note that $-\textrm{Tr}({{\boldsymbol{X}}}^{2})$ is a concave function of
${{\boldsymbol{X}}}$ and hence given a feasible solution
${\boldsymbol{X}}_{t}$, it is possible to construct the surrogate function
$\tilde{g}({\boldsymbol{X}}|{\boldsymbol{X}}_{t})$ for the above problem using
Lemma III.2

$\begin{array}[]{ll}\tilde{g}({\boldsymbol{X}}|{\boldsymbol{X}}_{t})=\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)+\textrm{Tr}({{\boldsymbol{X}}}^{2})-2\textrm{Tr}({{\boldsymbol{X}}}{{{\boldsymbol{X}}}}_{t})\end{array}$
(24)

Note that $\tilde{g}({\boldsymbol{X}}|{\boldsymbol{X}}_{t})$ is a surrogate to
a surrogate function. Nonetheless, since
$\tilde{g}({\boldsymbol{X}}|{\boldsymbol{X}}_{t})$ is a tighter surrogate for
$g({\boldsymbol{X}}|{\boldsymbol{X}}_{t})$, it can be viewed as a direct
surrogate for the objective function in Problem ($P_{2}$). Then, at any
iteration $t$, given ${\boldsymbol{X}}_{t}$, an optimal solution to the
following surrogate minimization problem

$\begin{array}[]{ll}{\boldsymbol{X}}_{t+1}=\underset{{{\boldsymbol{X}}}\in
Toep}{\rm arg\>min}\>\|{\boldsymbol{X}}-{\boldsymbol{B}}\|_{F}^{2}\\\ \hskip
42.67912pt{\rm subject\
to}\quad\boldsymbol{\hat{X}}={\mathbf{I}}\otimes{\boldsymbol{X}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\end{array}$
(25)

is required, where
${\boldsymbol{B}}={{{\boldsymbol{X}}}}_{t}-0.5({{\boldsymbol{X}}}_{t})^{-1}$.
An efficient way to solve Problem (25) is now devised. To this end, let us
cast (25) in terms of a single block diagonal matrix ${\mathbf{T}}$ as follows

$\begin{array}[]{ll}\underset{{{\mathbf{T}}}}{\rm
minimize}\>\|{\mathbf{T}}-{\mathbf{I}}\otimes{\boldsymbol{B}}\|_{F}^{2}\\\
{\rm subject\ to}\>\>{\mathbf{T}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\\\
\hskip 49.79231pt{{\mathbf{T}}\>\textrm{is a block diagonal matrix with the
same diagonal Toeplitz blocks}}\end{array}$ (26)

Problem (26) seeks for the point ${\mathbf{T}}$ belonging to the intersection
of the two sets - the set of block diagonal matrices with the same diagonal
Toeplitz blocks and the set defined by the Linear Matrix Inequality (LMI)
${\mathbf{T}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}$ \- closest to
${\mathbf{I}}\otimes{\boldsymbol{B}}$. Since both the sets are convex, one can
use the technique of alternating projection or Projection Onto the Convex Sets
(POCS) to solve Problem (26) [37, 38, 39]. To describe this algorithm, let us
denote by ${\mathcal{P}_{Toep}}({\boldsymbol{\Psi}})$ and
${\mathcal{P}_{LMI}}({\boldsymbol{\Psi}})$ the orthogonal projection of a
matrix ${\boldsymbol{\Psi}}$ onto the set of block diagonal matrices with the
same diagonal Toeplitz blocks and onto the set defined by
${\mathbf{T}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}$, respectively. Both
${\mathcal{P}_{Toep}}({\boldsymbol{\Psi}})$ and
${\mathcal{P}_{LMI}}({\boldsymbol{\Psi}})$ can be obtained as follows:

  1. 1.

Calculation of ${\mathcal{P}_{Toep}}({\boldsymbol{\Psi}})$:
Given a matrix ${\boldsymbol{\Psi}}$, the orthogonal projection of the matrix
onto the set of block diagonal matrices with the same diagonal Toeplitz blocks
can be obtained in two steps. First the averaging of the diagonal blocks of
the matrix ${\boldsymbol{\Psi}}$ is performed. Then, the elements along each
diagonal of the matrix obtained in the previous step are averaged.

  2. 2.

Calculation of ${\mathcal{P}_{LMI}}({\boldsymbol{\Psi}})$:
Given a matrix ${\boldsymbol{\Psi}}$, its orthogonal projection onto the set
defined by LMI ${\mathbf{T}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}$ can be
obtained as follows. First, calculate the EigenValue Decomposition (EVD) of
the matrix ${\boldsymbol{\Psi}}-\bar{{\mathbf{R}}}_{\textrm{SCM}}$, i.e.,
obtain
$[\bar{{\boldsymbol{U}}},\bar{{\boldsymbol{V}}}]={\rm{eig}}({\boldsymbol{\Psi}}-\bar{{\mathbf{R}}}_{\textrm{SCM}})$,
where $\bar{{\boldsymbol{U}}}$ and $\bar{{\boldsymbol{V}}}$ are matrices
containing the eigenvalues and eigenvectors of the spectral decomposition,
respectively. Then, the orthogonal projection
${\mathcal{P}_{LMI}}({\boldsymbol{\Psi}})$ is given by
$\bar{{\boldsymbol{V}}}{\rm{max}}(\bar{{\boldsymbol{U}}},{\boldsymbol{0}})\bar{{\boldsymbol{V}}}^{H}+\bar{{\mathbf{R}}}_{\textrm{SCM}}$.

According to POCS method, given an initial value
${\mathbf{T}}_{0}={\mathbf{I}}\otimes{\boldsymbol{B}}$, first compute
${\boldsymbol{Y}}_{k+1}={\mathcal{P}_{Toep}}({\mathbf{T}}_{k})$ and then,
using ${\boldsymbol{Y}}_{k+1}$, determine
${\mathbf{T}}_{k+1}={\mathcal{P}_{LMI}}({\boldsymbol{Y}}_{k+1})$, i.e., the
starting point of the next iteration. Hence, POCS finds a sequence of iterates
${\mathbf{T}}_{k}$ by alternatingly projecting between the two convex sets.
However, as reported in [40], POCS may suffer from slow convergence. A
refinement of POCS is Dykstra’s projection [41, 42] which finds a point
closest to ${\mathbf{I}}\otimes{\boldsymbol{B}}$ by adding correction vectors
${\mathbf{P}}_{k}$ and ${\boldsymbol{Q}}_{k}$ before every projection, which
in-turn ensures convergence of sequence ${\mathbf{T}}_{k+1}$ to the optimal
solution ${{\mathbf{T}}}^{*}$ [41]. The pseudocode of Dykstra’s algorithm is
shown in Algorithm 2.

Algorithm 2: Pseudocode of Dykstra’s algorithm
---
Input:
${\boldsymbol{B}}={{{\boldsymbol{X}}}}_{t}-0.5({{\boldsymbol{X}}}_{t})^{-1}$
Initialize: Initialize ${\mathbf{T}}_{0}={\mathbf{I}}\otimes{\boldsymbol{B}}$,
${{{\mathbf{P}}}_{0}}={\boldsymbol{0}}$ and
${{{\boldsymbol{Q}}}_{0}}={\boldsymbol{0}}$
Repeat:
1)
${\boldsymbol{Y}}_{k}={\mathcal{P}_{Toep}}({\mathbf{T}}_{k}+{\mathbf{P}}_{k})$
2) ${\mathbf{P}}_{k+1}={\mathbf{T}}_{k}+{\mathbf{P}}_{k}-{\boldsymbol{Y}}_{k}$
3)
${\mathbf{T}}_{k+1}={\mathcal{P}_{LMI}}({\boldsymbol{Y}}_{k}+{\boldsymbol{Q}}_{k})$
4)
${\boldsymbol{Q}}_{k+1}={\boldsymbol{Y}}_{k}+{\boldsymbol{Q}}_{k}-{\mathbf{T}}_{k+1}$
$k\leftarrow k+1$
until convergence
Output: ${{\mathbf{T}}}^{*}={\mathbf{T}}_{k+1}$.

Once the optimal solution ${\mathbf{T}}^{*}$ is obtained via Dykstra’s
projection, the matrix ${{\boldsymbol{X}}}_{t+1}$ can be constructed from its
upper left block of size $m\times m$. This process is repeated until the whole
algorithm, i.e., including the outer loop, converges. Finally, using the
relation in (7), the optimal Toeplitz matrix ${\mathbf{R}}^{*}$ of Problem
($P_{1}$) is attained. ATOM2 is summarized in Algorithm 3.

Algorithm 3: Pseudocode of ATOM2
---
Input: Samples
$\\{{\mathbf{y}}_{1},{\mathbf{y}}_{2},\cdots,{\mathbf{y}}_{n}\\}$
Initialize: Set _t_ = 0. Initialize ${\boldsymbol{X}}_{0}$.
Repeat: Given ${{{\boldsymbol{X}}}_{t}}$ perform the $t+1$-th step.
Run Dykstra’s algorithm with
${\boldsymbol{B}}={{{\boldsymbol{X}}}}_{t}-0.5({{\boldsymbol{X}}}_{t})^{-1}$
as input and obtain ${\boldsymbol{X}}_{t+1}$ from
the upper left block of ${\mathbf{T}}^{*}$.
$t\leftarrow t+1$
until convergence
Output: ${\mathbf{R}}_{{\rm{\bf{ATOM2}}}}=\dfrac{{\boldsymbol{X}}^{*}}{m}$,
where ${\boldsymbol{X}}^{*}$ denotes the value of ${\boldsymbol{X}}_{t}$ at
the completion of the outer loop.

ATOM2 requires the initialization of the matrix ${\boldsymbol{X}}$. In this
respect, a similar scheme as in ATOM1 is followed, i.e., at each outer
iteration, the initial guess required to determine ${\boldsymbol{X}}_{t+1}$ in
the inner loop is set equal to ${\boldsymbol{X}}_{t}$. Furthermore, at $t=0$,
the initial value of the matrix ${\boldsymbol{X}}_{0}$ is obtained according
to the initialization scheme discussed in [18].

###  III-D Computational complexity of ATOM2

Similar to ATOM1, ATOM2 is also an iterative algorithm with outer and inner
loops. The outer loop updates the Toeplitz matrix ${\boldsymbol{X}}_{t}$ and
the inner loop implements the Dykstra’s algorithm - which requires the
computation of the matrices $\bar{{\mathbf{R}}}_{\textrm{SCM}}$ and
$({{\boldsymbol{X}}_{t}})^{-1}$. Note that, the former matrix is iteration
independent and therefore can be pre-computed with a cost of
$\mathcal{O}\left(\left(\dfrac{m(m+1)}{2}\right)^{2}\right)$. On the other
hand, the matrix $({{\boldsymbol{X}}_{t}})^{-1}$ is outer loop iteration
dependent and therefore can be computed once in the outer loop. Consequently,
apart from the inner loop computations, the outer loop demands only the
computation of $({{\boldsymbol{X}}_{t}})^{-1}$ \- which can be computed
efficiently with complexity $\mathcal{O}(m\>{\rm{log}}m)$. Meanwhile, the
computational load of the inner loop stems from the computation of EVD of the
matrix $({\boldsymbol{Y}}_{k}+{\boldsymbol{Q}}_{k})$ plus a rank one matrix
$\bar{{\mathbf{R}}}_{\textrm{SCM}}$ \- which has a complexity of about
$\mathcal{O}(m^{6})$. However, as already pointed, at the expense of slightly
slower convergence, one can also implement the inner loop of ATOM2 using POCS
algorithm. In this case, it is demanded to compute the EVD of a structured
matrix, i.e., a Toeplitz structured matrix ${\boldsymbol{Y}}_{k}$ plus a rank
one matrix $\bar{{\mathbf{R}}}_{\textrm{SCM}}$. This can be computed
efficiently in $\mathcal{O}(m^{2}\>{\rm{log}}m+m^{2})$ [43, 44]. Therefore,
for large values of covariance matrix dimension $m$ one can implement the
inner loop of ATOM2 algorithm using POCS algorithm. Hence, the total
computational cost of ATOM2 and ATOM2 using POCS (hereafter refereed to as
ATOM2-POCS) is $\mathcal{O}(\eta(m^{6}))+\mathcal{O}(m^{4})$ and
$\mathcal{O}(\eta(m^{2}\>{\rm{log}}(m)))+\mathcal{O}(m^{4})$, respectively and
$\eta$ is used to represent the total number of inner loop iterations required
by the algorithm to converge. In Table I, the computational complexity of
ATOM1, ATOM2, and ATOM2-POCS is compared with that of the state-of-the-art
iterative algorithms [18, 15]. Unlike the proposed algorithms, the state-of-
the art methods are single loop iteration algorithms. Therefore, in the case
of [18, 15] $\eta$ is used to represent the number of iterations required by
the algorithm to converge. Inspection of Table I shows that ATOM1 and ATOM2
(whose inner loop is implemented using Dykstra’s algorithm) has the highest
complexity when compared to ATOM2-POCS, MELT and EM. As to ATOM2-POCS, the per
cycle computational cost is only $m$ times larger than MELT and EM.
Nevertheless, it is worth anticipating that this complexity increase is
complemented by a superior performance in terms of covariance matrix MSE and
achieved SINR.

ATOM1 |  ATOM2 |  ATOM2-POCS |  MELT [18] |  EM [15]
---|---|---|---|---
$\mathcal{O}(\eta{m^{6}})$ |  $\mathcal{O}(\eta{m^{6}})$ |  $\mathcal{O}(\eta(m^{2}\>{\rm{log}}(m)))+\mathcal{O}(m^{4})$ |  $\mathcal{O}\left(\eta(m{\rm{log}}(m))\right)$ |  $\mathcal{O}\left(\eta(m{\rm{log}}(m)\right))$
Table I: Comparison among computational complexity of ATOM1 and ATOM2 with
other state-of-the-art iterative algorithms.

###  III-E Proof of convergence

In this subsection, the proof of convergence of ATOM1 and ATOM2 is
established. In this regard, it is worth pointing out that both the algorithms
differ in the way they construct and optimize the surrogate function for the
Problem ($P_{2}$). Nonetheless, since ATOM1 and ATOM2 are based on the MM
framework, the proof of convergence based on the following lemma will hold for
both algorithms.
Before stating the lemma, let us first introduce the first-order optimality
condition for minimizing a function over a convex constraint set. A point
${\boldsymbol{X}}$ is a stationary point of $f(\cdot)$ if
$f^{\prime}({\boldsymbol{X}};{\boldsymbol{D}})\geq 0$ for all
${\boldsymbol{D}}$ such that
${\boldsymbol{X}}+{\boldsymbol{D}}\in\mathcal{C}$, where $\mathcal{C}$ is the
convex constraint set and $f^{\prime}({\boldsymbol{X}};{\boldsymbol{D}})$ is
the directional derivative of $f(\cdot)$ at point ${\boldsymbol{X}}$ in
direction ${\boldsymbol{D}}$ and is defined as [26]

$\begin{array}[]{ll}f^{\prime}({\boldsymbol{X}};{\boldsymbol{D}})=\underset{\lambda\downarrow
0}{\lim}\>\textrm{inf}\>\dfrac{f({\boldsymbol{X}}+\lambda{\boldsymbol{D}})-f({\boldsymbol{X}})}{\lambda}\end{array}$
(27)

Based on the following lemma, both ATOM1 and ATOM2 are guaranteed to converge
to a stationary point of Problem ($P_{2}$)

###### Lemma III.3

Denoting by $\\{{\boldsymbol{X}}_{t}\\}$ the sequence generated by either
ATOM1 or ATOM2, then the objective function of Problem ($P_{2}$) monotonically
decreases along the iterations. Besides, any positive definite cluster point
to $\\{{\boldsymbol{X}}_{t}\\}$ is a stationary point to Problem ($P_{2}$).

###### Proof:

See Appendix C for details. ∎

Before concluding this subsection, it is worth pointing out that $m\geq n/2$
is a sufficient condition to guarantee that $\\{{\boldsymbol{X}}_{t}\\}$ is a
bounded sequence whose cluster points are all positive definite matrices.
Indeed, under such condition, the sequence $\\{{\boldsymbol{X}}_{t}\\}$
belongs to a compact subset of positive definite matrices.

###  III-F Extensions of ATOM2

The extension of ATOM2 to handle additional constraints other than the
Toeplitz structure in the covariance estimation process is now addressed. In
particular, it is shown that ATOM2 can be generalized to account for the
following scenarios: Banded Toeplitz matrices, Toeplitz-block-Toeplitz
matrices, low rank matrices with Toeplitz structure plus a scalar matrix, and
finally, matrices with Toeplitz structure satisfying a condition number
constraint. However, as discussed at the end of Subsection III-B, ATOM1 cannot
be directly extended to tackle the general constraints as for instance an
upper bound requirement to the condition number.
1\. MLE of low rank Toeplitz covariance term plus a scalar matrix
In radar applications and in some array signal processing contexts [20, 18,
45], the covariance matrix ${\mathbf{R}}$ to be estimated often exhibits the
following structure

${\mathbf{R}}={\mathbf{R}}_{c}+\sigma{\mathbf{I}}$ (28)

where ${\mathbf{R}}_{c}$ denotes the covariance matrix of the interference
(clutter or jammers) and $\sigma{\mathbf{I}}$ is the covariance matrix of the
thermal noise (where $\sigma>0$ represents the spectral density level of the
white noise and is unknown). Also, ${\mathbf{R}}_{c}$ usually exhibits a low
rank Toeplitz structure with known rank $r\leq m$ [46]. As a result, the
covariance estimation problem can be cast as

$\begin{array}[]{ll}\underset{{\mathbf{R}}\in Toep,{\mathbf{R}}_{c}\succeq
0,\sigma\geq 0}{\rm
minimize}\>\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{\mathbf{R}}^{-1}{\mathbf{y}}_{i}+\log|{\mathbf{R}}|\\\
\hskip 5.69046pt{\rm subject\
to}\quad\>\>\>{{\mathbf{R}}}={\mathbf{R}}_{c}+\sigma{\mathbf{I}}\\\ \hskip
59.75095pt\quad{\rm Rank}({{\mathbf{R}}_{c}})\leq r\end{array}$ (29)

Following a similar line of reasoning as that used in the development of
ATOM2, the constrained ML estimate can be obtained by leveraging Lemma III.1
and solving the following equivalent problem

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}\in
Toep,{\boldsymbol{X}}_{c}\succeq 0,\boldsymbol{\hat{X}}},\sigma\geq 0}{\rm
minimize}\>\log|{\boldsymbol{X}}|\\\ \hskip 14.22636pt{\rm subject\
to}\quad\>\>\boldsymbol{\hat{X}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\\\
\quad\quad\quad\quad\quad\quad\quad{{\boldsymbol{X}}}={\boldsymbol{X}}_{c}+\sigma{\mathbf{I}}\\\
\quad\quad\quad\quad\quad\quad\quad{\rm Rank}({{\boldsymbol{X}}_{c}})\leq
r\end{array}$ (30)

In this respect, the MM framework is exploited and, using Lemma III.2, the
following surrogate optimization problem is obtained

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}\in
Toep},{\boldsymbol{X}}_{c}\succeq 0,\boldsymbol{\hat{X}},\sigma\geq 0}{\rm
minimize}\>\|{\boldsymbol{X}}-{\boldsymbol{B}}\|_{F}^{2}\\\ \hskip
14.22636pt{\rm subject\
to}\quad\>\>\boldsymbol{\hat{X}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\\\
\quad\quad\quad\quad\quad\quad\quad{{\boldsymbol{X}}}={\boldsymbol{X}}_{c}+\sigma{\mathbf{I}}\\\
\quad\quad\quad\quad\quad\quad\quad{\rm Rank}({{\boldsymbol{X}}_{c}})\leq
r\end{array}$ (31)

where
${\boldsymbol{B}}={{\boldsymbol{X}}}_{t}-0.5({{\boldsymbol{X}}}_{t})^{-1}$. It
can be tackled using alternating projection among three subsets at each
iteration:
$\mathcal{S}_{1}=\left\\{{\boldsymbol{X}}:\boldsymbol{\hat{X}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\right\\}$,
$\mathcal{S}_{2}=\left\\{{\boldsymbol{X}}:{{\boldsymbol{X}}}={\boldsymbol{X}}_{c}+\sigma{\mathbf{I}},{\boldsymbol{X}}_{c}\succeq
0,{\rm Rank}({{\boldsymbol{X}}_{c}})\leq r,\sigma\geq 0\right\\}$, and
$\mathcal{S}_{3}=Toep$. Although there are no convergence guarantees for
alternating projections in the presence of non-convex sets (the set of low
rank matrices is non-convex set), this strategy usually proves effective, for
instance see [47] for the frame design problem. The projection of a given
matrix onto $\mathcal{S}_{1}$ and $\mathcal{S}_{3}$ cane be performed as
discussed in Section III B. Let $\hat{{\boldsymbol{B}}}$ represent the matrix
obtained after projecting the output of the previous algorithm iteration onto
$\mathcal{S}_{1}$ (the algorithm is initialized with the input matrix
${\boldsymbol{B}}$). Then, the projection of $\hat{{\boldsymbol{B}}}$ onto
$\mathcal{S}_{2}$ requires the solution of the following problem

$\begin{array}[]{ll}\underset{{\boldsymbol{X}},{\boldsymbol{X}}_{c}\succeq
0,\sigma\geq 0}{\rm
minimize}\>\>\>\|{\boldsymbol{X}}-\hat{{\boldsymbol{B}}}\|_{F}^{2}\\\ \hskip
0.28436pt{\rm subject\
to}\>\>{{\boldsymbol{X}}}={\boldsymbol{X}}_{c}+\sigma{\mathbf{I}}\\\
\quad\quad\quad\quad\quad{\rm Rank}({{\boldsymbol{X}}_{c}})\leq r\end{array}$
(32)

Let ${\boldsymbol{\beta}}=[\beta_{1},\beta_{2},\cdots,\beta_{m}]^{T}$ and
${\boldsymbol{\gamma}}=[\gamma_{1},\gamma_{2},\cdots,\gamma_{m}]^{T}$ be the
eigenvalues of the matrix $\hat{{\boldsymbol{B}}}$ and ${\boldsymbol{X}}$,
sorted in decreasing order, respectively. Now, denoting by ${\boldsymbol{U}}$
a unitary matrix such that
$\hat{{\boldsymbol{B}}}={\boldsymbol{U}}\rm{diag}(\boldsymbol{\beta}){\boldsymbol{U}}^{H}$
and exploiting the unitary invariance of the objective and the constraint
functions of (32), it follows that an optimal solution to (32) is
$\bar{{\boldsymbol{X}}}={\boldsymbol{U}}\rm{diag}({\boldsymbol{\gamma}}^{*}){\boldsymbol{U}}^{H}$
with
$\boldsymbol{\gamma}^{*}=[\gamma_{1}^{*},\cdots,\gamma_{r}^{*},\sigma,...\sigma]^{T}\in\mathbb{R}^{m}$,
the optimal solution to

$\begin{array}[]{ll}\underset{\gamma_{1}\geq\gamma_{r}\geq\sigma\geq
0}{\rm{minimize}}\displaystyle\sum_{i=1}^{r}(\gamma_{i}-\beta_{i})^{2}+\displaystyle\sum_{i=r+1}^{m}(\sigma-\beta_{i})^{2}\end{array}$
(33)

The optimal solution to (33) is

$\displaystyle\gamma^{*}_{i}$ $\displaystyle=$
$\displaystyle\beta_{i},i=1,2\cdots,r,$ (34) $\displaystyle\sigma^{*}$
$\displaystyle=$
$\displaystyle\dfrac{1}{m-r}\displaystyle\sum_{i=r+1}^{m}\beta_{i}.$ (35)

Hence
$\bar{{\boldsymbol{X}}}=\beta_{1}{\boldsymbol{u}}_{1}{\boldsymbol{u}}_{1}^{H}+\cdots\beta_{r}{\boldsymbol{u}}_{r}{\boldsymbol{u}}_{r}^{H}+\sigma^{*}({\boldsymbol{u}}_{r+1}{\boldsymbol{u}}_{r+1}^{H}+\cdots+{\boldsymbol{u}}_{m}{\boldsymbol{u}}_{m}^{H})$.
2\. MLE of banded Toeplitz covariance matrix
The covariance matrix is constrained to exhibit a banded Toeplitz structure of
bandwidth $b$ (see [17, 48] for relevant applications). For instance, assuming
a bandwidth $b=2$ and dimension $m=5$ the covariance matrix enjoys the
following structure

${\mathbf{R}}=\begin{bmatrix}r_{1}&r_{2}&r_{3}&0&0\\\
r^{*}_{2}&r_{1}&r_{2}&r_{3}&0\\\ r^{*}_{3}&r^{*}_{2}&r_{1}&r_{2}&r_{3}\\\
0&r^{*}_{3}&r^{*}_{2}&r_{1}&r_{2}\\\
0&0&r^{*}_{3}&r^{*}_{2}&r_{1}\end{bmatrix}$

Then, the MLE problem for banded Toeplitz covariance matrix can be formulated
as

$\begin{array}[]{ll}\underset{{\mathbf{R}}\in Band-Toep,\>{\mathbf{R}}\succ
0}{\rm
minimize}\>\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{\mathbf{R}}^{-1}{\mathbf{y}}_{i}+\log|{\mathbf{R}}|\end{array}$
(36)

where $Band-Toep$ is used to denote the set of banded Toeplitz matrices.
Exploiting again Lemma III.1, the above problem can be cast in the following
equivalent form

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}\in Band-
Toep,\boldsymbol{\hat{X}}}}{\rm minimize}\>\log|{\boldsymbol{X}}|\\\ \hskip
17.07164pt{\rm subject\
to}\quad\>\boldsymbol{\hat{X}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\end{array}$
(37)

Hence (37) is handled via MM framework solving the following surrogate
minimization problem

$\begin{array}[]{ll}\underset{{{\mathbf{T}}}}{\rm
minimize}\>\|{\mathbf{T}}-{\mathbf{I}}\otimes{\boldsymbol{B}}\|_{F}^{2}\\\
{\rm subject\ to}\>\>{\mathbf{T}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\\\
\hskip 49.79231pt{{\mathbf{T}}\>\textrm{is a block diagonal matrix with the
same diagonal banded Toeplitz blocks}}\end{array}$ (38)

The above problem involves two convex sets: the set defined by the LMI
${\mathbf{T}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}$ and the set of block
diagonal matrices with each block being equal and having a banded Toeplitz
structure with bandwidth $b$. Consequently, Dykstra’s projection algorithm or
POCS can be used to solve Problem (38). The projection of a matrix onto the
LMI set can be calculated as discussed earlier in Subsection III-C. The
projection of a matrix $\hat{{\boldsymbol{\Psi}}}$ onto the set of block
banded Toeplitz matrices can be obtained in two steps. In the first step,
averaging the diagonal blocks of the matrix $\hat{{\boldsymbol{\Psi}}}$ is
performed. In the next step, the elements outside the bandwidth are replaced
with zeros. In contrast, the elements within the bandwidth are replaced with
the average of the respective diagonal elements of the matrix obtained in the
first step [49].
3\. MLE of Toeplitz-block-Toeplitz covariance matrix
In space-time adaptive processing radar applications, the covariance matrix
exhibits a Toeplitz-block-Toeplitz (TBT) structure. This refers to a block
Toeplitz matrix with each block having a Toeplitz structure [50, 51]. An
example of a TBT-structured covariance matrix with $p$ blocks is shown below

${\mathbf{R}}=\begin{bmatrix}{\mathbf{R}}_{1}&{\mathbf{R}}_{2}&\dots&{\mathbf{R}}_{p}\\\
{\mathbf{R}}^{*}_{2}&{\mathbf{R}}_{1}&\dots&{\mathbf{R}}_{2}\\\
\vdots&\ddots&\ddots&\vdots\\\
{\mathbf{R}}^{*}_{p}&\dots&{\mathbf{R}}^{*}_{2}&{\mathbf{R}}_{1}\end{bmatrix}$

The MLE problem of TBT covariance matrix is formulated as

$\begin{array}[]{ll}\underset{{\mathbf{R}}\in TBT,{\mathbf{R}}\succ 0}{\rm
minimize}\>\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{\mathbf{R}}^{-1}{\mathbf{y}}_{i}+\log|{\mathbf{R}}|\end{array}$
(39)

where the notation $TBT$ is used to indicate the set of Toeplitz-block-
Toeplitz matrices. The minimizer of (39) is obtained by solving at any given
step the following surrogate optimization problem

$\begin{array}[]{ll}\underset{{{\mathbf{T}}}}{\rm
minimize}\>\|{\mathbf{T}}-{\mathbf{I}}\otimes{\boldsymbol{B}}\|_{F}^{2}\\\
{\rm subject\ to}\>\>{\mathbf{T}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\\\
\hskip 49.79231pt{{\mathbf{T}}\>\textrm{is a block diagonal matrix with the
same diagonal TBT blocks}}\end{array}$ (40)

Problem (40) exhibits two constraints - $1$) a LMI constraint and $2$) a
structural constraint - where the optimization variable ${\mathbf{T}}$ is
confined to be a block diagonal matrix with each block being equal and having
a TBT structure. Since both the constraints are convex, Dykstra’s projection
or POCS can be applied to solve Problem (40). The projection of a matrix onto
the LMI set can be calculated as discussed earlier in Section III B. The
projection of a given matrix $\bar{{\boldsymbol{\Psi}}}$ onto the set of
matrices having the TBT constraint can be obtained as follows. First, a matrix
$\bar{{\boldsymbol{\Psi}}}_{\rm{Avg}}$ is obtained by averaging the diagonal
blocks of the matrix $\bar{{\boldsymbol{\Psi}}}$. Next $p$ matrices are
obtained by averaging the diagonal blocks of the matrix
$\bar{{\boldsymbol{\Psi}}}_{\rm{Avg}}$ from the preceding step. Finally, each
of the $p$ matrices are projected onto the Toeplitz set as described in
Section III C [49].
4\. MLE of Toeplitz covariance matrix with condition number constraint
In this extension, a condition number constraint $\kappa$ (assumed to be
known) in addition to the Toeplitz structural constraint on the covariance
matrix (see for applications [52, 53]) is considered. The MLE problem can be
formulated as

$\begin{array}[]{ll}\underset{{\mathbf{R}}\in Toep,{\mathbf{R}}\succ 0}{\rm
minimize}\>\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{\mathbf{R}}^{-1}{\mathbf{y}}_{i}+\log|{\mathbf{R}}|\\\
{\rm subject\
to}\>\>\>\dfrac{\lambda_{max}({\mathbf{R}})}{\lambda_{min}({\mathbf{R}})}\leq\kappa\end{array}$
(41)

To get a minimizer of the above problem, the steps used in ATOM2 are borrowed
to arrive at the following surrogate optimization problem

$\begin{array}[]{ll}\underset{{{\mathbf{T}}}}{\rm
minimize}\>\|{\mathbf{T}}-{\mathbf{I}}\otimes{\boldsymbol{B}}\|_{F}^{2}\\\
{\rm subject\ to}\>\>{\mathbf{T}}\succeq\bar{{\mathbf{R}}}_{\textrm{SCM}}\\\
\hskip 49.79231pt{{\mathbf{T}}\>\textrm{is a block diagonal matrix with the
same diagonal Toeplitz blocks}}\\\ \hskip 49.79231pt\textrm{whose condition
number is upper bounded by\>}\kappa\end{array}$ (42)

The above problem includes three constraints: the LMI constraint, the
structural constraint which confines the covariance matrix to follow a
Toeplitz structure, and the condition number constraint. All the three
constraints are convex and therefore, Dykstra’s projection or POCS can be
applied to solve Problem (42). The projection of a given matrix
$\tilde{{\boldsymbol{\Psi}}}$ onto
$\dfrac{\lambda_{max}(\tilde{{\boldsymbol{\Psi}}})}{\lambda_{min}(\tilde{{\boldsymbol{\Psi}}})}\leq\kappa$
is now discussed. Given a matrix $\tilde{{\boldsymbol{\Psi}}}$, let
${\boldsymbol{V}}{\boldsymbol{\Gamma}}{\boldsymbol{V}}^{H}$ denote its EVD
with
${\boldsymbol{\Gamma}}=\rm{diag}([\gamma_{1},\gamma_{2},\cdots,\gamma_{m}]^{T})$
where $\gamma_{1}\geq\gamma_{2}\geq\cdots\geq\gamma_{m}$ and the columns of
${\boldsymbol{V}}$ contain the corresponding eigenvectors. Then, the
projection of $\tilde{{\boldsymbol{\Psi}}}$ onto the condition number
constraint is given by [21]:
${\boldsymbol{V}}{\boldsymbol{\Lambda}}{\boldsymbol{V}}^{H}$ where
${\boldsymbol{\Lambda}}={\rm{diag}}({\boldsymbol{\lambda}}(u^{*}))$ and
${\boldsymbol{\lambda}}(u)=[\lambda_{1}(u),\lambda_{2}(u),\cdots,\lambda_{n}(u)]^{T}$
with

$\lambda_{i}(u)=\min\left(\kappa
u,\max\left(\gamma_{i},\max\left(0,u\right)\right)\right)\>\>i=1,2,\cdots,m$
(43)

and $u^{*}$ can be obtained by solving a convex optimization problem with a
linear computational complexity, see [52] for more details. The projection
onto the convex sets defined by the LMI and Toeplitz structure can be
evidently handled as discussed in Subsection III-C.

##  IV Cramer-Rao Lower bound calculation

In this section, the Cramer-Rao Lower Bound (CRLB) is derived for three
structured covariance matrices, namely, the Toeplitz, the banded Toeplitz, and
the TBT structures. The CRLB provides a lower bound on the variance of any
unbiased estimator [54]. To proceed further, let $\theta$ represent the real
value vector parametrizing a given covariance matrix structure of interest.
The specific definition of $\theta$ is provided in the next subsections for
each case study. Then, the CRLB is the inverse of the Fisher Information
matrix (FIM) whose $(i,k)^{th}$ element is

$\begin{array}[]{ll}[{\mathbf{F}}]_{i,k}=\textbf{E}\left[\frac{\partial^{2}\log
f({\mathbf{y}};{\mathbf{R}})}{\partial\theta_{i}\partial\theta_{k}}\right]\end{array}$
(44)

where $f({\mathbf{y}};{\mathbf{R}})$ is the likelihood function i.e.,
$f({\mathbf{y}};{\mathbf{R}})=\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{\mathbf{R}}^{-1}{\mathbf{y}}_{i}+\log|{\mathbf{R}}|$.
The $(i,k)^{th}$ element of the FIM can be computed using the _Slepian–Bangs
formula_ [2]

$[{\mathbf{F}}]_{i,k}=n\textrm{Tr}\left({\mathbf{R}}^{-1}\frac{\partial{\mathbf{R}}}{\partial\theta_{i}}{\mathbf{R}}^{-1}\frac{\partial{\mathbf{R}}}{\partial\theta_{k}}\right)$
(45)

In the following subsections, the lower bounds on
${\mathbf{R}}({\boldsymbol{\theta}})$ for some covariance structures
considered in the manuscript are derived.

###  IV-A Toeplitz matrix

As the entries of the TSC matrix are completely characterized by its first
row, i.e., $[r_{1},r_{2},\cdots r_{m}]^{T}$, the covariance matrix
${\mathbf{R}}\in\mathbb{H}^{m\times m}$ can be parameterized by
${\boldsymbol{\theta}}=[r_{1},\Re(r_{2}),\cdots\Re(r_{m}),\\\
\Im(r_{2}),...,\Im(r_{m})]^{T}\in\mathbb{R}^{2m-1}$ where $\Re(r_{i})$ and
$\Im(r_{i})$ denotes the real and imaginary parts of $r_{i}$, respectively.
Then, the covariance matrix ${\mathbf{R}}$ can be expressed in terms of
${\boldsymbol{\theta}}$ and basis matrices ${\boldsymbol{B}}^{\rm{Toep}}_{g}$
(defined as in (47)), $g=1,2,\cdots,m$ [16]

$\begin{array}[]{ll}{\mathbf{R}}=\displaystyle\sum_{g=1}^{m}\theta_{g}{\Re}({\boldsymbol{B}}^{\rm{Toep}}_{g})+j\displaystyle\sum_{g=m+1}^{2m-1}\theta_{g}{\Im}({\boldsymbol{B}}^{\rm{Toep}}_{g-m+1})\end{array}$
(46)

The $(i,k)^{th}$ element of the matrix ${\boldsymbol{B}}^{\rm{Toep}}_{g}$ is
given as

$\begin{array}[]{ll}[{\boldsymbol{B}}^{\rm{Toep}}_{g}]_{i,k}=\begin{cases}1&i-k=g-1=0\\\
1+j&k-i=g-1\neq 0\\\ 1-j&i-k=g-1\neq 0\\\
0&\rm{otherwise}\end{cases}\end{array}$ (47)

Using (46), $\frac{\partial{\mathbf{R}}}{\partial\theta_{i}}$ can be obtained
as

$\frac{\partial{\mathbf{R}}}{\partial\theta_{i}}=\begin{cases}\Re({\boldsymbol{B}}^{\rm{Toep}}_{i})&1\leq
i\leq m\\\ j\Im({\boldsymbol{B}}^{\rm{Toep}}_{i-m+1})&m+1\leq i\leq
2m-1\end{cases}$

Substituting $\frac{\partial{\mathbf{R}}}{\partial\theta_{i}}$ in (45), yields
the FIM for Toeplitz covariance matrix. Then, the CRLB on $r_{i}$ is

$\begin{cases}[{\mathbf{F}}^{-1}]_{i,i}&i=1\\\
[{\mathbf{F}}^{-1}]_{i,i}+[{\mathbf{F}}^{-1}]_{i+m-1,i+m-1}&i\neq
1\end{cases}$

###  IV-B Banded Toepltiz matrix

In the case of banded Toeplitz matrix with bandwidth $b$, the first row of the
covariance matrix ${\mathbf{R}}\in\mathbb{H}^{m\times m}$ has only $b+1$ non-
zero terms. Therefore, ${\mathbf{R}}$ can be parameterized via
${\boldsymbol{\theta}}=[r_{1},\Re(r_{2}),\cdots\Re(r_{b+1})\\\
,\Im(r_{2}),...,\Im(r_{b+1})]^{T}\in\mathbb{R}^{2b+1}$. Besides ${\mathbf{R}}$
can be expressed in terms of basis matrices ${\boldsymbol{B}}^{\rm{Toep}}_{g}$
and real coefficients ${\boldsymbol{\theta}}$

$\begin{array}[]{ll}{\mathbf{R}}=\displaystyle\sum_{g=1}^{b+1}\theta_{g}{\Re}({\boldsymbol{B}}^{\rm{Toep}}_{g})+j\displaystyle\sum_{g=b+2}^{2b+1}\theta_{g}{\Im}({\boldsymbol{B}}^{Toep}_{g-b})\end{array}$
(48)

and consequently

$\frac{\partial{\mathbf{R}}}{\partial\theta_{i}}=\begin{cases}{\Re}({\boldsymbol{B}}^{\rm{Toep}}_{i})&1\leq
i\leq b+1\\\ j{\Im}({\boldsymbol{B}}^{\rm{Toep}}_{i-b})&b+2\leq i\leq
2b+1\end{cases}$

Substituting $\frac{\partial{\mathbf{R}}}{\partial\theta_{i}}$ in (45), yields
the FIM for banded Toeplitz covariance matrix and the CRLB on $r_{i}$ can be
obtained as

$\begin{cases}[{\mathbf{F}}^{-1}]_{i,i}&i=1\\\
[{\mathbf{F}}^{-1}]_{i,i}+[{\mathbf{F}}^{-1}]_{i+b,i+b}&i\neq 1\end{cases}$

###  IV-C Toeplitz-block-Toeplitz matrix

A TBT matrix with $p$ blocks of size $l$ can be specified by its first row
i.e., $[r_{1},{r_{2}},\cdots,{r_{l}},\cdots,r_{(p-1)(2l-1)+1},\\\
{r_{(p-1)(2l-1)+2}},\cdots,{r_{(p-1)(2l-1)+l}}]^{T}$. As a consequence, a TBT
matrix can be expressed as follows

$\begin{array}[]{ll}{\mathbf{R}}^{\rm{TBT}}={\boldsymbol{C}}_{0}\otimes{\mathbf{R}}_{0}+\displaystyle\sum_{w=1}^{{p-1}}\left(\left({\boldsymbol{C}}_{w}\otimes{\mathbf{R}}_{w}\right)+\left({\boldsymbol{C}}_{w}^{T}\otimes{\mathbf{R}}^{*}_{w}\right)\right)\end{array}$
(49)

where $(i,k)^{th}$ element of the matrix
${\boldsymbol{C}}_{w}\in\mathbb{R}^{l\times l}$ is given by

$[{\boldsymbol{C}}_{w}]_{i,k}=\begin{cases}1&i-k=w\\\
0&\rm{otherwise}\end{cases}$

and ${\mathbf{R}}_{w}$ is

$\begin{array}[]{ll}{\mathbf{R}}_{w}=\displaystyle\sum_{g=1}^{l}\theta_{(g+qw)}{\Re}({\boldsymbol{B}}^{\rm{Toep}}_{g})+j\displaystyle\sum_{g={l+1}}^{2l-1}\theta_{(g+qw)}{\Im}({\boldsymbol{B}}^{Toep}_{g-l+1})\end{array}$
(50)

where
${\boldsymbol{\theta}}=[r_{1},\Re({r_{2}}),\cdots,\Re({r_{l}}),\Im({r_{2}}),\cdots,\Im({r_{l}}),\cdots,r_{(p-1)(2l-1)+1},\Re({r_{(p-1)(2l-1)+2}}),\cdots,\Re({r_{(p-1)(2l-1)+l}}),\\\
\Im({r_{(p-1)(2l-1)+2}}),\cdots,\Im({r_{(p-1)(2l-1)+l}})]^{T}\in\mathbb{R}^{(2l-1)p}$
and ${\boldsymbol{B}}^{\rm{Toep}}_{g}$ of size $\l\times l$ is defined in
(47). Let $\theta_{i}$ belong to the $z^{th}$ block and $q=2l-1$, then
$\frac{\partial{\mathbf{R}}}{\partial\theta_{i}}$ is given by

$\frac{\partial{\mathbf{R}}}{\partial\theta_{i}}=\begin{cases}{\boldsymbol{C}}_{0}\otimes\Re({\boldsymbol{B}}^{\rm{Toep}}_{i})&1\leq
i\leq{l},z=0\\\
j{\boldsymbol{C}}_{0}\otimes\Im({\boldsymbol{B}}^{\rm{Toep}}_{i-l+1})&{l}+1\leq
i\leq q,z=0\\\
{\boldsymbol{C}}_{z}\otimes\Re({\boldsymbol{B}}^{\rm{Toep}}_{i-qz})+{\boldsymbol{C}}_{z}^{T}\otimes{\Re}({\boldsymbol{B}}^{\rm{Toep}}_{i-qz})&1+qz\leq
i\leq{l}+qz,z\neq 0\\\
j{\boldsymbol{C}}_{z}\otimes{\Im}({\boldsymbol{B}}^{\rm{Toep}}_{i-qz-l+1})-j{\boldsymbol{C}}_{z}^{T}\otimes{\Im}({{\boldsymbol{B}}^{*}}^{{\rm{Toep}}}_{i-qz-l+1})&{l}+1+qz\leq
i\leq q+qz,z\neq 0\\\ \end{cases}$

Substituting $\frac{\partial{\mathbf{R}}}{\partial\theta_{i}}$ in (45), the
FIM for TBT covariance matrix is obtained. Then, the CRLB on $r_{i}$ is

$\begin{cases}[{\mathbf{F}}^{-1}]_{2i-(z+1),2i-(z+1)}&i=1,l+1,\cdots,(p-1)l+1\\\
[{\mathbf{F}}^{-1}]_{i+2z,i+2z}+[{\mathbf{F}}^{-1}]_{i+l-1+2z,i+l-1+2z}&i\neq
1,l+1,\cdots,(p-1)l+1\\\ \end{cases}$

##  V Numerical Simulations

In this section, the performance of the proposed Toeplitz covariance matrix
estimators ATOM1 and ATOM2 is analyzed in comparison with some state-of-the-
art algorithms via numerical simulations using the MSE metric. In particular,
ATOM1 and ATOM2 are compared with the EM-based algorithm [15], MELT [18], and
the SCM estimator. Furthermore, ATOM2 is also compared with techniques for
estimating the covariance matrix with banded Toeplitz, TBT, and condition
number constraints. Finally, the performance of the estimators is evaluated in
terms of maximum achievable SINR in a typical radar signal processing
scenario. All the algorithms are implemented in MATLAB using a PC with $2.40$
GHz processor and $16$ GB RAM.

###  V-A Assessment of iterative algorithms convergence for on-grid and off-
grid frequencies

In this simulation, the convergence of the ATOM1 and ATOM2 (whose inner loop
was implemented via Dykstra’s algorithm) is assessed in comparison with MELT
and EM algorithms. To this end, the data ${\mathbf{y}}_{k}$ are generated
according to the following model

$\begin{array}[]{ll}{\mathbf{y}}_{k}=\sqrt{{\mathbf{R}}}{\boldsymbol{n}}_{k},k=1,2,\cdots,n\end{array}$
(51)

where ${\boldsymbol{n}}_{k}$’s are drawn randomly from a zero-mean circularly
symmetric Gaussian distribution with independent and identically distributed
entries with unit mean square value and the size of the training data is
$n=460$. The convergence of the iterative algorithms for the two different
experimental setups are now evaluated. In the first setup, the true underlying
Toeplitz covariance matrix ${\mathbf{R}}$ of dimension $m=6$ is constructed by
choosing the $2$-nd, $3$-rd, $5$-th, $7$-th, $8$-th and the $11$-th columns of
the DFT matrix with $L=2m-1$ in (4), whose corresponding frequencies and
amplitudes are: $[0.5712,1.1424,2.2848,3.4272,3.9984,5.7120]$ rad/sec and
$[3,6,4,1,7,5]$, respectively. Fig. 1.a shows the Objective Value (OV) of
Problem $(P_{1})$, i.e., the negative log likelihood, versus the number of
iterations for the first experimental setup. It can be seen that all the
algorithms numerically decrease the underlying cost function and converge to
the same OV. For reference, the OV of Problem $(P_{2})$, i.e.,
$\log|{\boldsymbol{X}}|$ versus the number of iterations of ATOM1 and ATOM2 is
also shown in Fig. 1.b. This figure indicates that the proposed algorithms
monotonically decrease the OV - which is expected since these algorithms
optimize Problem $(P_{2})$ using the MM framework. In the second experimental
setup, the true underlying Toeplitz covariance matrix is constructed such that
one of the frequencies does not lie on the Fourier grid. The experimental
setting is same as in case study 1, with the exception that the Fourier
frequency $2.2848$ rad/sec is replaced with $2.5$ rad/sec. Fig. 2.a shows the
OV of Problem $(P_{1})$ versus the number of iterations for the second case
study. Fig. 2.b indicates the corresponding OV of Problem $(P_{2})$ versus the
number of iterations of the proposed methods. From Fig. 2.a, it can be seen
that while MELT and EM converges to a value of $6.16$, ATOM1 and ATOM2
converges to $5.81$. Therefore, when one of the frequencies does not lie on
the Fourier grid, the state-of-the-art iterative algorithms converge to a
larger value of the negative log-likelihood as compared to the proposed
algorithms. This is because unlike the counterparts, the proposed algorithms
estimate the Toeplitz covariance matrix without reparametrizing it via the CE
technique and thus covers the whole set of Toeplitz covariance matrices.

(a)

(b)

Figure 1: On-grid frequencies scenario for $m=6$ and $n=460$. a) Negative log-
likelihood $f({\mathbf{R}})$ vs. outer iterations; b) $\log|{\boldsymbol{X}}|$
vs. outer iterations

(a)

(b)

Figure 2: Off-grid frequencies scenario for $m=6$ and $n=460$. a) Negative
log-likelihood $f({\mathbf{R}})$ vs. outer iterations; b)
$\log|{\boldsymbol{X}}|$ vs. outer iterations

Next, the computational time of the proposed techniques is compared with the
other state-of-the-art iterative algorithms. The data samples
${\mathbf{y}}_{k}$ are generated using (51) and $n=50$ samples. The Toeplitz
covariance matrix ${\mathbf{R}}$ is generated according to the model in (4)
with $L=2m-1$ for four different values of $m$ \- $4$, $8$, $16$ and $32$. The
iterative algorithms have been run until the following condition is met

$\dfrac{\|{\mathbf{R}}_{t}-{\mathbf{R}}_{t+1}\|_{F}}{\|{\mathbf{R}}_{t}\|_{F}}\leq
10^{-4}$ (52)

or until the maximum number of iterations (set equal to $1000$) is reached.
Table. II compares the average computational time (averaged over $50$ Monte-
Carlo trials) of the different algorithms. The results show that ATOM2 has a
smaller execution time than ATOM1. This is because the inner loop of ATOM1
(implemented via ADMM) requires an higher number of iterations and hence a
longer run time to converge than ATOM2 inner loop involved in Dykstra’s
algorithm. From the table, it can also be seen that MELT has the least
computational time. Nevertheless, although the proposed algorithms have a
longer computational time than the counterparts, the obtained covariance
matrix, performs, in general, better in terms of MSE than the solution
provided by MELT and EM.

Table II: Comparison of the average run time (in seconds) of the iterative algorithms Covariance matrix dimension $m$ |  ATOM1 |  ATOM2 |  MELT[18] |  EM[15]
---|---|---|---|---
4 |  0.41 |  0.01 |  0.002 |  0.008
8 |  2.62 |  0.054 |  0.006 |  0.025
16 |  33.16 |  0.506 |  0.031 |  0.045
32 |  154.56 |  1.348 |  0.1043 |  0.2208

###  V-B MSE vs $n$ for Toeplitz covariance matrix

Data ${\mathbf{y}}_{k}\in\mathbb{C}^{15}$ are generated in the same way as in
the previous subsection, modeling ${\mathbf{R}}$ as in (1) with specific
frequencies and amplitudes. The number of samples $n$ ranges between $50$ and
$500$ in steps of $50$. As in the previous subsection, two different
experiments are conducted assuming that the true Toeplitz covariance matrix is
generated using on-grid and off-grid frequencies, respectively. The algorithms
are compared using MSE metric

$\textrm{MSE}({\mathbf{r}})=\textbf{E}\left[\dfrac{1}{m}\displaystyle\sum_{i=1}^{m}\left|r_{i}-\hat{r}_{i}\right|^{2}\right]$
(53)

where ${\mathbf{r}}$ and $\hat{{\mathbf{r}}}$ indicate the first row of the
true covariance matrix ${\mathbf{R}}$ and its estimate $\hat{{\mathbf{R}}}$,
respectively. The number of Monte-Carlo trials used in this experiment to
evaluate numerically (53) is equal to $100$. The sum of CRLB on each $r_{i}$
derived in Section IV is used as benchmark. Fig. 3.a shows that in the first
experiment222The frequencies used in the first experiment are:
$[0.2167,0.6500,1.0833,1.3,1.5166,1.9500,2.3833,2.8166,3.2499,3.6832,\\\
4.1166,4.5499,4.9832,5.4165,5.8499]$ rad/sec. Their corresponding amplitudes
increase linearly from $1$ to $15$ with a unit step., ATOM 1, ATOM 2, EM and
MELT achieve a similar performance and reach the least MSE as compared to SCM.
Fig. 3.b highlight that in the second experiment (where the frequency $0.2167$
rad/sec is replaced with off-grid frequency $0.5$ rad/sec), ATOM 1 and ATOM 2
perform better than MELT, EM and SCM and, unlike the counterparts, achieve the
CRLB when $n$ is sufficiently large. Furthermore, MELT and EM exhibit similar
MSE’s. Hence, unlike MELT and EM, the performance of ATOM 1 and ATOM 2 do not
depend on the Fourier gridding.

(a)

(b)

Figure 3: MSE vs. number of samples $n$ for Toeplitz covariance matrix. a) on-
grid frequencies; b) off-grid frequencies

###  V-C MSE vs $n$ for banded Toeplitz covariance matrix

In this subsection, the performance of ATOM2 is analyzed with the state-of-
the-art algorithms when the covariance matrix is constrained to be a banded
Toeplitz matrix. The ground truth banded Toeplitz matrix of dimension $m=15$
and bandwidth $b=6$ is constructed by alternately projecting a random
Hermitian matrix onto the set of banded Toeplitz matrices and the set of PSD
matrices. Hence, the data ${\mathbf{y}}_{k}$ are generated using (51). The
number of samples $n$ ranges between $50$ and $500$ in steps of $50$. The
number of Monte-Carlo experiments is the same as in the previous simulation
and different algorithms are compared with in terms of MSE using the CRLB as
benchmark. For each Monte-Carlo experiment the same ground truth covariance
matrix was employed. Fig. 4 highlights that the proposed algorithm provide a
lower MSE than the state-of-the-art algorithms. This is because unlike the
counterparts, the new method is developed without the CE. Also, as in the
previous subsection, from Fig. 4 it can be observed that MELT and EM achieve a
similar performance.

Figure 4: MSE(${\mathbf{r}}$) vs. number of samples $n$ for banded Toeplitz
covariance matrix

###  V-D MSE vs $n$ for TBT covariance matrix

ATOM2 is compared with the counterparts when the covariance matrix is
constrained to be a TBT matrix. The ground truth covariance matrix of
dimension $m=16$ with block size equal to $4$ is constructed by alternately
projecting a random Hermitian matrix onto the set of TBT matrices and the set
of PSD matrices. Using (51), the data ${\mathbf{y}}_{k}$ are generated using
the ground truth TBT covariance matrix. The number of samples $n$ ranges
between $50$ and $500$ in steps of $50$. The number of Monte-Carlo experiments
is the same as in the previous simulation and the algorithms are compared in
terms of the MSE metric and with CRLB as benchmark. Fig. 5 shows that ATOM2
uniformly achieves the least MSE. As already highlighted the superior
performance of the proposed method stems from the design criterion which does
not require reparametrizing the covariance matrix using the CE.

Figure 5: MSE(${\mathbf{r}}$) vs. number of samples $n$ for TBT covariance
matrix

###  V-E MSE vs $n$ for Toeplitz covariance matrix with condition number
constraint

ATOM2 is compared with the state-of-the-art algorithms when the Toeplitz
covariance matrix is forced to comply with a condition number constraint. The
data ${\mathbf{y}}_{k}$ and ground truth Toeplitz covariance matrix of
dimension $m=15$ are generated using the Fourier frequencies mentioned in
Subsection IV-B. The condition number $\kappa$ is set equal to the condition
number of the true underlying Toeplitz covariance matrix [21]. Fig. 6.a shows
the MSE versus the number of samples $n$ for the on-grid frequency case study.
Inspection of the figure highlights that the proposed algorithm, capitalizing
on the a-priori knowldge on the covariance condition number, outperforms the
state-of-the-art counterparts even when the frequencies lie on the Fourier
grid. In Fig. 6.b one of the Fourier frequencies is replaced with an off-grid
point $0.5$ rad/sec. Also in this last case ATOM2 yields the lowest MSE
getting closer and closer to the CRLB.

(a)

(b)

Figure 6: MSE(${\mathbf{r}}$) vs. number of samples $n$ for Toepltiz
covariance matrix with condition number constraint. a) on-grid frequencies; b)
off-grid frequencies

###  V-F Radar Application

In this subsection, the performance of the covariance estimation algorithms is
evaluated with reference to the maximum achievable SINR in a radar spatial
processing context. To this end, let us consider a radar system having uniform
linear array with $m=6$ sensors, pointing toward the boresight direction. The
distance between each sensors is equal to $d=\dfrac{\lambda}{2}$, where
$\lambda$ is the radar operating wavelength. The interference covariance
matrix is constructed as
${\mathbf{R}}={\mathbf{R}}_{s}+\sigma_{a}^{2}{\mathbf{I}}$ where
$\sigma_{a}^{2}$ is the power level of the white disturbance noise and
${\mathbf{R}}_{s}$ is obtained by summing the covariance matrices of the $J$
wide-band jammers i.e.,

$\begin{array}[]{ll}{\mathbf{R}}_{s}(p,q)=\displaystyle\sum_{i=1}^{J}\sigma_{i}^{2}\textrm{sinc}(0.5B_{f}(p-q)\phi_{i})e^{{j}(p-q)\phi_{i}}\hskip
28.45274ptp,q=1,2,\cdots m\end{array}$ (54)

where $\sigma_{i}^{2}$ denotes the power of the $i$-th interferer,
$B_{f}=\dfrac{B}{f}$ is the fractional bandwidth, $B$ is the instantaneous
Bandwidth, $f=\dfrac{c}{\lambda}$ and $c$ is the speed of light. The $i^{th}$
jammer phase angle is denoted by $\phi_{i}$ and is equal to
$\phi_{i}=\dfrac{2\pi d\sin(\theta_{i})}{\lambda}$. The covariance estimation
algorithms are compared using the average SINR metric

${\rm{SINR}}_{\rm{avg}}=\dfrac{1}{K}\displaystyle\sum_{i=1}^{K}\dfrac{|\hat{{\boldsymbol{w}}_{i}}^{H}{\boldsymbol{s}}(\theta)|^{2}}{\hat{{\boldsymbol{w}}}_{i}^{H}{\mathbf{R}}\hat{{\boldsymbol{w}}_{i}}}$
(55)

where $K$ is the number of Monte-Carlo trials, $s(\theta)$ is the $m$
dimensional steering vector and
${\hat{{\boldsymbol{w}}}_{i}}={\hat{{\mathbf{R}}}_{i}}^{-1}{\boldsymbol{s}}(\theta)$
is the estimate of the optimal weight vector for spatial processing with
${\hat{{\mathbf{R}}}_{i}}$ being the estimate of the interference covariance
matrix for the $i$-th trial. Previous assumptions imply that for the above
considered linear array, the steering vector is given by
$s(\theta)=[1,e^{j\pi\sin(\theta)},\cdots,e^{j\pi\sin(\theta)(m-1)}]^{T}$. The
simulation setup assumes $J=2$ jammers having the same power
$\sigma_{i}^{2}=20$ dB with $\theta_{1}=9.8^{\circ}$ and
$\theta_{2}=-8.8^{\circ}$. The fractional bandwidth of the two jammers and
white noise power level is equal to 0.3 and $10$ dB, respectively. The average
SINR as a function of theta is shown in Fig. 7.a, Fig. 7.b, and Fig. 7.c for
$n=m$, $n=2m$, and $n=3m$, respectively. These figures also indicate the SINR
bound calculated as
${\boldsymbol{s}}(\theta)^{H}{\mathbf{R}}^{-1}{\boldsymbol{s}}(\theta)$.
Inspection of the plots highlights that as the number of samples $n$ increases
ATOM1 and ATOM2 gets nearer and nearer to the SINR bound.

(a)

(b)

(c)

Figure 7: Average SINR vs $\theta$ in the presence of multiple jammers. a)
$m=n$ b) $n=2m$, and c) $n=3m$

##  VI Conclusion

In this paper, MLE of TSC matrices is considered. Precisely leveraging the MM
framework, two iterative algorithms ATOM1 and ATOM2 are developed. Both
inherit the key properties of MM i.e., they monotonically decrease the
underlying cost function with guaranteed convergence to a stationary point of
the equivalent MLE problem. Subsequently, ATOM2 is extended to handle
covariance matrix MLE forcing other Toeplitz-related structures, such as:
banded Toeplitz, TBT, low rank Toeplitz plus a scalar matrix, and Toeplitz
structure satisfying a condition number constraint. Simulation results
indicate that the proposed algorithms can perform better than some state-of-
the-art techniques in terms of MSE and the SINR metrics.
Some of the possible future research directions are now outlined. In the
present work, for low rank TSC matrix estimation, the rank ($r$) of the
covariance matrix is assumed to be known. However, in many scenarios a perfect
prior knowledge of $r$ could not be available. Therefore, the proposed
algorithm could be extended and analyzed to the case of unknown $r$. Another
possible extension of the proposed technique could be MLE of a Toeplitz
covariance matrix assuming a compound Gaussian distribution for the
underlining data which has a significant application in low-grazing angle
target detection [55, 56].

## Appendix A

## Equivalence of ($P_{1}$) and ($P_{2}$)

Let us cast Problem ($P_{1}$) as

$\begin{array}[]{ll}\underset{{{\mathbf{R}}}\in Toep,{\mathbf{R}}\succ
0,q>0}{\rm minimize}\>q+\log|{\mathbf{R}}|\>\>{\rm
s.t.}\quad\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{{\mathbf{R}}}^{-1}{\mathbf{y}}_{i}\leq
q\end{array}$ (56)

Now, with variable substitution ${\boldsymbol{X}}={\mathbf{R}}q$, Problem (56)
is equivalent to

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}}\in Toep,{\boldsymbol{X}}\succ
0,q>0}{\rm
min}\>q+\log|\dfrac{{\mathbf{I}}}{q}|+\log|{\boldsymbol{X}}|\>\>{\rm
s.t.}\quad\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{{\boldsymbol{X}}}^{-1}{\mathbf{y}}_{i}\leq
1\end{array}$ (57)

The minimizer w.r.t. $q$ is constant and is given as $q^{*}=m$. Next, an
optimal solution w.r.t. ${\boldsymbol{X}}$ can be obtained by solving the
following optimization problem

$\begin{array}[]{ll}{\boldsymbol{X}}^{*}=\underset{{{\boldsymbol{X}}}\in
Toep,{\boldsymbol{X}}\succ 0}{\rm arg\>min}\>\log|{\boldsymbol{X}}|\>\>{\rm
s.t.}\quad\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}{\mathbf{y}}_{i}^{H}{{\boldsymbol{X}}}^{-1}{\mathbf{y}}_{i}\leq
1\end{array}$ ($P_{2}$)

As a result, a minimizer to Problem ($P_{1}$) is
${\mathbf{R}}^{*}=\dfrac{{\boldsymbol{X}}^{*}}{q^{*}}=\dfrac{{\boldsymbol{X}}^{*}}{m}$.

## Appendix B

## SDP Formulation of Problem (10)

Problem (10) can be written as

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}\in Toep,{\boldsymbol{X}}\succ
0}}{\rm
minimize}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)\\\
{\rm subject\
to}\quad\bar{{\mathbf{y}}}^{H}\boldsymbol{\hat{X}}^{-1}\bar{{\mathbf{y}}}\leq{n}\\\
\quad\quad\quad\quad\quad\boldsymbol{\hat{X}}={\mathbf{I}}\otimes{\boldsymbol{X}}\end{array}$
(58)

where $\boldsymbol{\hat{X}}$ is a block diagonal matrix with the diagonal
blocks equal to ${\boldsymbol{X}}$ and $\bar{{\mathbf{y}}}$ is a vector
obtained by stacking the vectors
${\mathbf{y}}_{1},{\mathbf{y}}_{2},\cdots,{\mathbf{y}}_{n}$. Being
${\boldsymbol{X}}\succ 0$, using the Schur complement theorem [28], Problem
(58) can be formulated as an SDP

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}\in
Toep,{\boldsymbol{\hat{X}}}}}{\rm
minimize}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)\\\
{\rm subject\ to}\>\begin{bmatrix}n&\bar{{\mathbf{y}}}^{H}\\\
\bar{{\mathbf{y}}}&\boldsymbol{\hat{X}}\end{bmatrix}\succeq 0\end{array}$ (59)

Once again leveraging the Schur complement theorem, the following optimization
problem is obtained

$\begin{array}[]{ll}\underset{{{\boldsymbol{X}}\in
Toep,{\boldsymbol{\hat{X}}}}}{\rm
minimize}\>\textrm{Tr}\left(({{\boldsymbol{X}}}_{t})^{-1}{{\boldsymbol{X}}}\right)\\\
{\rm subject\
to}\quad\boldsymbol{\hat{X}}\succeq\dfrac{1}{n}\bar{{\mathbf{y}}}\bar{{\mathbf{y}}}^{H}\end{array}$
(60)

which is also an SDP.

## Appendix C

## Proof of Lemma III.3

To begin with, let us denote by $g({\boldsymbol{X}}|{\boldsymbol{X}}_{t})$ the
objective function involved in the surrogate optimization problem exploited by
either ATOM1 or ATOM2. This function, regardless of the method, satisfies the
following two inequalities

$g({\boldsymbol{X}}_{t}|{\boldsymbol{X}}_{t})=\tilde{f}({\boldsymbol{X}}_{t})$
(61) $g({\boldsymbol{X}}_{t+1}|{\boldsymbol{X}}_{t})\geq
g\tilde{f}({\boldsymbol{X}}_{t+1})$ (62)

where $\tilde{f}({\boldsymbol{X}})=\log|{\boldsymbol{X}}|$. Leveraging the
above inequalities, it follows that

$\tilde{f}({\boldsymbol{X}}_{t+1})\overset{(a)}{\leq}g({\boldsymbol{X}}_{t+1}|{\boldsymbol{X}}_{t})\overset{(b)}{\leq}g({\boldsymbol{X}}_{t}|{\boldsymbol{X}}_{t})\overset{(c)}{=}\tilde{f}({\boldsymbol{X}}_{t})$
(63)

In (63), the inequality (a) and equality (c) stem from (62) and (61),
respectively; besides, the inequality (b) is obtained by exploiting the fact
that ATOM1 and ATOM2 globally solve the corresponding convex surrogate
optimization problem. Therefore, (63) implies that the sequence of objective
value of Problem $(P_{2})$ generated by the proposed algorithms is
monotonically decreasing , i.e.,

$\tilde{f}({\boldsymbol{X}}_{0})\geq\tilde{f}({\boldsymbol{X}}_{1})\geq\tilde{f}({\boldsymbol{X}}_{2})\geq\cdots$
(64)

Next, let us denote by ${\boldsymbol{Z}}$ a cluster point to
$\\{{\boldsymbol{X}}_{t}\\}$ and let $\\{{\boldsymbol{X}}_{r_{t}}\\}$ be a
subsequence of $\\{{\boldsymbol{X}}_{t}\\}$ converging to ${\boldsymbol{Z}}$.
Then, from (61), (62), and (64)

$\begin{array}[]{ll}g\left({\boldsymbol{X}}_{r_{t+1}}|{\boldsymbol{X}}_{r_{t+1}}\right)=\tilde{f}\left({\boldsymbol{X}}_{t_{j+1}}\right)\leq\tilde{f}\left({\boldsymbol{X}}_{r_{t}+1}\right)\leq
g\left({\boldsymbol{X}}_{r_{t}+1}|{\boldsymbol{X}}_{r_{t}}\right)\leq
g\left({\boldsymbol{X}}|{\boldsymbol{X}}_{r_{t}}\right),\forall\,\,{\mbox{feasible}}\,\,{\boldsymbol{X}}.\end{array}$
(65)

Thus, letting $t\rightarrow\infty$

$g({\boldsymbol{Z}}|{\boldsymbol{Z}})\leq
g({\boldsymbol{X}}|{\boldsymbol{Z}}),$ (66)

which implies that
$g^{\prime}({\boldsymbol{Z}}|{\boldsymbol{Z}};{\boldsymbol{D}})\geq 0$ where
$g^{\prime}(\cdot|{\boldsymbol{Z}};{\boldsymbol{D}})$ is the directional
derivative of the surrogate function at point ${\boldsymbol{Z}}$ in a feasible
direction ${\boldsymbol{D}}$. Finally, by Proposition 1 in [26], the surrogate
function $g(\cdot|{\boldsymbol{Z}})$ and the objective function
$\tilde{f}(\cdot)$ have the same first order behavior at ${\boldsymbol{Z}}$.
Therefore, $g^{\prime}({\boldsymbol{Z}}|{\boldsymbol{Z}};{\boldsymbol{D}})\geq
0$ implies that $\tilde{f}^{\prime}({\boldsymbol{Z}};{\boldsymbol{D}})\geq 0$.
Hence, ${\boldsymbol{Z}}$ is a stationary point of the objective function
$\tilde{f}(\cdot)$ of $(P_{2})$.

## References

  * [1] D. R. Fuhrmann, E. J. Kelly, and R. Nitzberg, “A cfar adaptive matched filter detector,” _IEEE Trans. Aerosp. Electron. Syst_ , vol. 28, no. 1, pp. 208–216, 1992.
  * [2] P. Stoica, R. L. Moses _et al._ , “Spectral analysis of signals,” 2005.
  * [3] I. S. Reed, J. D. Mallett, and L. E. Brennan, “Rapid convergence rate in adaptive arrays,” _IEEE Transactions on Aerospace and Electronic Systems_ , no. 6, pp. 853–863, 1974.
  * [4] A. Farina, “Antenna-based signal processing techniques for radar systems(book),” _Norwood, MA: Artech House, 1992._ , 1992.
  * [5] B. Armstrong, H. Griffiths, C. Baker, and R. White, “Performance of adaptive optimal doppler processors in heterogeneous clutter,” _IEE Proceedings-Radar, Sonar and Navigation_ , vol. 142, no. 4, pp. 179–190, 1995\.
  * [6] B. Himed and W. L. Melvin, “Analyzing space-time adaptive processors using measured data,” in _Conference Record of the Thirty-First Asilomar Conference on Signals, Systems and Computers_ , vol. 1. IEEE, 1997, pp. 930–935.
  * [7] D. R. Fuhrmann, “Application of toeplitz covariance estimation to adaptive beamforming and detection,” _IEEE Transactions on Signal Processing_ , vol. 39, no. 10, pp. 2194–2198, 1991.
  * [8] W. J. Roberts and Y. Ephraim, “Hidden markov modeling of speech using toeplitz covariance matrices,” _Speech Communication_ , vol. 31, no. 1, pp. 1–14, 2000.
  * [9] C.-I. Chang, “Hyperspectral imaging: techniques for spectral detection and classification,” vol. 1, 2003.
  * [10] C. Caratheodory and L. Fejér, “On the connection between the extremes of harmonic functions and their coefficients and on picard-landau’s theorem,” _Rendiconti del Circolo Matematico di Palermo (1884-1940)_ , vol. 32, no. 1, pp. 218–239, 1911.
  * [11] U. Grenander and G. Szegö, “Toeplitz forms and their applications,” 1958.
  * [12] V. F. Pisarenko, “The retrieval of harmonics from a covariance function,” _Geophysical Journal International_ , vol. 33, no. 3, pp. 347–366, 1973.
  * [13] A. Dembo, C. L. Mallows, and L. A. Shepp, “Embedding nonnegative definite toeplitz matrices in nonnegative definite circulant matrices, with application to covariance estimation,” _IEEE Transactions on Information Theory_ , vol. 35, no. 6, pp. 1206–1212, 1989.
  * [14] S. Haykin, “Nonlinear methods of spectral analysis,” vol. 34, 2006.
  * [15] M. I. Miller and D. L. Snyder, “The role of likelihood and entropy in incomplete-data problems: Applications to estimating point-process intensities and toeplitz constrained covariances,” _Proceedings of the IEEE_ , vol. 75, no. 7, pp. 892–907, 1987.
  * [16] M. J. Turmon and M. I. Miller, “Maximum-likelihood estimation of complex sinusoids and toeplitz covariances,” _IEEE transactions on signal processing_ , vol. 42, no. 5, pp. 1074–1086, 1994.
  * [17] L. P. Christensen, “An em-algorithm for band-toeplitz covariance matrix estimation,” in _2007 IEEE International Conference on Acoustics, Speech and Signal Processing-ICASSP’07_ , vol. 3. IEEE, 2007, pp. III–1021.
  * [18] P. Babu, “Melt—maximum-likelihood estimation of low-rank toeplitz covariance matrix,” _IEEE Signal Processing Letters_ , vol. 23, no. 11, pp. 1587–1591, 2016.
  * [19] H. Li, P. Stoica, and J. Li, “Computationally efficient maximum likelihood estimation of structured covariance matrices,” _IEEE Transactions on Signal Processing_ , vol. 47, no. 5, pp. 1314–1323, 1999.
  * [20] B. Kang, V. Monga, and M. Rangaswamy, “Computationally efficient toeplitz approximation of structured covariance under a rank constraint,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 51, no. 1, pp. 775–785, 2015.
  * [21] X. Du, A. Aubry, A. De Maio, and G. Cui, “Toeplitz structured covariance matrix estimation for radar applications,” _IEEE Signal Processing Letters_ , vol. 27, pp. 595–599, 2020.
  * [22] F. Barbaresco, “Information geometry of covariance matrix: Cartan-siegel homogeneous bounded domains, mostow/berger fibration and frechet median,” in _Matrix Information Geometry_. Springer, 2013, pp. 199–255.
  * [23] B. Balaji, F. Barbaresco, and A. Decurninge, “Information geometry and estimation of toeplitz covariance matrices,” in _2014 International Radar Conference_. IEEE, 2014, pp. 1–4.
  * [24] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 3, pp. 794–816, 2016.
  * [25] D. R. Hunter and K. Lange, “A tutorial on mm algorithms,” _The American Statistician_ , vol. 58, no. 1, pp. 30–37, 2004.
  * [26] M. Razaviyayn, M. Hong, and Z.-Q. Luo, “A unified convergence analysis of block successive minimization methods for nonsmooth optimization,” _SIAM Journal on Optimization_ , vol. 23, no. 2, pp. 1126–1153, 2013.
  * [27] A. Aubry, A. De Maio, A. Zappone, M. Razaviyayn, and Z.-Q. Luo, “A new sequential optimization procedure and its applications to resource allocation for wireless systems,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 24, pp. 6518–6533, 2018.
  * [28] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [29] Y. Nesterov, “Introductory lectures on convex optimization: A basic course,” _Springer Science & Business Media_, vol. 87, 2003.
  * [30] F. Xu and P. Pan, “A new algorithm for positive semidefinite matrix completion,” _Journal of Applied Mathematics_ , vol. 2016, 2016.
  * [31] M. Fazel, T. K. Pong, D. Sun, and P. Tseng, “Hankel matrix rank minimization with applications to system identification and realization,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 34, no. 3, pp. 946–977, 2013\.
  * [32] S. Boyd, N. Parikh, and E. Chu, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Now Publishers Inc_ , 2011.
  * [33] R. Nishihara, L. Lessard, B. Recht, A. Packard, and M. Jordan, “A general analysis of the convergence of admm,” in _International Conference on Machine Learning_. PMLR, 2015, pp. 343–352.
  * [34] G. H. Golub and C. F. Van Loan, “Matrix computations,” 1996.
  * [35] Z. Wen, D. Goldfarb, and W. Yin, “Alternating direction augmented lagrangian methods for semidefinite programming,” _Mathematical Programming Computation_ , vol. 2, no. 3-4, pp. 203–230, 2010.
  * [36] C. Chen, B. He, Y. Ye, and X. Yuan, “The direct extension of admm for multi-block convex minimization problems is not necessarily convergent,” _Mathematical Programming_ , vol. 155, no. 1-2, pp. 57–79, 2016.
  * [37] H. H. Bauschke and J. M. Borwein, “On projection algorithms for solving convex feasibility problems,” _SIAM review_ , vol. 38, no. 3, pp. 367–426, 1996\.
  * [38] W. Cheney and A. A. Goldstein, “Proximity maps for convex sets,” _Proceedings of the American Mathematical Society_ , vol. 10, no. 3, pp. 448–450, 1959.
  * [39] C. L. Byrne, “Alternating minimization and alternating projection algorithms: A tutorial,” _Sciences New York_ , pp. 1–41, 2011.
  * [40] D. Henrion and J. Malick, “Projection methods for conic feasibility problems: applications to polynomial sum-of-squares decompositions,” _Optimization Methods & Software_, vol. 26, no. 1, pp. 23–46, 2011.
  * [41] J. P. Boyle and R. L. Dykstra, “A method for finding projections onto the intersection of convex sets in hilbert spaces,” in _Advances in order restricted statistical inference_. Springer, 1986, pp. 28–47.
  * [42] R. L. Dykstra, “An algorithm for restricted least squares regression,” _Journal of the American Statistical Association_ , vol. 78, no. 384, pp. 837–842, 1983.
  * [43] O. E. Livne and A. Brandt, “$n$ roots of the secular equation in $\mathcal{O}(n)$ operations,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 24, no. 2, pp. 439–453, 2002.
  * [44] F. T. Luk and S. Qiao, “A fast eigenvalue algorithm for hankel matrices,” _Linear Algebra and Its Applications_ , vol. 316, no. 1-3, pp. 171–182, 2000\.
  * [45] A. P. Shikhaliev, L. C. Potter, and Y. Chi, “Low-rank structured covariance matrix estimation,” _IEEE Signal Processing Letters_ , vol. 26, no. 5, pp. 700–704, 2019.
  * [46] Y. Sun, A. Breloy, P. Babu, D. P. Palomar, F. Pascal, and G. Ginolhac, “Low-complexity algorithms for low rank clutter parameters estimation in radar systems,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 8, pp. 1986–1998, 2015.
  * [47] J. A. Tropp, I. S. Dhillon, R. W. Heath, and T. Strohmer, “Designing structured tight frames via an alternating projection method,” _IEEE Transactions on information theory_ , vol. 51, no. 1, pp. 188–209, 2005.
  * [48] T. A. Barton and S. T. Smith, “Structured covariance estimation for space-time adaptive processing,” in _1997 IEEE International Conference on Acoustics, Speech, and Signal Processing_ , vol. 5. IEEE, 1997, pp. 3493–3496.
  * [49] ——, “Structured covariance estimation for space-time adaptive processing,” in _1997 IEEE International Conference on Acoustics, Speech, and Signal Processing_ , vol. 5. IEEE, 1997, pp. 3493–3496.
  * [50] Y. I. Abramovich, B. A. Johnson, and N. K. Spencer, “Two-dimensional multivariate parametric models for radar applications—part i: Maximum-entropy extensions for toeplitz-block matrices,” _IEEE transactions on signal processing_ , vol. 56, no. 11, pp. 5509–5526, 2008.
  * [51] ——, “Two-dimensional multivariate parametric models for radar applications—part ii: Maximum-entropy extensions for hermitian-block matrices,” _IEEE transactions on signal processing_ , vol. 56, no. 11, pp. 5527–5539, 2008.
  * [52] A. Aubry, A. De Maio, and L. Pallotta, “A geometric approach to covariance matrix estimation and its applications to radar problems,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 4, pp. 907–922, 2017.
  * [53] L. Pallotta, A. Aubry, A. De Maio, and A. Farina, “Estimation of a structured covariance matrix with a condition number constraint for radar applications,” in _2012 IEEE Radar Conference_. IEEE, 2012, pp. 0778–0783.
  * [54] S. M. Kay, _Fundamentals of statistical signal processing: estimation theory_. Prentice-Hall, Inc., 1993.
  * [55] K. Ward, C. Baker, and S. Watts, “Maritime surveillance radar. part 1: Radar scattering from the ocean surface,” in _IEE Proceedings_ , vol. 137, no. 2. IET, 1990, pp. 51–62.
  * [56] F. Gini and M. Greco, “Suboptimum approach to adaptive coherent radar detection in compound-gaussian clutter,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 35, no. 3, pp. 1095–1104, 1999.