Convex-structured covariance estimation via the entropy loss under the majorization-minimization algorithm framework

1.   Introduction

Structural covariance/correlation matrices appear in diverse fields. For example, in array signal processing, the covariance matrix of received array signal is noise covariance plus a weighted sum of rank-one matrices ^[1,2,3]; in time-series analysis, a common form is the Toeplitz structural autocovariance/autocorrelation matrix ^[4]; and in certain high-dimensional problems, the covariance matrix can often be sparse or exhibit a banded structure ^[5]. Covariance estimation using structure constraints is common in research articles, e.g., ^[5,6,7,8]. It is believed that structural estimates are more accurate than that excluding structural information because exploiting structure information in the estimation process usually implies a reduction in the number of parameters to be estimated. Structural covariance/correlation estimation plays an important role in several applications, such as low-rank signal detection ^[9,10], direction-of-arrival (DOA) estimation ^[2,3], and spectral analysis of time-series data ^[4].

The Weighted sum of known Rank-one matrices with unknown Weights (W-Rank1-W) structural matrices are often used to describe the covariance of received array signal ^[2]. In ^[11], a weighted covariance fitting criterion was proposed, providing a W-Rank1-W structural covariance estimation of the received array signals. Based on the maximum likelihood principle, ^[12] estimated the W-Rank1-W structural covariance of the received array signals that were assumed to follow a circular Gaussian distribution. The weight estimates in ^[11,12] are widely known to play an important role in DOA estimation. Subsequently, many researchers have noted that the weights of W-Rank1-W are sparse in DOA estimation, and similar DOA estimation ideas were implemented ^[2,13], where the W-Rank1-W structural covariance estimation was addressed specifically under sparse weight constraints. We estimate the W-Rank1-W structural covariance by minimizing the relative entropy (also known as Kullback-Leibler divergence) between two circular Gaussian distributions ^[14], yielding a low-complexity algorithm based on the Majorization-Minimization (MM) algorithm framework ^[15]. The entropy loss measures the difference between pairs of distributions with two different covariance matrices. The proposed MM algorithm is based on the result that the entropy loss function is bounded above by a sequence of separable functions. Our numerical experiments reveal that weight estimates based on minimizing the entropy loss are more accurate for DOA estimation than the estimates reported in ^[11,12] using the weighted covariance fitting criteria and the maximum-likelihood principle.

An auto-covariance matrix in a stationary time series always exhibits a Toeplitz structure ^[4]. Moreover, in certain cases, the autocorrelation coefficient is inversely related to the time lag, e.g., in auto-regressive time series, some moving average time series, and some auto-regressive moving average time series. ^[16] approached Toeplitz covariance estimation by minimizing the F-norm loss of the sample covariance matrix, and ^[6,17] proposed maximum likelihood Toeplitz estimates under circularly-symmetric complex normal distribution and complex elliptical distribution assumptions. ^[18] provided a banded Toeplitz estimator by projecting a given positive definite matrix onto the convex cone of nonnegative definite banded Toeplitz matrices. ^[19,20] proposed Toeplitz estimates by minimizing entropy loss and modified entropy loss, respectively. ^[21] proposed Toeplitz estimation based on entropy loss, considering the possible sparsity of the covariance matrix. This study considers the inverse relationship between the autocorrelation coefficient and time lag in time-series analysis and proposes a Toeplitz autocorrelation matrix estimation by minimizing the entropy loss between two normal distributions. As Toeplitz matrices can be expressed in the W-Rank1-W structure, an MM algorithm similar to that for the W-Rank1-W structural estimation is proposed; however, in each iteration of the algorithm, a second-order cone programming (SOCP) problem is required to be solved. Numerical experiments reveal that our estimates are more accurate than those that do not consider the inverse relationship between the autocorrelation coefficient and time lag.

The remainder of this paper is organized as follows: In Section 2, we formulate the estimation problem of interest. Section 3 presents two different algorithms for the W-Rank1-W structural covariance estimation and Section 4 presents a MM algorithm for the autocorrelation matrix estimation involving the inverse relationship between the autocorrelation coefficient and time lag. Finally, Section 5 describes the results of our numerical experiments, revealing the characteristics and applications of the proposed algorithms.

Notations: $\mathbb{C}^{ {{{v}}}} (\mathbb{R}^{ {{{v}}}})$ denotes the set of ${v}$-dimensional vectors in complex (real) field. $\mathbb{S}^ {v}_{+} (\mathbb{S}^ {v}_{++})$ denotes the set of ${v} \times {v}$ symmetric (real field) and hermitian (complex field) positive semi-definite (definite) matrices. ${\mathbf{0}}$ denotes the zero matrix and ${\mathbf{I}}_{ {v}}$ denotes the ${v} \times {v}$ identity matrix. ${\mathbf{X}}_1 - {\mathbf{X}}_2 \succeq 0 \; (\succ 0)$ indicates that the matrix ${\mathbf{X}}_1 - {\mathbf{X}}_2$ is positive semi-definite (definite). ${\mathbf{x}} \succeq 0 (\succ 0)$ indicates that all elements in the vector ${\mathbf{x}}$ are non-negative (positive). ${\mathbf{X}}_{i, j}$ denotes the (i, j)-th entry of the matrix ${\mathbf{X}}$, where the row and column indices begin from $0$. ${\rm{diag}}({\mathbf{x}})$ denotes a diagonal matrix with the vector ${\mathbf{x}}$ defining its diagonal elements. The superscripts $(\cdot)^{ \mathrm{T}}$ and $(\cdot)^{ \mathrm{H}}$ represent the transpose and conjugate transpose matrices, respectively. $\|\cdot\|$ denotes the $l_2$-norm of a vector, and $|\cdot|$ denotes the determinant of a scale.

2.   Problem formulation

Let $\mathcal{CN}_1$ and $\mathcal{CN}_2$ denote circularly-symmetric complex normal distributions with zero mean and ${\bf{\Sigma}}_1 \in \mathbb{S}_{++}^{ {{{v}}}}$ and ${\bf{\Sigma}}_2 \in \mathbb{S}_{++}^{ {{{v}}}}$ as covariance matrices, respectively. The probability densities of $\mathcal{CN}_1$ and $\mathcal{CN}_2$ are

The relative entropy from $\mathcal{CN}_1$ to $\mathcal{CN}_2$ is defined ^[14] as follows:

where $\mathrm{E}_{ {\mathbf{z}} \sim \mathcal{CN}_1}(\cdot)$ denotes the expectation operator, given that ${\mathbf{z}}$ follows the distribution, $\mathcal{CN}_1$. By a brief proof (Appendix A), we have (2.2) rewritten as

This relative entropy measures the difference between the distributions, $\mathcal{CN}_1$ and $\mathcal{CN}_2$.

In applications such as array signal processing and time-series analysis, the covariance is structural, and structural covariance estimation is a common problem. Let ${\bf{\Sigma}}_1 \triangleq {\bf{\Sigma}}$ be an unknown structural covariance matrix and ${\bf{\Sigma}}_2 \triangleq {\mathbf{S}}$ be a given covariance matrix. By minimizing the entropy loss, $D(\mathcal{CN}_2 \| \mathcal{CN}_1)$,

where $\Omega$ denotes a nonempty set that is the intersection of the closed set that characterizes covariance structure and a positive semi-definite cone $\mathbb{S}_{+}^{ {{{v}}}}$, we achieve structural covariance estimation.

Throughout this study, we adopt the following assumptions: the loss function $L_0({\bf{\Sigma}}_k) \rightarrow +\infty$ when the sequence $\{ {\bf{\Sigma}}_k\}$ tends towards the boundary of the positive semi-definite cone, $\mathbb{S}^{ {{{v}}}}_{+}$. Under this assumption, the case in which ${\bf{\Sigma}}$ is singular can be excluded from the algorithm analysis.

3.   W-Rank1-W structural covariance matrix estimation

Estimation of the W-Rank1-W structural covariance matrix is described in this section, i.e., the estimation of ${\bf{\Sigma}}$ in the set $\Omega_1$:

where ${\mathbf{a}}_0, \ldots, {\mathbf{a}}_m$ denote the vectors known in advance. W-Rank1-W structural covariance commonly occurs in array signal processing and can often be applied to DOA estimation.

Given a covariance matrix, which is often the sample covariance matrix, we solve the problem (2.4) with $\Omega = \Omega_1$:

yielding a W-Rank1-W structural covariance estimate. Any ${\bf{\Sigma}} \in \Omega_1$ can be rewritten as:

where ${\mathbf{A}} = ({\mathbf{a}}_0, {\mathbf{a}}_1, \cdots, {\mathbf{a}}_m)$ and ${\mathbf{p}} = (p_0, p_1, \cdots, p_m)^{ \mathrm{T}}$. Because ${\bf{\Sigma}}$ is a function of ${\mathbf{p}}$ (${\mathbf{A}}$ is known in advance), the objective function in (3.2) can be rewritten as the following function of ${\mathbf{p}}$:

Using the objective function expression, $L({\mathbf{p}})$, the optimization problem (3.2) can be rewritten as follows:

The objective function in (3.2) is a convex function in terms of ${\bf{\Sigma}}$ ^[22], and ${\bf{\Sigma}}$ is an affine map of ${\mathbf{p}}$; therefore, it is convex in ${\mathbf{p}}$, i.e., $L({\mathbf{p}})$ is convex in ${\mathbf{p}}$. Moreover, the constraint set $\{ {\mathbf{p}} \mid {\mathbf{p}} \succeq 0 \}$ is convex. Thus, (3.5) represents a convex optimization problem.

The following algorithms for W-Rank1-W covariance estimation are used to solve problem (3.5). When ${\mathbf{p}}^*$ denotes the minimum point of the problem (3.5):

the W-Rank1-W structural covariance estimation is given by

We present two algorithms in the following subsections—the primal-dual interior-point method and another designed within the MM algorithm framework. These two algorithms have different computational complexities, and we briefly discuss and compare them.

3.1.   Primal-dual interior-point algorithm

Problem (3.5) is a convex optimization problem with linear inequality constraints that can be solved effectively using the primal-dual interior-point method ^[23].

Algorithm 1 outlines the primal-dual interior-point method in the context of (3.5). In the algorithm, the primal-dual search direction $(\Delta {\mathbf{p}}, \Delta \mathit{\boldsymbol{\lambda}})$ at the point $({\mathbf{p}}, \mathit{\boldsymbol{\lambda}})$ is defined as the solution of the linear equation:

where $\tau > 0$ denotes a preset parameter, and ${\mathbf{h}} = (h_i)$ and ${\mathbf{H}} = (H_{ij})$ with $i, j = 0, 1, \ldots, m$ denote the gradient and Hessian matrix of $L({\mathbf{p}})$ at the point $({\mathbf{p}}, \mathit{\boldsymbol{\lambda}})$, respectively. Specifically,

3.2.   MM algorithm

In each iteration of the primal-dual interior-point method given by Algorithm 1, the Hessian matrix computation has computational complexity = $O(m^5)$. This section presents a novel method for decreasing the computational complexity of each iteration using the MM algorithm framework ^[15].

In the MM algorithm framework, if a continuously differentiable function $g({\mathbf{p}} \mid {\mathbf{p}}_{k})$ satisfies the following condition:

where equality is achieved at ${\mathbf{p}} = {\mathbf{p}}_{k}$, the sequence $\{ {\mathbf{p}}_k\}$ generated by

converges to a stationary point of problem (3.5). Because the optimization problem (3.5) is convex, its stationary point is actually the optimal point, ${\mathbf{p}}^*$.

$\log \det ({\mathbf{X}})$ is a concave function that is bounded above by its first-order Taylor expansion at any ${\mathbf{X}}_k$ ^[15]:

with the equality achieved at ${\mathbf{X}} = {\mathbf{X}}_k$. Let ${\mathbf{X}} \triangleq ({\mathbf{A}} {\rm{diag}}({\mathbf{p}}) {\mathbf{A}}^{ \mathrm{H}})^{-1}$ and ${\mathbf{X}}_k \triangleq ({\mathbf{A}} {\rm{diag}}({\mathbf{p}}_k) {\mathbf{A}}^{ \mathrm{H}})^{-1}$, and then inserting them into (3.13) achieves

where the equality in (3.14b) is achieved at ${\mathbf{p}} = {\mathbf{p}}_k$ and $s_0$ is a constant. Applying the Schur complement (a brief proof is provided in Appendix B), for any ${\mathbf{p}}_k \succ 0$, we obtain

where

and equality is achieved at ${\mathbf{p}} = {\mathbf{p}}_k$. By substituting (3.15) into (3.14b), we obtain

where equality is achieved at ${\mathbf{p}} = {\mathbf{p}}_k$. Because $\log \det\left(({\mathbf{A}} {\rm{diag}}({\mathbf{p}}) {\mathbf{A}}^{ \mathrm{H}})^{-1}\right) = - \log \det\left(({\mathbf{A}} {\rm{diag}}({\mathbf{p}}) {\mathbf{A}}^{ \mathrm{H}})\right)$, then by (3.17) we have

where equality is achieved at ${\mathbf{p}} = {\mathbf{p}}_k$. Let us denote

by

By substituting (3.18) into (3.4), for any ${\mathbf{p}}_k \succ 0$, we obtain

with equality achieved at ${\mathbf{p}} = {\mathbf{p}}_k$.

By applying (3.21) under the MM algorithm framework, we derive the following proposition:

Proposition 3.1. For any value of ${\mathbf{p}}_0 \succ 0$: sequentially solving

generates a sequence $\{ {\mathbf{p}}_k\}$ that converges to the global optimal point ${\mathbf{p}}^*$ of the problem (3.5).

Proof. Because (3.21) holds only for any ${\mathbf{p}}_k \succ 0$, and not for any ${\mathbf{p}}_k \succeq 0$ as necessary in (3.11), the general convergence result of MM algorithms cannot be applied here. Therefore, we need to analyze the convergence further.

Consider the following $\epsilon$-approximation of the problem (3.5):

with $\epsilon > 0$.

Now, applying (3.22) to $\tilde{ {\mathbf{p}}} \triangleq {\mathbf{p}} + \epsilon \mathit{\boldsymbol{1}} \succ 0$, we conclude that the sequence $\{ {\mathbf{p}}_k^{\epsilon}\}$ converges to the optimal point $({\mathbf{p}}^{\epsilon})^*$ of (3.23), and

for any feasible direction ${\mathbf{d}}$, where $\nabla L^{\epsilon}(({\mathbf{p}}^{\epsilon})^*)$ denotes a gradient of the objective function $L^{\epsilon}({\mathbf{p}})$ for ${\mathbf{p}}$ in (3.23) at $({\mathbf{p}}^{\epsilon})^*$.

Let $\epsilon_{k'}$ be a positive sequence, with $\lim_{k' \rightarrow +\infty} \epsilon_{k'} = 0$. Then, as $\nabla L^{\epsilon}(({\mathbf{p}}^{\epsilon})^*)$ is continuous around $({\mathbf{p}}^{\epsilon})^*$ and $\epsilon$, the limit point, ${\mathbf{p}}^*$, of the sequence, $\{({\mathbf{p}}^{\epsilon_{k'}})^*\}$, is the optimal point for problem (3.5).

In practice, $\epsilon$ can be chosen as to be an arbitrarily small number. Directly applying (3.22) or adapting it to solve the $\epsilon$-approximation problem is essentially idetnical. Therefore, the sequence $\{ {\mathbf{p}}_k\}$ in (3.22) converges to the optimal point of problem (3.5). □

Problem (3.22), which generates $\{ {\mathbf{p}}_k\}$, is a simple optimization problem. Thus, obtaining an analytical solution to it is easy. We denote the $i$-th element of ${\mathbf{p}}_{k+1}$ by $({\mathbf{p}}_{k+1})_i$. We have

where $({\mathbf{W}}_k)_{i, i}$ and ${\mathbf{M}}_{i, i}$ denote the $i$-th diagonal elements of ${\mathbf{W}}_k$ in (3.16) and ${\mathbf{M}}$ in (3.20), respectively.

We outline the MM algorithm for (3.5) in Algorithm 2.

The computational complexity of each iteration is $O(m^3)$, which is mainly generated by matrix inversion and multiplication. To demonstrate the difference in terms of computational complexity between Algorithms 1 and 2, we compare their iteration numbers and computation times in the numerical experiments presented in Section 5.

4.   Real Toeplitz correlation matrix estimation

This section describes the estimation of Toeplitz matrices. Complex Toeplitz covariance matrices appear commonly in linear array signal processing ^[17,24]. During the application of sparse signal models, the complex Toeplitz covariance matrices can always be formulated as W-Rank1-W structural matrices that have been discussed in Section 3 ^[2]. Thus, in this section, we consider only real Toeplitz matrices—specifically the Toeplitz auto-correlation matrix with decreasing correlation coefficients with respect to increasing time lags in the time series analysis. Note that, the mentioned auto-correlation matrix is composed of auto-correlation coefficients at different lags.

Let us denote

and define ${\rm{toep}}({\mathbf{c}})$ to be a symmetric Toeplitz matrix with ${\mathbf{c}}$ as its first column:

This section estimates the Toeplitz correlation matrix ${\mathbf{C}} \in \Omega_2$, where

Note that the constraint $c_0 = 1$ can be rewritten as follows:

where ${\mathbf{e}}$ denotes a row vector with $1$ as its first element, $0$ as its other elements. The constraint $|c_0| \geq |c_1| \geq \cdots \geq |c_{ {v}-1}|$ can be rewritten as linear constraints:

where ${\mathbf{E}}_1$ denotes a $({v}-1) \times {v}$ matrix with diagonal elements = $-1$ and the first subdiagonal elements = $1$, and ${\mathbf{E}}_2$ denotes a $({v}-1) \times {v}$ matrix with diagonal elements = $-1$ and the first subdiagonal elements = $-1$:

We denote a normal distribution with correlation matrix ${\mathbf{R}} \in \mathbb{S}_{++}^{ {{{v}}}}$ by $\mathcal{N}_0$ and a normal distribution with Toeplitz correlation matrix ${\mathbf{C}} \in \Omega_2$ by $\mathcal{N}$. Assuming the variances to be equal, the relative entropy from $\mathcal{N}$ to $\mathcal{N}_0$ is:

Given a correlation matrix estimate of ${\mathbf{R}}$, such as the sample correlation matrix, minimizing $D(\mathcal{N}_0 \| \mathcal{N})$ in terms of ${\mathbf{C}}$:

we deduce the Toeplitz correlation estimate ${\mathbf{C}}^* \in \Omega_2$. The optimization problem (4.8) is convex and can be resolved using the primal-dual interior-point algorithm with the SDT3 toolbox in Matlab. This problem involves a nonlinear objective function and a positive-definiteness constraint of ${\mathbf{C}}$. Application of the primal-dual interior-point algorithm requires the computation of the Hessian matrix and the positive-definiteness to be ensured in each iteration. In the following section, we present an estimation algorithm that avoids these time-consuming computations.

4.1.   Estimation of the real Toeplitz correlation matrix in $\Omega_2$

In this subsection, we first transform the estimate of the real Toeplitz correlation matrix to an estimate of the W-Rank1-W structural matrix by special variable substitution and then deduce an estimation algorithm within the MM algorithm framework for this estimation.

Given a circulant matrix ${\mathbf{T}} = {\rm{toep}}({\mathbf{t}})$ with

and $d_i = d_{l+1-i}$, $i = 1, \cdots, l$ (specifically, ${\mathbf{t}} = (c_0, c_1, \cdots, c_{ {v}-1}, c_{ {v}-1}, c_{ {v}-2}, \cdots, c_1)^{ \mathrm{T}}$ if $l = 0$), i.e.,

where

Then, the upper-left ${v}$-order matrix of ${\mathbf{T}}$ is the Toeplitz matrix:

The circulant matrix ${\mathbf{T}}$ has a Fourier feature decomposition ^[25]:

where ${\mathbf{F}}$ is $m_l = l+2 {v}-1$ order Fourier matrix and

Specifically, ${\mathbf{F}}_{i, j} = \frac{1}{\sqrt{m_l}} \mathrm{e}^{ij \frac{2 \pi \mathrm{\iota}}{m_l}}$ with $\iota$ as the imaginary unit, $i, j = 0, 1, \cdots, m_l-1$, and ${\mathbf{T}} \succeq 0$ if and only if ${\mathbf{p}}_{(l)} \succeq 0$.

By substituting (4.14) into (4.13), we derive

where ${\mathbf{A}}_l = ({\mathbf{I}}_{p}, \mathit{\boldsymbol{0}}) {\mathbf{F}}^{ \mathrm{H}}$.

The Corollary C3 in ^[25] states that

where ${\mathbf{p}}_{(l)}$ is a function of ${\mathbf{t}}$ in (4.15) and ${\mathbf{t}}$ is a function of ${\mathbf{c}}$ in (4.9). Moreover, the item 2 of the Notes in ^[25] states that ${\rm{toep}}({\mathbf{c}}) \succ 0$ implies ${\rm{toep}}({\mathbf{c}}) \in \{ {\mathbf{C}} \mid {\mathbf{C}} = {\rm{toep}}({\mathbf{c}}), {\mathbf{p}}_{(l)} \succeq 0\}$ for a sufficiently large $l$, meaning that if $l$ is sufficiently large,

In the following, we simplify the constraint ${\rm{toep}}({\mathbf{c}}) \succ 0$ in $\Omega_2$ to ${\mathbf{p}}_{(l)} \succ 0$ for certain $l$. If $l$ is sufficiently large, compared to $\{ {\mathbf{C}} \mid {\mathbf{C}} = {\rm{toep}}({\mathbf{c}}), {\rm{toep}}({\mathbf{c}}) \succ 0\}$ containing only positive-definite matrices, by (4.18), the set $\{ {\mathbf{C}} \mid {\mathbf{C}} = {\rm{toep}}({\mathbf{c}}), {\mathbf{p}}_{(l)} \succeq 0\}$ contains more positive semi-definite matrices. Moreover, as ${\mathbf{C}}$ tends toward the boundary of the positive semi-definite cone, the value of the objective function in (4.8) tends to $+\infty$. Thus, although we change the constraint set, if $l$ is sufficiently large, the modification ensures that the optimal point remains constant. The numerical experiments reported in Section 5 demonstrate that the estimation accuracy remains good for different $l$, and we outline some numerical results there.

Now, we present an algorithm for the problem with the constraint ${\rm{toep}}({\mathbf{c}}) \succ 0$ modified to ${\mathbf{p}}_{(l)} \succeq 0$ for certain $l$:

where

In particular, by applying the expressions (4.4), (4.5), (4.10), (4.11), (4.15) and (4.16), the problem to be resolved, i.e., problem (4.19), becomes

It is easy to see ${\mathbf{F}}^{ \mathrm{H}} {\mathbf{F}} = {\mathbf{I}}$ and (4.21c) can be rewritten as follows:

By substituting (4.22) into (4.21e), we obtain (4.21e) as follows:

Now, we apply the equality constraints (4.21c) and (4.21e), i.e., (4.22) and (4.23), and transform the optimization problem (4.21) into a new optimization problem involving only the variable, ${\mathbf{p}}_{(l)}$.

It is easy to see that ${\mathbf{I}}_{(-1)} {\mathbf{F}}^{ \mathrm{H}} = {\mathbf{F}}^{ \mathrm{H}} {\mathbf{I}}_{(-1)}$. Then, by substituting (4.22) into (4.21d), we obtain the constraint (4.21d) as follows:

By substituting (4.23) into (4.21f)–(4.21h), we obtain the constraints (4.21f)–(4.21h) as follows:

Most importantly, applying (4.24)–(4.27), the optimization problem can be rewritten as a problem involving only the variable, ${\mathbf{p}}_{(l)}$, i.e., as an estimation problem of W-Rank-W structural matrices with an unknown weight vector, ${\mathbf{p}}_{(l)}$:

Denoting the optimal point of (4.28) by ${\mathbf{p}}_{(l)}^*$, we have ${\mathbf{C}}^* = {\mathbf{A}}_l {\rm{diag}}({\mathbf{p}}_{(l)}^*) {\mathbf{A}}_l^{ \mathrm{H}}$ as an estimate of the correlation matrix.

In the following subsection, we present an MM algorithm for problem (4.28) by applying the design concept of Algorithm 2.

4.1.1.   MM algorithm for the real Toeplitz correlation matrix estimation

By substituting (4.28b) into the objective function in (4.28a), we obtain: Applying (3.18) to the objective function (4.28a), we obtain the following result. For any $({\mathbf{p}}_{(l)})_k \succ 0$, we have

where $s_1$ is a constant and

and

and

and the equality in (4.29) is achieved at ${\mathbf{p}}_{(l)} = ({\mathbf{p}}_{(l)})_k$.

Applying (4.30) and convergence analysis similar to that for Proposition 3.1, we have the following proposition under the MM algorithm framework.

Proposition 4.1. Given $({\mathbf{p}}_{(l)})_0 \succ 0$, sequentially solving the optimization problem

generates a sequence $\{({\mathbf{p}}_{(l)})_k\}$ that converges to the optimal point ${\mathbf{p}}_{(l)}^*$ of the problem (4.28).

Denoting $r_i = 1 / q_i$ and using the epigraph form, we can rewrite problem (4.35) as a SOCP problem:

This SOCP problem can be solved using the Matlab toolbox, SDPT3.

In conclusion, the MM algorithm for problem (4.28) is presented as Algorithm 3, and is referred to as the MM-T algorithm in the following. In each iteration, Algorithm 3 involves SOCP optimization that can be resolved using the Matlab Toolbox, SDP3.

5.   Numerical experiments

In this section, we evaluate the performance of Algorithms 1–3 numerically. First, Algorithms 1 and 2 are utilized for DOA estimation. Then, Algorithm 3 is used to estimate Toeplitz autocorrelation matrices in time-series analysis and its performance is compared with those of other estimation methods. All computations were performed using Matlab2018a on a system equipped with an Intel(R) Xeon(R) Platinum 8372HC CPU at 3.40GHz.

5.1.   DOA estimation

The observation ${\mathbf{x}}(t) \in \mathbb{C}^{ {{{v}}}}$ received by a uniform linear array at time $t$ is always modeled as ^[2]:

where $\theta_i = (i/2)^{\circ}$ denotes the grid point in $180$-degree coverage: $0^{\circ} \backsim 180^{\circ}$, $s_i(t) \in \mathbb{C}$ denotes the source signal from direction $\theta_i$, ${\mathbf{a}}(\theta_i) = (1, \mathrm{e}^{\imath \pi \cos(\theta_i)}, \cdots, \mathrm{e}^{({v}-1) \imath \pi \cos(\theta_i)})^{ \mathrm{T}}$ denotes the steering vector, and ${\mathbf{n}}(t) \in \mathbb{C}^{ {{{v}}}}$ denotes the noise vector. In addition, let us assume that ${\mathbf{s}}(t) = (s_i(t))$ is $\mathcal{CN}$-distributed with zero mean and covariance ${\rm{diag}}({\mathbf{p}}_{{d}})$, ${\mathbf{n}}(t)$ is $\mathcal{CN}$-distributed with zero mean and covariance ${\rm{diag}}(\mathit{\boldsymbol{\sigma}})$, and ${\mathbf{s}}(t)$ and ${\mathbf{n}}(t)$ are both temporally and mutually independent. In particular, the noise energy vector is $\mathit{\boldsymbol{\sigma}} = \sigma^2 \mathit{\boldsymbol{1}}$ with $\mathit{\boldsymbol{1}}$ as a vector with all elements being $1$, and $\sigma$ can be designed to vary the signal-to-noise ratio (SNR) as follows:

where ${\mathbf{A}}_d = ({\mathbf{a}}(\theta_0), {\mathbf{a}}(\theta_1), \cdots, {\mathbf{a}}(\theta_{360}))$. The covariance matrix of ${\mathbf{x}}(t)$ is given by the W-Rank1-W structure:

where ${\mathbf{A}} = [ {\mathbf{A}}_d, {\mathbf{I}}]$ and ${\mathbf{p}} = (({\mathbf{p}}_{{d}})^{ \mathrm{T}}, \mathit{\boldsymbol{\sigma}}^{ \mathrm{T}})^{ \mathrm{T}}$. We denote ${\mathbf{p}}_{{d}} = (p'_i)$, where $p'_i$ is the signal energy in the direction, $\theta_i$. Peak detection is performed on the energy vector, ${\mathbf{p}}_{{d}}$, where the peak energy occurs at the source signal, enabling the DOA estimation.

In this experimental scenario, $4$ source signals arrived from the directions: $\xi_1 = 30^{\circ}$, $\xi_2 = 60^{\circ}$, $\xi_3 = 120^{\circ}$, and $\xi_4 = 160^{\circ}$, with corresponding signal energies of $15$, $5$, $5$, and $15$, respectively. Let us consider the case in which the linear array consists of ${v} = 15$ sensors and $n$ independent array observations are simulated by the model (5.1) to obtain the sample covariance matrix. By applying Algorithms 1 and 2 to estimate ${\bf{\Sigma}}$ in (5.3), we obtain an estimate of ${\mathbf{p}}_{{d}}$, denoted by ${\mathbf{p}}_{{d}}^*$. By performing peak detection on ${\mathbf{p}}_{{d}}^*$, we obtained the DOA estimates: $\xi^*_1, \cdots, \xi_4^*$, assuming the number of true DOAs is known in advance.

Table 1 presents a performance comparison of Algorithms 1 and 2 in terms of computational complexity and estimation accuracy, where

where ${\mathbf{p}}_{(i)}^*$ denotes the estimate of ${\mathbf{p}}$ in (5.3) in the $i$-th Monte Carlo simulation. These results demonstrate that the MM algorithm (Algorithm 2) exhibits better computational efficiency than the PDIP (primal-dual interior-point) algorithm (Algorithm 1), and that their computational accuracies are comparable. In Figure 1, the DOA estimation performances of the MM method (Algorithm 2) and the classical DOA estimation methods: Likelihood-based estimation of sparse parameters (LIKES) in ^[12] and sparse iterative covariance-based estimation (SPICE) in ^[11], are compared, where

and $(\xi_j^*)_{(i)}$ is the estimate of $\xi_j$ in the $i$-th Monte Carlo simulation. It is observed that the MM method offers higher DOA estimation accuracy. In Figure 2, the normalized spectrum is defined to be ${\mathbf{p}}^* / \max({\mathbf{p}}^*)$, where $\max({\mathbf{p}}^*)$ denotes the maximal element in the vector, ${\mathbf{p}}^*$. The black dashed lines represent the true DOAs and the normalized spectrum curves represent the results of a randomly selected realization. The location of the peak spectrum lies very close to the true DOA, and a higher SNR is observed corresponding to a higher angular resolution.

5.2.   Toeplitz correlation matrix estimation

In a $q$-order moving average time series (MA($\hat{q}$)),

where $\phi_i = \phi^i$, $i = 1, \cdots, \hat{q}$, and $\varepsilon_t, \cdots, \varepsilon_{t-\hat{q}}$ are independently and identically distributed with zero mean and variance 1. A simple calculation yields that the ${v} \times {v}$ autocorrelation matrix ${\mathbf{C}}$ exhibits a Toeplitz structure with decreasing autocorrelation coefficients as the time lag increases. Denote the first column of ${\mathbf{C}}$ as ${\mathbf{c}}$, and use $\hat{q} = 5$. This experiment assessed the performance of the MM-T algorithm in estimating this autocorrelation matrix, focusing on its estimation accuracy and computation time, and compared them with these of other methods:

● PDIP-T: the primal-dual interior-point algorithm solving (4.8);

● Samp: the estimator of Toeplitz covariance matrices mentioned in ^[16] that minimizes the F-norm loss between the sample covariance matrices;

● HY: the estimation of Toeplitz covariances via the entropy loss function using the alternating direction method of multipliers algorithm given in ^[21].

Note that the sample autocorrelation matrix is used as the given correlation matrix ${\mathbf{R}}$ in the MM-T and PDIP-T algorithms, and the sample autoconvariance matrix is used as the given covariance matrix in the HY method, and that the sample size is $n = 1000$.

Table 2 lists the RMSE$_{\text{toep}}$ corresponding to different $l$ values under different ${v}$ and $\phi$. The RMSE$_{\text{toep}}$ is defined as

where ${\mathbf{c}}_{(i)}^*$ denotes the estimate of ${\mathbf{c}}$ achieved during the $i$-th Monte Carlo repetition. For different $l$, $\text{RMSE}_{\text{toep}}$ values remain stable. Therefore, we assumed $l = 10$ to complete this experiment.

Figure 3(a) depicts the computation times of the MM-T and PDIP-T algorithms. The results indicate that as the order of the autocorrelation matrix increases (i.e., the considered time lag ${v}-1$ increases), the rate of increase in the computational time of MM-T becomes slower than that of PDIP-T, and when ${v} \geq 50$, MM-T requires less time.

Figure 3(b) compares the estimation accuracies of the different methods. The MM-T and PDIP-T algorithms exhibit lower $\text{RMSE}_{\text{toep}}$ values, which is reasonable as they utilize a decreasing correlation structure, unlike the Samp and HY methods. The estimation accuracies of the MM-T and PDIP-T algorithms are essentially equal; however, when ${v}$ is moderate to large, the computation time required by MM-T is shorter. Therefore, the MM-T algorithm is recommended over PDIP-T.

5.3.   Real data analysis

We collected the closing stock prices of 511 companies in the technology field on Fridays during the period between 2023-9-15 and 2023-11-17 (8 weeks). We compute the volatility series as follows:

where $\text{Close}_{\text{current Friday}}$ and $\text{Close}_{\text{previous Friday}}$ denote the closing stock prices on current and previous Fridays, respectively, and $\text{Volatility}_{\text{current Friday}}$ denotes the volatility of closing stock prices on Fridays.

Subsequently, to investigate the temporal dependencies in the volatility series, we assume that the autocorrelation matrix of the volatility series exhibits a Toeplitz structure: We apply the MM-T and HY algorithms to estimate the autocorrelation matrix. Table 3 presents the first column of the autocorrelation matrix estimation; i.e., the autocorrelation coefficients corresponding to different time lags.

Table 3 shows that the results of the MM-T and HY methods are similar. Moreover, the autocorrelation coefficients corresponding to different lags are nearly identical, but lie around a small value of $0.22$. The small autocorrelation coefficients indicate no obvious dependency between volatilities on different Fridays, which is consistent with the unpredictable nature of stock price fluctuations.

6.   Conclusions

We focus on estimating of the W-Rank1-W structural covariance matrix, which appears commonly in signal processing, and the auto-correlation matrix characterized by decreasing correlation coefficients, which appears commonly in time-series analysis. To this end, novel algorithms are proposed using the MM algorithm framework. Numerical experiments reveal that the proposed MM algorithms exhibit superior computational efficiencies. The proposed strategy for dealing with the surrogate function under the MM framework can be applied to other similar convex optimization problems.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported in part by the National Nature Science Foundation of China under Grants 62203451 and the Sichuan Province Science and Technology Support Program under Grants 2022JDRC0068 and 2021JDRC0080. The Fundamental Research Funds for the Central Universities 24CAFUC03055.

Conflict of interest

The authors declare that they have no competing interests.