High-Dimensional Covariance Estimation via Constrained L_q-Type Regularization

Abstract

High-dimensional covariance matrix estimation is one of the fundamental and important problems in multivariate analysis and has a wide range of applications in many fields. In practice, it is common that a covariance matrix is composed of a low-rank matrix and a sparse matrix. In this paper we estimate the covariance matrix by solving a constrained $L_q$-type regularized optimization problem. We establish the first-order optimality conditions for this problem by using proximal mapping and the subspace method. The proposed stationary point degenerates to the first-order stationary points of the unconstrained $L_q$ regularized sparse or low-rank optimization problems. A smoothing alternating updating method is proposed to find an estimator for the covariance matrix. We establish the convergence of the proposed calculation method. The numerical simulation results show the effectiveness of the proposed approach for high-dimensional covariance estimation.

1. Introduction

High-dimensional covariance matrix estimation is an important topic in statistical inference and learning. The objective is to estimate the covariance matrix of p variables based on sample data of size n. Much statistical analysis of high-dimensional data involves covariance estimation. In a high-dimensional setting (usually $p\gg n$), researchers usually assume that the covariance matrix is sparse and have proposed various regularization techniques [1,2,3,4,5,6,7,8,9] to consistently estimate covariance matrix $\Sigma$. However, this assumption is not appropriate in many applications such as financial analysis, depending on the equity market risks, gene expressions stimulated by cytokines, etc. A more reasonable assumption based on the factor model [10,11,12,13,14] is that the high-dimensional covariance matrix is the sum of a low-rank matrix and a sparse matrix. Specifically, a factor model can be formulated as

$$\begin{array}{c}\hfill {\mathit{Y}}_i=\mathit{B}{\mathit{f}}_i+{\mathit{u}}_i,i=1,\cdots ,n,\end{array}$$

(1)

where ${\mathit{Y}}_i={(Y_{1i},\cdots ,Y_{pi})}^T$ is the observed variable, $\mathit{B}={({\mathit{b}}_1,\cdots ,{\mathit{b}}_p)}^T$ is the matrix of factor loadings, ${\mathit{f}}_i$ is a $K\times 1$ vector of common factors, and ${\mathit{u}}_i={(u_{1i},\cdots ,u_{pi})}^T$ is the idiosyncratic error component, uncorrelated with ${\mathit{f}}_i$. Under this model, the covariance matrix of ${\mathit{Y}}_i$ is given by

$$\begin{array}{ccc}\hfill \Sigma & =& \mathit{B}\mathrm{cov}\left({\mathit{f}}_i\right){\mathit{B}}^T+{\Sigma }_u\hfill \\ & :=& \mathit{L}+\mathit{S},\hfill \end{array}$$

(2)

where $\mathit{L}$ is the covariance matrix of $\mathit{B}{\mathit{f}}_i$ and $\mathit{S}$ is the covariance matrix of ${\mathit{u}}_i$. Since K is usually less than p, $\mathit{L}$ is a low-rank matrix. More generally, there is a large class of problems where one is interested in estimating the covariance matrix of ${\mathit{X}}_i$ from the noisy observations ${\mathit{Y}}_i={\mathit{X}}_i+{\mathit{u}}_i$, where the noise ${\mathit{u}}_i$ is uncorrelated with the signal ${\mathit{X}}_i$. In addition, if ${\mathit{X}}_i$ is generated from a linear system so that some of its elements are linear functions of the others, then the covariance matrix of ${\mathit{X}}_i$ has low-rank. Based on the structural decomposition (2), this paper proposes a novel approach to estimate the covariance matrix by solving a constrained low-rank and sparse regularized optimization problem.

Under structural decomposition (2), researchers have proposed some approaches to estimate the covariance matrix. When the common factors in Model (1) are known and observable and $\mathit{S}$ is diagonal, Reference [10] proposed a least-squares-based covariance estimation (LS-CE). The core idea of LS-CE is first to estimate $\mathit{B}$ by the least-squares method and, then, obtain the final estimator assuming that $\mathit{S}$ is diagonal. Compared to the sample covariance matrix estimator, the LS-CE estimator is always invertible, even in a high-dimensional setting. Moreover, this estimator is shown to have substantial gains when estimating the inverse of the covariance matrix, however, it does not have much advantage when estimating the covariance matrix. Instead of assuming diagonal $\mathit{S}$, Reference [11] assumed that $\mathit{S}$ is sparse and applied the adaptive thresholding technique of [4] to obtain an estimator for $\mathit{S}$. Then, a least-squares-based covariance estimator with conditional sparsity (LS-CS-CE) under Model (1) was studied in [11]. More recently, a robust approach based on Huber loss minimization was proposed in [13]. This approach estimates the initial joint covariance matrix of the observed data and the factors and, then, uses it to recover the covariance matrix of the observed data. Moreover, the proposed estimator in [13] only requires the condition that the fourth moment of the data exists. This method is applicable to a wider range of data, including sub-Gaussian and elliptical distributions. When the common factors are unknown and unobservable, Reference [12] proposed a non-parametric estimator via the thresholding principal orthogonal complements, which is called the principal orthogonal complement thresholding (POET) method. The POET uses the K principal components to estimate $\mathit{L}$ and applies the adaptive thresholding technique [4] on the rest of the $(p-K)$ components to estimate $\mathit{S}$. However, as stated in [9], using the simple thresholding approaches in [1,2,3,4,5,9] to estimate $\mathit{S}$ may lead to the covariance estimators in [11,12,13], which are necessarily positive semi-definite.

Moreover, as suggested in [14], an intuitive optimization-type approach is to solve the problem:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{\mathit{L},\mathit{S}\in {\mathbb{R}}^{p\times p}}{min}& \hfill \frac{1}{2}\parallel \mathit{L}+\mathit{S}-\tilde{\Sigma }{\parallel }_F^2+{\lambda }_L\mathrm{rank}\left(\mathit{L}\right)+{\lambda }_S\sum _{(u,v)\in D^c}I\left(\right|S_{uv}\left|\right),\hfill \end{array}\end{array}$$

(3)

where $\tilde{\Sigma }$ is the sample covariance matrix, $\mathrm{rank}\left(\mathit{L}\right)$ denotes the rank of matrix $\mathit{L}$, $I\left(\right|S_{uv}\left|\right)=1$ if $|S_{uv}|\ne 0$ and $I\left(\right|S_{uv}\left|\right)=0$ if $|S_{uv}|=0$, $D^c:=\{(u,v)\in \mathbb{N}\times \mathbb{N}:1\le u\le p,1\le v\le p,u\ne v\}$ is the index set of non-diagonal elements of $\mathit{S}$, and ${\lambda }_L$ and ${\lambda }_S$ are tuning parameters. To facilitate the calculation, Reference [14] proposed the low-rank and sparse covariance estimator (LOREC) via solving the following convex relaxation optimization problem:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{\mathit{L},\mathit{S}\in {\mathbb{R}}^{p\times p}}{min}& \hfill \frac{1}{2}\parallel \mathit{L}+\mathit{S}-\tilde{\Sigma }{\parallel }_F^2+{\lambda }_R{\parallel \mathit{L}\parallel }_{\ast }+{\lambda }_S\sum _{(u,v)\in D^c}\left|S_{uv}\right|,\hfill \end{array}\end{array}$$

(4)

where ${\parallel \mathit{L}\parallel }_{\ast }:={\sum }_{i=1}^p{\sigma }_i\left(\mathit{L}\right)$ denotes the nuclear norm of $\mathit{L}$ and ${\sigma }_i\left(\mathit{L}\right)$ is the ith largest singular value of $\mathit{L}$. However, the estimators for $\mathit{L}$ and $\mathit{S}$ proposed in [14] are not guaranteed to be positive semi-definite. Clearly, the optimization problems (3) and (4) can be considered as special cases of the robust principal component analysis (RPCA) in [15] and its convex relaxation, respectively.

In this paper, we propose an optimization model similar to some models of RPCA, which is formulated as

$$\begin{array}{c}\hfill \underset{\mathit{L},\mathit{S}\in {\mathbb{R}}^{p\times z}}{min}\mathrm{rank}\left(\mathit{L}\right)+{\lambda }_S{\parallel \mathit{S}\parallel }_0,\mathrm{s}.\mathrm{t}.\mathit{L}+\mathit{S}=\tilde{\Sigma },\end{array}$$

where ${\parallel \mathit{S}\parallel }_0$ denotes the number of the non-zero elements of $\mathit{S}$. Reference [15] showed that the principal component pursuit estimate solving

$$\begin{array}{c}\hfill \underset{\mathit{L},\mathit{S}\in {\mathbb{R}}^{p\times z}}{min}{\parallel \mathit{L}\parallel }_{\ast }+{\lambda }_S\sum _{u=1}^p\sum _{v=1}^p\left|S_{uv}\right|,\mathrm{s}.\mathrm{t}.\mathit{L}+\mathit{S}=\tilde{\Sigma },\end{array}$$

(5)

exactly recovers the low-rank $\mathit{L}$ and the sparse $\mathit{S}$. Furthermore, Reference [16] suggested that $\mathit{L}$ and $\mathit{S}$ can be estimated by solving the following unconstrained optimization problem:

$$\begin{array}{c}\hfill \underset{\mathit{L},\mathit{S}\in {\mathbb{R}}^{p\times z}}{min}{\parallel \mathit{L}\parallel }_{\ast }+{\lambda }_S\sum _{u=1}^p\sum _{v=1}^p|S_{uv}|+\frac{{\lambda }_L}{2}{\parallel \mathit{L}+\mathit{S}-\tilde{\Sigma }\parallel }_F^2.\end{array}$$

(6)

More references about the RPCA and related works based on convex and non-convex regularization can be found in [17,18,19,20,21].

The main contributions of this paper can be summarized as follows:

• We construct a constrained adaptive $L_q$-type regularized optimization problem for covariance estimation, which is a non-smooth and non-convex problem. This optimization problem is based on the structural decomposition (2) and the idea of RPCA in [16], but is not limited to the factor model (1).
• We study the first-order optimality condition of this problem and define a class of hybrid first-order stationary points. We establish the relationship between the first-order stationary point and the global minimizer of this problem.
• We propose a smoothing alternating updating method to solve the constructed non-smooth, non-convex optimization problem and established its convergence. The simulation results show that the proposed smoothing alternating updating method is efficacious for the constrained adaptive $L_q$-type regularized optimization problem.

The rest of this paper is organized as follows. In Section 2, we give some notations and preliminaries. In Section 3, we construct a constrained adaptive $L_q$-type regularized optimization problem for covariance estimation. In Section 4, we study the first-order optimality condition of this problem and define a class of hybrid first-order stationary points. In Section 5, we propose a smoothing alternating updating method to solve the non-smooth, non-convex optimization problem and establish its convergence. The simulation results are given in Section 6. The conclusions and discussion are given in Section 7, while the mathematical proofs are given in Section 8.

2. Notations and Preliminaries

Throughout this paper, we use the following notations. The sets of $p\times p$ positive semi-definite symmetric and positive definite symmetric matrices are denoted by ${\mathbb{S}}_{+}^p$ and ${\mathbb{S}}_{++}^p$, respectively. It is known that ${\mathbb{S}}_{+}^p$ and ${\mathbb{S}}_{++}^p$ are cones.

Given matrices $\mathit{X}:={\left(X_{uv}\right)}_{u,v=1}^p\in {\mathbb{R}}^{p\times p}$ and $\mathit{Y}:={\left(Y_{uv}\right)}_{u,v=1}^p\in {\mathbb{R}}^{p\times p}$, where $X_{uv}$ is the $(u,v)$th element of $\mathit{X}$, $\langle \mathit{X},\mathit{Y}\rangle :={\sum }_{u=1}^p{\sum }_{v=1}^p{X}_{uv}{Y}_{uv}$, and $\mathit{X}\circ \mathit{Y}$ denotes the Hadamard product ${(\mathit{X}\circ \mathit{Y})}_{uv}=:X_{uv}{Y}_{uv}$ for all $u,v=1,\cdots ,p$. Define the elementwise maximum norm (or max norm) ${\left|\mathit{X}\right|}_{\max }:={max}_{u,v}\left|X_{uv}\right|$, the Frobenius norm ${\parallel \mathit{X}\parallel }_F:={\left({\sum }_{u=1}^p{\sum }_{v=1}^p{X}_{uv}^2\right)}^{1/2}$, and the spectral norm ${\parallel \mathit{X}\parallel }_2:={sup}_{{\parallel \mathit{x}\parallel }_2\le 1}{\parallel \mathit{X}\mathit{x}\parallel }_2$, where ${\parallel \mathit{x}\parallel }_2:=\sqrt{{\sum }_{i=1}^p{x}_i^2}$ is 2-norm of vector $\mathit{x}$. Moreover, given a matrix $\mathit{X}\in {\mathbb{R}}^{p\times p}$ and a index set $\mathsf{\Gamma }\subseteq \left\{\right(u,v):1\le u\le p,1\le v\le p\}$, ${\mathit{X}}_{\mathsf{\Gamma }}$ denotes the array of the same dimension such that ${\left[{\mathit{X}}_{\mathsf{\Gamma }}\right]}_{uv}=X_{uv}$ for $(u,v)\in \mathsf{\Gamma }$ and ${\left[{\mathit{X}}_{\mathsf{\Gamma }}\right]}_{uv}$ are undefined for $(u,v)\in {\mathsf{\Gamma }}^c$. For convenience, we will write the $(u,v)$th element of ${\mathit{X}}_{\mathsf{\Gamma }}$ as $X_{uv}$ regardless if it is defined or not.

Given a vector $\mathit{x}\in {\mathbb{R}}^p$, $\mathrm{Diag}\left(\mathit{x}\right)$ denotes a square matrix with vector $\mathit{x}$ on its diagonal and zeros elsewhere. Given a matrix $\mathit{X}\in {\mathbb{S}}_{+}^p$ and letting $D_i\left(\mathit{X}\right)$ denote the ith largest eigenvalue of $\mathit{X}$, $i=1,\cdots ,p$, write $D\left(\mathit{X}\right)={(D_1\left(\mathit{X}\right),\cdots ,D_p\left(\mathit{X}\right))}^T\in {\mathbb{R}}^p$. Define sign operator:

$$\mathrm{sign}\left(t\right)=\left\{\begin{array}{cc}\hfill 1,\hfill & \hfill \mathrm{if}t>0\hfill \\ \hfill 0,\hfill & \hfill \mathrm{if}t=0\hfill \\ \hfill -1,\hfill & \hfill \mathrm{if}t<0,\hfill \end{array}\right.$$

and set $D:=\left\{\right(u,v)\in \mathbb{N}\times \mathbb{N}:1\le u\le p,1\le v\le p,u=v\}$.

Let $f:{\mathbb{R}}^{p\times p}\to \overline{\mathbb{R}}$ be a convex function, and let $\overline{\mathit{X}}\in \mathrm{dom}f$. The subdifferential of f at $\overline{\mathit{X}}$ is defined by

$$\partial f\left(\overline{\mathit{X}}\right):=\left\{\mathit{\xi }\in {\mathbb{R}}^{p\times p}:\underset{t\downarrow 0}{lim}\frac{f(\overline{\mathit{X}}+t\mathit{Z})-f\left(\overline{\mathit{X}}\right)}{t}\ge \langle \mathit{\xi },\mathit{Z}\rangle \mathrm{for}\mathrm{all}\mathit{Z}\in {\mathbb{R}}^{p\times p}\right\},$$

where $\mathrm{dom}f:=\{\mathit{X}\in {\mathbb{R}}^{p\times p}:f\left(\mathit{X}\right)<\infty \}$ denotes the domain of f, $\overline{\mathbb{R}}:=(-\infty ,\infty ]$. In other words, $\partial f\left(\overline{\mathit{X}}\right)$ is the closed convex hull of all points of the form ${lim}_{i\to \infty }\nabla f\left({\mathit{X}}_i\right)$, where $\left\{{\mathit{X}}_i\right\}$ is any sequence that converges to $\overline{\mathit{X}}$ while avoiding the set for which f fails to be differentiable. Moreover, the singular subdifferential of f at $\overline{\mathit{X}}$ is defined by

$${\partial }^{\infty }f\left(\overline{\mathit{X}}\right):=\left\{\mathit{\xi }\in {\mathbb{R}}^{p\times p}:(\mathit{\xi },0)\in N\left((\overline{\mathit{X}},f\left(\overline{\mathit{X}}\right));\mathrm{epi}f\right)\right\},$$

where $\mathrm{epi}f:=\{(\mathit{X},t)\in {\mathbb{R}}^{p\times p}\times \mathbb{R}:\mathit{X}\in \mathrm{dom}f,t\ge f\left(\mathit{X}\right)\}$ denotes the epigraph of f.

Finally, we give the definition of the proximal mapping in this paper. Given a closed set $\mathsf{\Omega }$ and a function $f:{\mathbb{R}}^{p\times p}\to \overline{\mathbb{R}}$, the proximal mapping on $\mathsf{\Omega }$ of f is the operator given by

$${\mathrm{Prox}}_{f,\mathsf{\Omega }}\left(\mathit{X}\right)=\underset{\mathit{U}\in \mathsf{\Omega }}{\mathrm{argmin}}\left\{f\left(\mathit{U}\right)+\frac{1}{2}{\parallel \mathit{U}-\mathit{X}\parallel }_F^2\right\}\mathrm{for}\mathrm{any}\mathit{X}\in {\mathbb{R}}^{p\times p}.$$

3. Covariance Estimation with Adaptive ${\mathit{L}}_{\mathit{q}}$-Type Regularization

Considering that sparse covariance matrix estimation is intrinsically a heteroscedastic problem ([4]) and $\mathit{L}$ and $\mathit{S}$ are all positive semi-definite matrices, we estimate the covariance matrix by solving the following adaptive $L_q$-type regularized optimization problem:

$$\begin{array}{cc}\hfill \underset{\mathit{L},\mathit{S}\in {\mathbb{S}}_{+}^p}{min}& \hfill F(\mathit{L},\mathit{S}):=\frac{1}{2}\parallel \mathit{L}+\mathit{S}-\tilde{\Sigma }{\parallel }_F^2+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q+\sum _{(u,v)\in D^c}{\lambda }_{uv}^S{\left|S_{uv}\right|}^q,\hfill \end{array}$$

(7)

where $q\in (0,1)$, ${\parallel \mathit{L}\parallel }_q^q:={\sum }_{i=1}^p{\sigma }_i{\left(\mathit{L}\right)}^q$ is the Schatten-q quasi-norm of matrix $\mathit{L}$ and ${\lambda }_L$ and ${\lambda }_{uv}^S$ are tuning parameters. Note that although we use the same q for $\mathit{L}$ and $\mathit{S}$ here for simplicity, it can be different in general. Our approach is called $L_q$ regularized covariance ($L_q$-CSC) estimation with conditional sparsity. The final estimator $\widehat{\Sigma }:=\widehat{\mathit{L}}+\widehat{\mathit{S}}$ is called the $L_q$ regularized covariance estimator with conditional sparsity ($L_q$-CSCE).

To deal with heteroscedasticity in sparse covariance matrix estimation, we adopt the adaptive thresholding approach of Reference [4] and [12] to calculate the entry-dependent threshold ${\lambda }_{uv}^S$ in (7) for $S_{uv}$ as

$$\begin{array}{c}\hfill {\lambda }_{uv}^S:=\lambda {\left(\frac{{\tilde{\sigma }}_{uu}{\tilde{\sigma }}_{vv}logp}{n}\right)}^{1/2},\end{array}$$

(8)

where $\lambda >0$ is a constant. Note that these entry-dependent tuning parameters for $\mathit{S}$ in Problem (7) depend on the variance of the corresponding variables. The parameter ${\tilde{\sigma }}_{uu}$ can be taken as the u-th diagonal elements of the estimator of the POET for $\mathit{S}$. It is worth mentioning that, when the data matrix is normalized, the covariance matrix is equivalent to its correlation matrix, in which case, all tuning parameters are $\lambda \sqrt{logp/n}$, and the adaptive $L_q$ regularization function acting on $\mathit{S}$ becomes the ordinary $L_q$ regularization function.

We now discuss the selection of $\tilde{\Sigma }$. The common one is the ordinary sample covariance matrix

$$\frac{1}{n-1}\sum _{i=1}^n({\mathit{Y}}_i-\overline{\mathit{Y}}){({\mathit{Y}}_i-\overline{\mathit{Y}})}^T,$$

where $\overline{\mathit{Y}}=n^{-1}{\sum }_{i=1}^n{\mathit{Y}}_i$, ${\mathit{Y}}_i$ is the i-th sample vector, $i=1,\cdots ,n$. The ordinary sample covariance is suitable for the data with exponential-type tails or polynomial-type tails; see [1,4,8]. However, in many practical problems data tend to be heavy-tailed. References [9,13,22] defined a Huber minimization sample covariance matrix, which is suitable for the data with finite fourth moment. Moreover, Reference [9] also defined a rank-based covariance matrix and a median of means covariance matrix for data with other distribution conditions.

Finally, we discuss some basic properties of Problem (7).

Proposition 1. 

The following statements of $F(\mathit{L},\mathit{S}):{\mathbb{R}}^{p\times p}\times {\mathbb{R}}^{p\times p}\to {\mathbb{R}}_{+}$ in Problem (7) hold:

(i)

$F(\mathit{L},\mathit{S})$ is proper, closed, and coercive, i.e., ${lim}_{{\parallel (\mathit{R},\mathit{S})\parallel }_F\to \infty }F(\mathit{L},\mathit{S})=\infty$.

(ii)

$F(\mathit{L},\mathit{S})$ attains its minimal value over ${\mathbb{S}}_{+}^p\times {\mathbb{S}}_{+}^p$.

(iii)

Let $f(\mathit{L},\mathit{S})=\frac{1}{2}{\parallel \mathit{L}+\mathit{S}-\tilde{\Sigma }\parallel }_F^2$; the gradient $\nabla f(\mathit{L},\mathit{S})$ is Lipschitz continuous with a constant $\ell =2$, that is, for any $({\mathit{L}}_1,{\mathit{S}}_1)$, $({\mathit{L}}_2,{\mathit{S}}_2)\in {\mathbb{R}}^{p\times p}\times {\mathbb{R}}^{p\times p}$, it holds that

$$\parallel \nabla f({\mathit{L}}_1,{\mathit{S}}_1)-\nabla f({\mathit{L}}_2,{\mathit{S}}_2){\parallel }_F\le 2{\parallel ({\mathit{L}}_1,{\mathit{S}}_1)-({\mathit{L}}_2,{\mathit{S}}_2)\parallel }_F.$$

(iv)

Given a matrix ${\mathit{L}}^{\ast }\in {\mathbb{R}}^{p\times p}$, the partial gradient function ${\nabla }_{S}f({\mathit{L}}^{\ast },\mathit{S})$ is Lipschitz continuous with a constant $\ell =1$, that is, for any ${\mathit{S}}_1$, ${\mathit{S}}_2\in {\mathbb{R}}^{p\times p}$, it holds that

$$\parallel {\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}_1)-{\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}_2){\parallel }_F\le {\parallel {\mathit{S}}_1-{\mathit{S}}_2\parallel }_F.$$

Moreover, given a matrix ${\mathit{S}}^{\ast }\in {\mathbb{R}}^{p\times p}$, the partial gradient function ${\nabla }_{L}f(\mathit{L},{\mathit{S}}^{\ast })$ is Lipschitz continuous with a constant $\ell =1$.

4. Necessary Optimality Conditions

In this section, we define a class of first-order stationary points for Problem (7), which is a hybridization of the first-order stationary points in [23,24,25] and the fixed point in [25]. Furthermore, we prove that any global minimizer of Problem (7) is a first-order stationary point of Problem (7). We also define an approximate variant of this stationary point. It is worth mentioning that the constraints in Problem (7) are treated as ordinary closed convex set constraints in this paper.

For the simple unconstrained $L_q$ regularized sparse optimization problem [23,24]:

$$\underset{\mathit{x}\in {\mathbb{R}}^p}{min}f\left(\mathit{x}\right)+\lambda {\parallel \mathit{x}\parallel }_q^q,$$

where ${\parallel \mathit{x}\parallel }_q^q:={\sum }_{i=1}^p{\left|x_i\right|}^q$, its the first-order stationary point is defined as the point satisfying

$$\mathrm{Diag}\left({\mathit{x}}^{\ast }\right)\nabla f\left({\mathit{x}}^{\ast }\right)+\lambda q{\left|{\mathit{x}}^{\ast }\right|}^q=\mathit{0}.$$

For the simple unconstrained $L_q$ regularized rank optimization problem [25]:

$$\underset{\mathit{X}\in {\mathbb{R}}^{m\times n}}{min}f\left(\mathit{X}\right)+\lambda {\parallel \mathit{X}\parallel }_q^q,$$

its first-order stationary point (or fixed point) is defined as the point satisfying

$${\mathit{X}}^{\ast }\in {\mathrm{Prox}}_{{\lambda \parallel \mathit{X}\parallel }_q^q/\ell ,{\mathbb{R}}^{p\times p}}\left({\mathit{X}}^{\ast }-\frac{1}{\ell }\nabla f\left({\mathit{X}}^{\ast }\right)\right),$$

where ℓ is a constant greater than the Lipschitz constant of $\nabla f\left(\mathit{X}\right)$. Inspired by these works, we can define a class of hybrid first-order stationary points of Problem (7).

Definition 1. 

$({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })\in {\mathbb{S}}_{+}^p\times {\mathbb{S}}_{+}^p$ is said to be the first-order stationary point of Problem (7) if $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ satisfies

$$\begin{array}{ccc}& & \hfill {\mathit{L}}^{\ast }\in {\mathrm{Proxi}}_{({\lambda }_L/\ell ){\parallel \mathit{L}\parallel }_q^q,{\mathbb{S}}_{+}^p}\left({\mathit{L}}^{\ast }-\frac{1}{\ell }({\mathit{L}}^{\ast }+{\mathit{S}}^{\ast }-\tilde{\Sigma })\right),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& & \hfill \mathit{0}\in \left\{{\mathit{S}}^{\ast }\circ ({\mathit{L}}^{\ast }+{\mathit{S}}^{\ast }-\tilde{\Sigma })+{\mathit{Q}}^{\ast }+{\mathit{S}}^{\ast }\circ \mathit{Y}:{\mathit{Y}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\in \partial h_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\left({\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\right)\right\},\hfill \end{array}$$

(10)

where ${\mathit{Q}}^{\ast }:={\left(Q_{uv}^{\ast }\right)}_{u,v=1}^p$ and

$$\begin{array}{c}\hfill Q_{uv}^{\ast }=\left\{\begin{array}{cc}\hfill q{\lambda }_{uv}^S{\left|S_{uv}^{\ast }\right|}^q,\hfill & \hfill \mathrm{if}(u,v)\in D^c,\hfill \\ \hfill 0,\hfill & \hfill \mathrm{if}(u,v)\in D,\hfill \end{array}\right.\end{array}$$

${\mathsf{\Gamma }}^{\ast }=\{(u,v):S_{uv}^{\ast }\ne 0\}$, $h_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\left({\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\right):={\delta }_{{\mathbb{S}}_{+}^p}\left([{\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}:\mathit{0}]\right)$, $\mathit{Y}=[{\mathit{Y}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}:{\mathit{Y}}_{{({\mathsf{\Gamma }}^{\ast }\bigcup D)}^c}]$, ${\mathit{Y}}_{{({\mathsf{\Gamma }}^{\ast }\bigcup D)}^c}$ can be any finite array.

Next, we investigate the relationship between the global minimizer of Problem (7) and the first-order stationary point of Problem (7) without any qualification conditions.

Theorem 1. 

Suppose that $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })\in {\mathbb{S}}_{+}^p\times {\mathbb{S}}_{+}^p$ is a global minimizer of Problem (7), and let ${\mathsf{\Gamma }}^{\ast }=\{(u,v)\in \mathbb{N}\times \mathbb{N}:S_{uv}^{\ast }\ne 0\}$. Then, $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ is a first-order stationary point with $\ell >1$ of Problem (7).

Again, under some qualification conditions, the form of the first-order stationary point of Problem (7) can be further analyzed. We now discuss the specific form of the subdifferential of $h_{{\mathsf{\Gamma }}^{\ast }\bigcup D}$ at $S_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }$. Note that

$$h_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\left({\mathit{X}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\right)={\delta }_{{\mathbb{S}}_{+}^p}\left([{\mathit{X}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}:\mathit{0}]\right)={\delta }_{{\mathbb{S}}_{+}^p}\left({\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\left({\mathit{X}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\right)\right),$$

where ${\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\left({\mathit{X}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\right):=[{\mathit{X}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}:\mathit{0}]$ is the mapping acting on array ${\mathit{X}}_{{\mathsf{\Gamma }}^{\ast }}$. The adjoint mapping ${\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }:{\mathbb{R}}^{p\times p}\to {\mathbb{R}}^{p\times p}$ is defined by

$$\langle {\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\left({\mathit{X}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\right),\mathit{Y}\rangle =\langle {\mathit{X}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D},{\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\left(\mathit{Y}\right)\rangle \mathrm{for}\mathrm{all}{\mathit{X}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\mathit{Y}\in {\mathbb{R}}^{p\times p}.$$

Specifically, ${\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\left(\mathit{Y}\right)={\mathit{Y}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}$. Clearly, ${\delta }_{{\mathbb{S}}_{+}^p}$ is convex, and ${\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}$ and ${\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }$ are linear mappings. It follows from Theorem 23.9 in [26] and Theorem 2.51 in [27] that, if ${\mathit{S}}^{\ast }\in {\mathbb{S}}_{+}^p$, then

$$\begin{array}{c}\hfill \{{\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\left(\mathit{Y}\right):\mathit{Y}\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)\}\subseteq \partial h_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\left({\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\right).\end{array}$$

(11)

Moreover, suppose that the qualification condition $\mathrm{Ker}{\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\bigcap {\partial }^{\infty }{\delta }_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)=\left\{\mathit{0}\right\}$ holds, where $\mathrm{Ker}{\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }$ is the kernel of linear mapping ${\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }$, then

$$\begin{array}{c}\hfill \{{\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\left(\mathit{Y}\right):\mathit{Y}\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)\}=\partial h_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\left({\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\right).\end{array}$$

(12)

Based on these discussions, the following statement holds: If $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ satisfies

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \hfill {\mathit{L}}^{\ast }\in {\mathrm{Proxi}}_{{\lambda }_L{\parallel \mathit{L}\parallel }_q^q,{\mathbb{S}}_{+}^p}(\tilde{\Sigma }-{\mathit{S}}^{\ast }),\hfill \\ \hfill & \hfill \mathit{0}\in \left\{{\mathit{S}}^{\ast }\circ ({\mathit{L}}^{\ast }+{\mathit{S}}^{\ast }-\tilde{\Sigma })+{\mathit{Q}}^{\ast }+{\mathit{S}}^{\ast }\circ \mathit{Y}:\mathit{Y}\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)\right\},\hfill \end{array}\end{array}$$

where ${\mathit{Q}}^{\ast }$ is defined in Definition 1 and the qualification condition $\mathrm{Ker}{\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\bigcap {\partial }^{\infty }{\delta }_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)=\left\{\mathit{0}\right\}$ holds, then $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ is a first-order stationary point of Problem (7).

Obviously, the qualification condition $\mathrm{Ker}{\mathcal{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\bigcap {\partial }^{\infty }{\delta }_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)=\left\{\mathit{0}\right\}$ holds automatically if ${\delta }_{{\mathbb{S}}_{+}^p}$ is locally Lipschitz continuous at ${\mathit{S}}^{\ast }$. In other words, if $S^{\ast }\in \mathrm{int}\left({\mathbb{S}}_{+}^p\right)={\mathbb{S}}_{++}^p$, this qualification condition holds. Indeed, the true covariance matrix ${\mathit{S}}^o$ is often assumed to be positive definite [7,8]. When ${\mathit{S}}^{\ast }\in {\mathbb{S}}_{++}^p$, we can establish the lower bound theory of the first-order stationary point of Problem (7).

Theorem 2. 

Let $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ be a first-order stationary point of Problem (7) satisfying $F({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })\le F({\mathit{L}}^0,{\mathit{S}}^0)+ϵ$ for some ${\mathit{L}}^0,{\mathit{S}}^0$ and $ϵ\ge 0$, ${\mathsf{\Gamma }}^{\ast }:=\{(u,v):S_{uv}^{\ast }\ne 0\}$. If ${\mathit{S}}^{\ast }\in {\mathbb{S}}_{++}^p$, then

$$|S_{uv}^{\ast }|\ge {\left(\frac{qmin\left\{{\lambda }_{uv}^S\right\}}{\sqrt{2[F({\mathit{L}}^0,{\mathit{S}}^0)+ϵ]}}\right)}^{\frac{1}{1-q}}\forall (u,v)\in D^c\bigcap {\mathsf{\Gamma }}^{\ast }.$$

Here, we only give the lower bound theory of the first-order stationary point of Problem (7). With reference to the definitions of the first- and second-order stationary point established by [23,24] for the general unconstrained $L_q$ regularized optimization, a class of second-order necessary optimality conditions of Problem (7) and its lower bound theory can be considered.

Note that Problem (7) is actually a very challenging optimization problem due to the positive definite constraints and the non-Lipschitz regularization terms, even finding a first-order stationary point. Here, we define a class of approximate first-order stationary points of Problem (7).

Definition 2. 

We say that $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })\in {\mathbb{S}}_{+}^p\times {\mathbb{S}}_{+}^p$ is an ε-approximate first-order stationary point ($\epsilon \ge 0$) of Problem (7) if there exist a set ${\mathsf{\Gamma }}^{\ast }\subseteq \{(u,v):1\le u,v\le p,S_{uv}^{\ast }\ne 0\}$ such that

$$\begin{array}{ccc}& & \hfill \mathrm{dist}\left({\mathit{L}}^{\ast },{\mathrm{Proxi}}_{({\lambda }_L/\ell ){\parallel \mathit{L}\parallel }_q^q,{\mathbb{S}}_{+}^p}\left({\mathit{L}}^{\ast }-\frac{1}{\ell }({\mathit{L}}^{\ast }+{\mathit{S}}^{\ast }-\tilde{\Sigma })\right)\right)\le \epsilon ,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & \hfill \mathrm{dist}\left(\mathit{0},{\mathit{L}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }+{\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }-{\tilde{\Sigma }}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}+{\mathit{Q}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }+{\mathsf{\Xi }}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}:{\mathsf{\Xi }}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\in \partial h_{{\mathsf{\Gamma }}^{\ast }\bigcup D}\left({\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }\bigcup D}^{\ast }\right)\right)\le \epsilon ,\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & \hfill \parallel {\mathit{S}}_{{({\mathsf{\Gamma }}^{\ast }\bigcup D)}^c}^{\ast }{\parallel }_F\le \epsilon .\hfill \end{array}$$

(15)

where ${\mathit{Q}}^{\ast }$ is defined in Definition 1.

Moreover, if $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ is a first-order stationary point of Problem (7), it is an $\epsilon$-approximate first-order stationary point of Problem (7) with $\epsilon =0$. If $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ is an $\epsilon$-approximate first-order stationary point of Problem (7) with $\epsilon =0$, it is a first-order stationary point of Problem (7) with $\epsilon =0$. In practice, the stopping point of any algorithm for an optimization problem is usually an approximate stationary point, rather than a stationary point.

5. Numerical Optimization

In this section, we propose a smoothing alternating updating (SAU) method to solve Problem (7), which is a combination of the smoothing projection gradient step in [28,29] and the proximal gradient step in [30] under the framework of the ordinary alternating minimization method. The convergence of the SAU method is established.

5.1. A Smooth Approximation Optimization Problem

Referring to the smoothing technique in [23,29], a class of smooth approximation optimizations of Problem (7) can be given by

$$\begin{array}{cc}\hfill \underset{\mathit{L},\mathit{S}\in {\mathbb{S}}_{+}^p}{min}& \hfill F_{\mu }(\mathit{L},\mathit{S}):=\frac{1}{2}\parallel \mathit{L}+\mathit{S}-\tilde{\Sigma }{\parallel }_F^2+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q+\sum _{(u,v)\in D^c}{\lambda }_{uv}^S{\psi }_{\mu }{\left(S_{uv}\right)}^q,\hfill \end{array}$$

(16)

where $\mu \in (0,\infty )$ and

$${\psi }_{\mu }\left(S_{uv}\right)=\left\{\begin{array}{cc}|S_{uv}|,\hfill & \hfill |S_{uv}|>\mu ,\hfill \\ \frac{S_{uv}^2}{2\mu }+\frac{\mu }{2},\hfill & \hfill |S_{uv}|\le \mu .\hfill \end{array}\right.$$

Obviously, ${\psi }_{\mu }\left(S_{uv}\right)$ is first-order continuously differential and

$$\begin{array}{c}\hfill {\left({\psi }_{\mu }{\left(S_{uv}\right)}^q\right)}^{\prime }=\left\{\begin{array}{cc}\hfill q|S_{uv}{|}^{q-1}\mathrm{sign}\left(S_{uv}\right),\hfill & \hfill |S_{uv}|>\mu ,\hfill \\ \hfill q{\left(\frac{S_{uv}^2}{2\mu }+\frac{\mu }{2}\right)}^{q-1}\frac{S_{uv}}{\mu },\hfill & \hfill |S_{uv}|\le \mu .\hfill \end{array}\right.\end{array}$$

(17)

Furthermore,

$${\nabla }_S{F}_{\mu }(\mathit{L},\mathit{S})=\mathit{L}+\mathit{S}-\tilde{\Sigma }+\mathsf{\Psi }\left(\mathit{S}\right),$$

where

$$\begin{array}{c}\hfill {\left(\mathsf{\Psi }\left(\mathit{S}\right)\right)}_{uv}=\left\{\begin{array}{cc}{\lambda }_{uv}^S{\left({\psi }_{\mu }{\left(S_{uv}\right)}^q\right)}^{\prime },\hfill & \hfill (u,v)\in D^c,\hfill \\ 0,\hfill & \hfill (u,v)\in D.\hfill \end{array}\right.\end{array}$$

(18)

Next, we define a class of hybrid first-order stationary points of Problem (16).

Definition 3. 

We say that $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })\in {\mathbb{S}}_{+}^p\times {\mathbb{S}}_{+}^p$ is a first-order stationary point of Problem (16) if $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ satisfies

$$\begin{array}{ccc}& & \hfill {\mathit{L}}^{\ast }\in {\mathrm{Proxi}}_{({\lambda }_L/\ell ){\parallel \mathit{L}\parallel }_q^q,{\mathbb{S}}_{+}^p}\left({\mathit{L}}^{\ast }-\frac{1}{\ell }({\mathit{L}}^{\ast }+{\mathit{S}}^{\ast }-\tilde{\Sigma })\right),\hfill \end{array}$$

(19)

$$\begin{array}{ccc}& & \hfill \mathit{0}\in \left\{{\nabla }_S{F}_{\mu }({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })+\mathit{Y},\mathit{Y}\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)\right\}.\hfill \end{array}$$

(20)

Similar to the proof of Theorem 1, it is easy to prove that any global solution of Problem (16) is a first-order stationary point of Problem (16) with $\ell >1$. The following conclusions establish the relationship between the first-order stationary points of Problems (7) and (16), which is the core of the convergence of the SAU method in next subsection.

Theorem 3. 

Let $\left\{{\mu }_k\right\}$ denote a sequence with ${\mu }_k>0$, $k=1,2,\cdots ,$ and ${\mu }_k\to 0$ as $k\to 0$ and $({\mathit{L}}_{{\mu }_k}^{\ast },{\mathit{S}}_{{\mu }_k}^{\ast })$ be a first-order stationary point of Problem (16) with $\mu ={\mu }_k$. Then, the following statements holds:

(i)

If ${\mu }_k\le \epsilon /p^2$, then $({\mathit{L}}_{{\mu }_k}^{\ast },{\mathit{S}}_{{\mu }_k}^{\ast })$ is an ε-approximate first-order stationary point of Problem (7).

(ii)

Any accumulation point of $({\mathit{L}}_{{\mu }_k}^{\ast },{\mathit{S}}_{{\mu }_k}^{\ast })$ is a first-order stationary point of Problem (7).

5.2. Smoothing Alternating Updating Method

Theorem 3 provides a feasible procedure to obtain an approximate solution of Problem (7) by solving Problem (16). This subsection gives the framework of the SAU method for Problem (7) and establishes its convergence.

First, we give the framework of the alternating updating (AU) method for solving Problem (16).

Algorithm 1. 

Given ${\alpha }_0>0$, $\beta \in (0,1)$, $\ell >1$, choose an the initial point $({\mathit{L}}^0,{\mathit{S}}^0)\in {\mathbb{S}}_{+}^p\times {\mathbb{S}}_{+}^p$. For $m=0,1,2,\cdots$, repeat the following two steps until convergence:

Step 1.Update $\mathit{L}$ by solving the following optimization problem:

$$\begin{array}{c}\hfill {\mathit{L}}^{m+1}=\underset{\mathit{L}\in {\mathbb{S}}_{+}^p}{\mathrm{argmin}}\langle {\nabla }_{L}f({\mathit{L}}^m,{\mathit{S}}^m),\mathit{L}-{\mathit{L}}^m\rangle +\frac{\ell }{2}\parallel \mathit{L}-{\mathit{L}}^m{\parallel }_F^2+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q.\end{array}$$

(21)

Step 2.Compute

$$\begin{array}{c}\hfill {\mathit{S}}^{m+1}=\underset{\mathit{S}\in {\mathbb{S}}_{+}^p}{\mathrm{argmin}}{\parallel \mathit{S}-({\mathit{S}}^m-{\alpha }_m{\nabla }_S{F}_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^m))\parallel }_F^2,\end{array}$$

(22)

where ${\alpha }_m={\alpha }_0{\beta }^{{\ell }_m}$ and ${\ell }_m$ is the smallest non-negative integer ℓ such that

$$F_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^m\left({\alpha }_0{\beta }^{\ell }\right)))\le F_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^m))-\frac{c}{2}\parallel {\mathit{S}}^m\left({\alpha }_0{\beta }^{\ell }\right){)-{\mathit{S}}^m\parallel }_F^2,$$

and ${\mathit{S}}^m\left(\alpha \right)={min}_{\mathit{S}\in {\mathbb{S}}_{+}^p}{\parallel \mathit{S}-({\mathit{S}}^m-\alpha {\nabla }_S{F}_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^m))\parallel }_F^2$.

The AU method proposed in this paper applies the proximal gradient step [30] to update $\mathit{L}$ and uses the smoothing projection gradient step [28,29] to update $\mathit{S}$. A key component at each iteration of the proposed AU method is to solve a singular-value minimization Problem (21), which can be rewritten as follows:

$$\begin{array}{c}\hfill \underset{\mathit{L}\in {\mathbb{S}}_{+}^p}{min}\frac{1}{2}\parallel \mathit{L}-\left({\mathit{L}}^m-\frac{1}{\ell }{\nabla }_{L}f({\mathit{L}}^m,{\mathit{S}}^m)\right){\parallel }_F^2+\frac{{\lambda }_L}{\ell }{\parallel \mathit{L}\parallel }_q^q.\end{array}$$

(23)

Theorem 3.4 of [30] provides a solution to Problem (23). To facilitate the expression, define vector q-thresholding operator $\mathcal{H}\left(\mathit{d}\right):={(h_v\left(d_1\right),\cdots ,h_v\left(d_p\right))}^T$ with $v=2{\lambda }_L/\ell$, where

$$\begin{array}{c}\hfill h\left(d_i\right)=\left\{\begin{array}{cc}{h}_{v,q}\left(d_i\right)\hfill & \hfill d_i>d^{\ast },\hfill \\ {\left(v(1-q)\right)}^{1/(2-q)}\mathrm{or}0\hfill & \hfill d_i=d^{\ast },\hfill \\ 0,\hfill & \hfill d_i<d^{\ast },\hfill \end{array}\right.\end{array}$$

(24)

$d^{\ast }=\frac{2-q}{2(1-q)}{\left(v(1-q)\right)}^{1/(2-q)}$, and $h_{v,q}\left(d_i\right)$ is the unique root of the single-variable minimization problem in $(\overline{d},+\infty )$:

$$\underset{x\ge 0}{min}{g}_{d_i}\left(x\right):=x^2-2d_{i}x+v{x}^q,$$

where $\overline{d}={(vq(1-q)/2)}^{1/(2-q)}$. Moreover, $h_{v,q}\left(d_i\right)$ is differentiable and strictly increasing on $[d^{\ast },+\infty )$. As per the suggestion in [25], $h_{v,q}\left(d_i\right)$ can be simply obtained by the numerical methods, for example the Newton method:

$$x^{z+1}=x^z-\frac{qv{\left(x^z\right)}^{q-1}/2+x^z-d_i}{q(q-1)v{\left(x^z\right)}^{q-2}/2+1}$$

with the initial point $d_0=1.5\overline{d}$.

Theorem 4 

(Theorem 3.4 in [30]). Let $\mathit{U}\mathrm{Diag}\left(\mathit{d}\right){\mathit{U}}^T$ be the eigenvalue decomposition of ${\mathit{L}}^m-(1/\ell ){\nabla }_{L}f({\mathit{L}}^m,{\mathit{S}}^m))$. Then, ${\mathit{L}}^{\ast }=\mathit{U}\mathrm{Diag}\left(\mathcal{H}\left(\mathit{d}\right)\right){\mathit{U}}^T$ is an optimal solution to Problem (23), where $\mathcal{H}\left(\mathit{d}\right)$ is defined as (24) with $v=2{\lambda }_L/\ell$.

When $q=1/2$, $\mathcal{H}\left(\mathit{d}\right)$ has a closed-form expression, that is $h_{v,q}\left(d_i\right)=\frac{2}{3}{d}_i(1+cos(\frac{2\pi }{3}-\frac{2\phi \left(d_i\right)}{3}))$, $\phi \left(d_i\right)=arccos\left(\frac{v}{8}{\left(\frac{d_i}{3}\right)}^{-3/2}\right)$, $d^{\ast }=\frac{\sqrt[3]{54}}{4}{v}^{2/3}$. Refer to [30] for the details.

Another key component of at each iteration of the proposed AU method is a projection problem over ${\mathbb{S}}_{+}^p$, which can be written as follows:

$$\begin{array}{c}\hfill \underset{\mathit{S}\in {\mathbb{S}}_{+}^p}{min}{\parallel \mathit{S}-{\mathit{C}}^m\parallel }_F^2,\end{array}$$

(25)

where $C^m:={\mathit{S}}^m-{\alpha }_m{\nabla }_S{F}_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^m)$. It is well known that, if $C^m$ has the eigen-decomposition $C^m:={\sum }_{i=1}^p{\lambda }_i{\mathit{v}}_i{\mathit{v}}_i^T$, then $S^{\ast }={\sum }_{i=1}^{p}max\{{\lambda }_i,0\}{\mathit{v}}_i{\mathit{v}}_i^T$ is the optimal solution of Problem (25).

By (23) and (25), it is easy to find that, when the initial matrices ${\mathit{L}}^0$ and ${\mathit{S}}^0$ in the AU method are symmetric, the matrices for $\mathit{L}$ and $\mathit{S}$ output by the AU method are also symmetric. The next theorem gives the convergence of the AU method.

Theorem 5. 

Let $\left\{({\mathit{L}}^m,{\mathit{S}}^m)\right\}$ be the sequence generated by the AU method. Then, the following statements hold:

(i)

$F_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^{m+1})\le F_{\mu }({\mathit{L}}^m,{\mathit{S}}^m)-\frac{c}{2}\parallel {\mathit{S}}^{m+1}-{\mathit{S}}^m{\parallel }_F^2-\frac{\ell -1}{2}{\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m\parallel }_F^2$.

(ii)

$\left\{({\mathit{L}}^m,{\mathit{S}}^m)\right\}$ is bounded.

(iii)

${lim}_{m\to \infty }{\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m\parallel }_F=0$, ${lim}_{m\to \infty }{\parallel {\mathit{S}}^{m+1}-{\mathit{S}}^m\parallel }_F=0$.

(iv)

Any accumulation point of $\left\{({\mathit{L}}^m,{\mathit{S}}^m)\right\}$ is a first-order stationary point of Problem (16).

We now give the framework of the SAU method for solving Problem (16).

Algorithm 2. 

Give $\tau \in (0,1)$. Choose an initial smoothing parameter ${\mu }_0$. For $k=0,1,2,\cdots$, repeat the following two steps until convergence:

Step 1.Solve Problem (16) with ${\mu }_k$ by the AU method, and obtain $({\mathit{L}}_{{\mu }_k}^{\ast },{\mathit{S}}_{{\mu }_k}^{\ast })$.

Step 2.Set ${\mu }_{k+1}\le \tau {\mu }_k$, and return to Step 1.

Theorem 3 shows that any accumulation point of $({\mathit{L}}_{{\mu }_k}^{\ast },{\mathit{S}}_{{\mu }_k}^{\ast })$ is a first-order stationary point of Problem (7), which provides the guarantee of the convergence of the SAU method. In practice, k cannot be taken as infinite, and ${\mu }_k$ cannot be taken as 0. Then, the $L_q$-CSCE we obtained by the SAU method is often an $\epsilon$-approximate first-order stationary point of Problem (7). Indeed, we stop the SAU method if it reaches a maximum iteration $k_{max}$ or ${\mu }_k<\epsilon$, where $\epsilon >0$ is a small constant. Moreover, in order to improve the efficiency of the SAU method, a trick is to take the output point $({\mathit{L}}_{{\mu }_k}^{\ast },{\mathit{S}}_{{\mu }_k}^{\ast })$ as the initial point of the AU method for Problem (16) with ${\mu }_{k+1}$. Thus, if the initial point of the AU method for solving Problem (16) with ${\mu }_0$, then the final matrices for $\mathit{L}$ and $\mathit{S}$ generated by the SAU method are also symmetric.

6. Simulations

In this section, we study the performance of the $L_q$-CSC method and the associated algorithms. The following models are considered in this section:

• Model 1 (banded matrix with ordering). ${\Sigma }^{\ast }=\mathrm{blockdiag}({\mathit{A}}_1,{\mathit{A}}_2)$, where ${\mathit{A}}_1=\left({\sigma }_{uv}\right)$, $1\le u,v\le p/2$, ${\sigma }_{uv}={(1-\frac{|i-j|}{10})}_{+}$, ${\mathit{A}}_2=4{\mathit{I}}_{p/2\times p/2}$. ${\Sigma }^{\ast }$ is a two-block diagonal matrix, ${\mathit{A}}_1$ is a banded and sparse covariance matrix, and ${\mathit{A}}_2$ is a diagonal matrix with four along the diagonal.
• Model 2 (simple factor). ${\Sigma }^{\ast }=\mathit{U}\mathit{D}{\mathit{U}}^T+\mathit{I}$, where $\mathit{U}\in R^{p\times 3}$ with orthonormal columns uniformly random and $\mathit{D}=8\mathit{I}$.
• Model 3 (random factor). ${\Sigma }^{\ast }=\mathit{U}\mathit{D}{\mathit{U}}^T+{\Sigma }_{\mathrm{Sparse}}$, where $\mathit{U}\in R^{p\times 3}$ with orthonormal columns uniformly random and $\mathit{D}=8\mathit{I}$, where ${\Sigma }_{\mathrm{Sparse}}$ is the sparse matrix generated by Model 1.

The matrix ${\Sigma }^{\ast }$ in Model 1 is sparse, while the matrix in Model 2 and Model 3 satisfies structural decomposition (2).

6.1. Comparison with Representative Methods

This subsection studies the advantages and disadvantages of the $L_q$-CSCE and existing methods. Four representative estimators (the adaptive thresholding estimator (ATE) with the hard thresholding rule in [4], the optimization-type estimator (OTE) in [8], the POET in [12], and the LOREC in [14]) were used for the comparison. The first two estimators are based on the assumption of sparse covariance. While the ATE is a good simple thresholding estimator, the OTE is a good optimization-type estimator. Our estimators POET and LOREC are based on the structural decomposition (2).

Under each model, we set $n=100$, and all random sample are generated from the normal distribution with mean 0 and covariance matrix $\Sigma$, for $p=30,100,200$. The losses are measured by three matrix norms: the spectral norm, the max norm, and the Frobenius norm. The standard errors are given in parentheses. It is known that $q=0.5$ is a good choice for the variable selection problem and the low-rank problem [24,25]. For convenience, the parameter q in Problem (7) is taken as $0.5$. In this subsection, we do not consider different options of q.

For Model 1, the tuning parameters ${\lambda }_L$ in Problems (4) and (7) are taken as ${10}^8$. Moreover, we set $K=0$ in the POET. Under these setting, the LOREC, POET, and $L_q$-CSCE degenerate into the simple soft thresholding estimator [1,2], the ATE, and the constrained adaptive thresholding estimator, respectively.

Table 1. Simulation results for Model 1 over 100 replications. The standard errors are given in parentheses.

p	SE	ATE	OTE	POET	LOREC	${\mathit{L}}_{\mathit{q}}$-CSCE
(Spectral norm)
30	3.76 (0.58)	1.73 (0.50)	1.89 (1.04)	1.74 (0.48)	3.04 (0.43)	1.74 (0.48)
100	8.58 (0.67)	2.69 (0.44)	2.70 (0.50)	2.61 (0.49)	4.94 (0.38)	2.60 (0.50)
200	14.10 (0.82)	4.73 (0.52)	3.27 (0.42)	3.02 (0.45)	6.17 (0.38)	3.02 (0.46)
(Max norm)
30	1.25 (0.23)	1.17 (0.28)	1.16 (0.28)	1.16 (0.28)	1.22 (0.22)	1.16 (0.28)
100	1.59 (0.26)	1.59 (0.26)	1.48 (0.32)	1.51 (0.33)	1.50 (0.25)	1.51 (0.33)
200	1.73 (0.21)	1.73 (0.21)	1.63 (0.25)	1.67 (0.25)	1.63 (0.19)	1.67 (0.25)
(Frobenius norm)
30	7.69 (0.52)	3.17 (0.48)	3.32 (1.03)	3.22 (0.43)	6.35 (0.47)	3.12 (0.44)
100	25.30 (0.64)	8.10 (0.66)	7.24 (0.52)	6.73 (0.49)	16.06 (0.58)	6.41 (0.52)
200	50.50 (0.79)	21.41 (0.80)	12.11 (0.48)	9.80 (0.47)	25.59 (0.62)	9.25 (0.49)

reports the means and standard errors of the five estimators under the three losses. In our analysis, the performance of the POET is almost the same as the ATE. Note that compared with the ATE, OTE, POET, and $L_q$-CSCE, the LOREC has a large bias under the spectral norm loss and the Frobenius norm loss, which is caused by the bias of the $L_1$ penalty. Moreover, the $L_q$-CSCE performs best under the spectral norm loss and the Frobenius norm loss. All methods perform similarly under the max norm.

For Models 2 and 3, we choose the tuning parameters in the LOREC and $L_q$-CSCE by setting the grid and setting $K=3$ in the POET.

Table 2. Simulation results for Model 2 over 100 replications. The standard errors are given in parentheses.

p	SE	ATE	OTE	POET	LOREC	${\mathit{L}}_{\mathit{q}}$-CSCE
(Spectral norm)
30	3.30 (0.75)	3.31 (0.10)	3.28 (0.03)	3.77 (0.79)	3.23 (0.73)	3.27 (0.77)
100	4.99 (0.75)	7.60 (0.08)	7.62 (0.13)	5.48 (0.85)	4.60 (0.62)	4.60 (0.63)
200	6.85 (0.59)	8.36 (0.79)	8.65 (0.49)	7.40 (0.28)	6.30 (0.25)	6.27 (0.63)
(Max norm)
30	0.73 (0.18)	0.74 (0.17)	0.73 (0.18)	0.75 (0.19)	0.72 (0.17)	0.72 (0.16)
100	0.57 (0.08)	0.62 (0.05)	0.62 (0.05)	0.56 (0.09)	0.55 (0.07)	0.50 (0.07)
200	0.54 (0.07)	0.54 (0.07)	0.52 (0.08)	0.54 (0.07)	0.52 (0.07)	0.51 (0.07)
(Frobenius norm)
30	5.59 (0.58)	5.70 (0.59)	5.58 (0.57)	5.47 (0.76)	5.17 (0.58)	4.98 (0.63)
100	12.62 (0.40)	13.58 (0.05)	13.49 (0.24)	10.26 (0.57)	10.06 (0.40)	7.91 (0.48)
200	22.54 (0.39)	14.21 (0.09)	13.94 (0.02)	16.04 (0.53)	18.71 (0.39)	14.37 (0.48)

and

Table 3. Simulation results for Model 3 over 100 replications. The standard errors are given in parentheses.

p	SE	ATE	OTE	POET	LOREC	${\mathit{L}}_{\mathit{q}}$-CSCE
(Spectral norm)
30	5.59 (1.04)	5.57 (0.85)	5.53 (1.03)	5.73 (1.09)	5.57 (1.04)	5.33 (0.99)
100	10.09 (0.77)	9.71 (0.76)	10.01 (0.77)	10.19 (0.78)	9.52 (0.76)	9.21 (0.96)
200	15.00 (0.93)	13.17 (0.84)	14.22 (0.89)	15.69 (0.95)	14.41 (0.93)	11.28 (1.11)
(Max norm)
30	1.63 (0.30)	1.70 (0.28)	1.63 (0.30)	1.65 (0.29)	1.63 (0.30)	1.59 (0.29)
100	1.71 (0.29)	1.71 (0.29)	1.71 (0.29)	1.71 (0.29)	1.67 (0.29)	1.72 (0.37)
200	1.76 (0.22)	1.76 (0.22)	1.75 (0.22)	1.76 (0.22)	1.75 (0.21)	1.60 (0.16)
(Frobenius norm)
30	10.44 (0.88)	11.09 (0.94)	10.34 (0.87)	10.54 (0.91)	10.39 (0.88)	10.10 (0.84)
100	27.99 (0.77)	27.48 (0.76)	27.78 (0.77)	27.76 (0.78)	25.69 (0.76)	19.66 (0.96)
200	52.86 (0.75)	49.25 (0.74)	50.16 (0.73)	49.76 (0.76)	49.09 (0.74)	32.23 (0.78)

report the means and standard errors of the five estimators under the three losses. Note that the ATE and OTE performed poorly for Models 2 and 3. In contrast, the LOREC, POET, and $L_q$-CSCE were more effective. The $L_q$-CSCE performs best under the spectral norm loss and the Frobenius norm loss. The same as the result of Model 1, all methods perform similarly under the max norm. It is worth mentioning that the adaptive technology is invalid for Model 2.

Moreover,

Table 4. Record of the number of positive semi-definite matrices over 100 replications.

Method	$\mathit{p}=30$		$\mathit{p}=100$		$\mathit{p}=200$
	$\mathit{L}$	$\mathit{S}$	$\mathit{L}$	$\mathit{S}$	$\mathit{L}$	$\mathit{S}$
Model1
ATE	-	-	-	-	-	-
OTE	-	-	-	-	-	-
POET	100	0	100	0	100	0
LOREC	100	86	100	0	100	0
$L_q$-CSCE	100	100	100	100	100	100
Model2
ATE	-	-	-	-	-	-
OTE	-	-	-	-	-	-
POET	100	100	100	2	100	0
LOREC	0	100	53	100	100	100
$L_q$-CSCE	100	100	100	100	100	100
Model3
ATE	-	-	-	-	-	-
OTE	-	-	-	-	-	-
POET	100	0	100	0	100	0
LOREC	31	0	0	25	0	0
$L_q$-CSCE	100	100	100	100	100	100

records the number of times that the estimators of the POET, LOREC, and $L_q$-CSCE for $\mathit{L}$ and $\mathit{S}$ are positive semi-definite matrices. Clearly, all estimators of the $L_q$-CSC estimation for $\mathit{L}$ and $\mathit{S}$ are positive semi-definite. It is not guaranteed that the estimators of the LOREC and POET for $\mathit{S}$ are positive semi-definite matrices.

6.2. Choice of q

In the previous simulations, we only consider the setting where $q=0.5$. This subsection supplements some simulations to study the influence of the choice of q on the performance of $L_q$-CSCE. Here, we consider the cases $q=0.1$, $0.3$, $0.5$$0.7$, and $0.9$.

Table 5. Simulation results of the $L_q$-CSCE with different q for Model 1 over 100 replications. The standard errors are given in parentheses.

p	$\mathit{q}=0.1$	$\mathit{q}=0.3$	$\mathit{q}=0.5$	$\mathit{q}=0.7$	$\mathit{q}=0.9$
(Spectral norm)
30	1.70 (0.49)	1.71 (0.48)	1.66 (0.51)	1.68 (0.46)	2.12 (0.43)
100	2.61 (0.54)	2.61 (0.54)	2.62 (0.55)	2.80 (0.46)	3.71 (0.32)
200	3.10 (0.52)	3.11 (0.52)	3.10 (0.53)	3.42 (0.49)	4.10 (0.46)
(Max norm)
30	1.18 (0.31)	1.18 (0.32)	1.13 (0.30)	1.21 (0.32)	1.21 (0.33)
100	1.53 (0.33)	1.53 (0.32)	1.53 (0.34)	1.54 (0.31)	1.52 (0.28)
200	1.65 (0.26)	1.65 (0.27)	1.65 (0.33)	1.66 (0.32)	1.66 (0.30)
(Frobenius norm)
30	3.12 (0.52)	3.15 (0.51)	3.10 (0.54)	3.73 (0.53)	4.45 (0.46)
100	6.38 (0.54)	6.34 (0.54)	6.34 (0.52)	6.43 (0.64)	7.21 (0.55)
200	9.31 (0.53)	9.22 (0.55)	9.21 (0.57)	9.70 (0.51)	11.12 (0.49)

,

Table 6. Simulation results of the $L_q$-CSCE with different q for Model 2 over 100 replications. The standard errors are given in parentheses.

p	$\mathit{q}=0.1$	$\mathit{q}=0.3$	$\mathit{q}=0.5$	$\mathit{q}=0.7$	$\mathit{q}=0.9$
(Spectral norm)
30	3.41 (0.75)	3.41 (0.74)	3.40 (0.72)	3.39 (0.68)	3.38 (0.62)
100	4.81 (0.75)	4.75 (0.10)	4.77 (0.03)	4.74 (0.79)	4.83 (0.79)
200	6.35 (0.75)	6.30 (0.10)	6.31 (0.03)	6.40 (0.79)	6.41 (0.79)
(Max norm)
30	0.73 (0.25)	0.74 (0.18)	0.74 (0.18)	0.73 (0.18)	0.74 (0.19)
100	0.57 (0.11)	0.53 (0.10)	0.53 (0.07)	0.53 (0.04)	0.54 (0.05)
200	0.53 (0.10)	0.50 (0.07)	0.51 (0.09)	0.52 (0.07)	0.54 (0.08)
(Frobenius norm)
30	5.58 (0.62)	5.36 (0.65)	5.33 (0.64)	5.32 (0.61)	5.52 (0.56)
100	8.15 (0.62)	8.09 (0.52)	8.11 (0.50)	8.30 (0.52)	9.14 (0.57)
200	14.81 (0.65)	14.49 (0.57)	14.62 (0.53)	14.74 (0.56)	15.55 (0.49)

and

Table 7. Simulation results of the $L_q$-CSCE with different q for Model 3 over 100 replications. The standard errors are given in parentheses.

p	$\mathit{q}=0.1$	$\mathit{q}=0.3$	$\mathit{q}=0.5$	$\mathit{q}=0.7$	$\mathit{q}=0.9$
(Spectral norm)
30	5.80 (1.01)	5.63 (0.99)	5.51 (0.98)	5.71 (0.95)	5.65 (1.00)
100	8.91 (1.04)	8.89 (0.90)	8.97 (0.93)	8.85 (0.84)	9.21 (0.79)
200	12.40 (1.13)	12.20 (1.04)	12.21 (1.01)	12.34 (1.06)	13.37 (0.99)
(Max norm)
30	1.60 (0.35)	1.56 (0.29)	1.56 (0.28)	1.56 (0.31)	1.60 (0.29)
100	1.73 (0.43)	1.73 (0.31)	1.74 (0.34)	1.76 (0.33)	1.69 (0.31)
200	1.71 (0.29)	1.71 (0.22)	1.72 (0.19)	1.73 (0.19)	1.78 (0.21)
(Frobenius norm)
30	10.12 (0.91)	9.98 (0.87)	10.03 (0.89)	10.14 (0.84)	10.42 (0.86)
100	20.60 (1.03)	20.01 (0.99)	20.04 (1.01)	20.45 (0.95)	21.36 (0.83)
200	31.56 (0.99)	31.40 (0.90)	31.29 (0.84)	31.49 (0.79)	34.47 (0.72)

report the means and standard errors of the $L_q$-CSCE with different q under the three losses. Note that the options for q did not affect the result of the $L_q$-CSCE under the max norm loss. Under other losses, for different models, the performance of the $L_q$-CSCE depends on the choice of q deeply or shallowly. Specifically, there was little difference between $q=0.1$, $q=0.3$, and $q=0.5$ under the spectral norm loss and the max norm loss. The $L_q$-CSCE with $q=0.7$ and $q=0.9$ performed poorly under the Frobenius norm, especially when $q=0.9$. Note that, when $q=0.1$, the $L_q$-CSCE is slightly worse than that of $q=0.3$ and $q=0.5$. We guess that a smaller q would make the minimizing functional more non-convex and, thus, more difficult to solve. The estimator induced by the $L_q$ function with a large q had a large deviation.

Problem (23) has an explicit solution when $q=0.5$, which makes the iteration of the SAU method simple. When $q=0.1$, $0.3$, $0.7$, and $0.9$, it is necessary to embed other numerical methods to find the approximate solution of Problem (23). Therefore, in practical applications, we recommend the setting where $q=0.5$.

7. Discussion and Conclusions

In this paper, we considered the problem of estimating a high-dimensional covariance matrix from noisy data and studied a class of estimators based on the structural decomposition (2). We estimated the high-dimensional covariance matrix by solving a constrained adaptive $L_q$-type regularized optimization Problem (7). This paper provided a systematic study for Problem (7). We first defined a class of hybrid first-order stationary points of Problem (7) and analyzed the relationship between these stationary points and the global solution of Problem (7). Secondly, we proposed an SAU method for solving Problem (7) and established its convergence. The final simulation results showed that the proposed SAU method can solve Problem (7) effectively and that the final $L_q$-CSCE had good numerical performance.

8. Proof

This section provides the detailed proof of the conclusions in this paper.

Proof of Proposition 1. 

$\left(i\right)\left(ii\right)$ The properness of $F(\mathit{L},\mathit{S})$ is obvious. Note that F is continuous over its domain and $\mathrm{dom}\left(F\right)$ is closed. Theorem 2.8 in [31] shows that F is closed. Furthermore, whether ${\parallel \mathit{L}\parallel }_F\to \infty$ or ${\parallel \mathit{S}\parallel }_F\to \infty$, at least one of the three terms in F tends to infinity. Thus, the coerciveness of $F(\mathit{L},\mathit{S})$ holds. Then, it follows from Theorem 2.14 in [31] that $F(\mathit{L},\mathit{S})$ attains its minimal value over ${\mathbb{S}}_{+}^p\times {\mathbb{S}}_{+}^p$.

$\left(iii\right)$ By calculation, $\nabla f(\mathit{L},\mathit{S})={[{\nabla }_{\mathit{L}}^{T}f(\mathit{L},\mathit{S}),{\nabla }_{\mathit{S}}^{T}f(\mathit{L},\mathit{S})]}^T$, where

$${\nabla }_{\mathit{L}}f(\mathit{L},\mathit{S})={\nabla }_{\mathit{S}}f(\mathit{L},\mathit{S})=\mathit{L}+\mathit{S}-\tilde{\Sigma }.$$

For any $({\mathit{L}}_1,{\mathit{S}}_1)$, $({\mathit{L}}_2,{\mathit{S}}_2)\in {\mathbb{R}}^{p\times p}\times {\mathbb{R}}^{p\times p}$, it holds that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \parallel \nabla f({\mathit{L}}_1,{\mathit{S}}_1)-\nabla f({\mathit{L}}_2,{\mathit{S}}_2){\parallel }_F^2& =\parallel {\nabla }_{\mathit{L}}f({\mathit{L}}_1,{\mathit{S}}_1)-{\nabla }_{\mathit{L}}f({\mathit{L}}_2,{\mathit{S}}_2){\parallel }_F^2+{\parallel {\nabla }_{\mathit{S}}f({\mathit{L}}_1,{\mathit{S}}_1)-{\nabla }_{\mathit{S}}f({\mathit{L}}_2,{\mathit{S}}_2)\parallel }_F^2\hfill \\ & =2\parallel {\mathit{L}}_1-{\mathit{L}}_2+{\mathit{S}}_1-{\mathit{S}}_2{\parallel }_F^2\hfill \\ & =2\parallel {\mathit{L}}_1-{\mathit{L}}_2{\parallel }_F^2+2{\parallel {\mathit{S}}_1-{\mathit{S}}_2\parallel }_F^2+4\langle {\mathit{L}}_1-{\mathit{L}}_2,{\mathit{S}}_1-{\mathit{S}}_2\rangle \hfill \\ & \le 4\parallel {\mathit{L}}_1-{\mathit{L}}_2{\parallel }_F^2+4{\parallel {\mathit{S}}_1-{\mathit{S}}_2\parallel }_F^2\hfill \\ & =4\parallel ({\mathit{L}}_1,{\mathit{S}}_1)-({\mathit{L}}_2,{\mathit{S}}_2){\parallel }_F^2,\hfill \end{array}\end{array}$$

which implies that the gradient $\nabla f(\mathit{L},\mathit{S})$ is Lipschitz continuous with a constant $\ell =2$.

$\left(iv\right)$ By calculation, ${\nabla }_{S}f({\mathit{L}}^{\ast },\mathit{S})={\mathit{L}}^{\ast }+\mathit{S}-\tilde{\Sigma }$. For any ${\mathit{S}}_1$, ${\mathit{S}}_2\in {\mathbb{R}}^{p\times p}$, it holds that ${\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}_1)-{\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}_2)={\mathit{S}}_1-{\mathit{S}}_2$. Then, the partial gradient function ${\nabla }_{S}f({\mathit{L}}^{\ast },\mathit{S})$ is Lipschitz continuous with a constant $\ell =1$. The same result holds for the partial gradient function ${\nabla }_{L}f(\mathit{L},{\mathit{S}}^{\ast })$. □

Proof of Theorem 1. 

Since $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ is a global minimizer of Problem (7), then ${\mathit{L}}^{\ast }$ is also a global minimizer of the problem:

$$\begin{array}{c}\hfill \underset{\mathit{L}\in {\mathbb{S}}_{+}^p}{min}f(\mathit{L},{\mathit{S}}^{\ast })+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q.\end{array}$$

By $\left(iv\right)$ of Proposition 1, we have that

$$f({\mathit{L}}_1,\mathit{S})\le f({\mathit{L}}_2,\mathit{S})+\langle {\nabla }_{L}f({\mathit{L}}_2,\mathit{S}),{\mathit{L}}_1-{\mathit{L}}_2\rangle +\frac{1}{2}{\parallel {\mathit{L}}_1-{\mathit{L}}_2\parallel }_F^2,\forall {\mathit{L}}_1,{\mathit{L}}_2\in {\mathbb{R}}^{p\times p}.$$

It follows that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })+{\lambda }_L{\parallel {\mathit{L}}^{\ast }\parallel }_q^q& \le f(\mathit{L},{\mathit{S}}^{\ast })+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q\hfill \\ & \le f(\mathit{L},{\mathit{S}}^{\ast })+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q+\frac{\ell -1}{2}{\parallel \mathit{L}-{\mathit{L}}^{\ast }\parallel }_F^2\hfill \\ & \le f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })+\langle {\nabla }_{L}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }),\mathit{L}-{\mathit{L}}^{\ast }\rangle +\frac{1}{2}\parallel \mathit{L}-{\mathit{L}}^{\ast }{\parallel }_F^2+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q+\frac{\ell -1}{2}{\parallel \mathit{L}-{\mathit{L}}^{\ast }\parallel }_F^2\hfill \\ & =f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })+\langle {\nabla }_{L}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }),\mathit{L}-{\mathit{L}}^{\ast }\rangle +\frac{\ell }{2}\parallel \mathit{L}-{\mathit{L}}^{\ast }{\parallel }_F^2+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q.\hfill \end{array}\end{array}$$

Thus,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{L}}^{\ast }& \in \underset{\mathit{L}\in {\mathbb{S}}_{+}^p}{\mathrm{argmin}}\left\{f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })+\langle {\nabla }_{L}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }),\mathit{L}-{\mathit{L}}^{\ast }\rangle +\frac{\ell }{2}\parallel \mathit{L}-{\mathit{L}}^{\ast }{\parallel }_F^2+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q\right\}\hfill \\ & ={\mathrm{Proxi}}_{({\lambda }_L/\ell ){\parallel \mathit{L}\parallel }_q^q,{\mathbb{S}}_{+}^p}\left({\mathit{L}}^{\ast }-\frac{1}{\ell }({\mathit{L}}^{\ast }+{\mathit{S}}^{\ast }-\tilde{\Sigma })\right).\hfill \end{array}\end{array}$$

Hence, (9) holds.

We now prove that (10) holds for $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$. Observe that ${\mathit{S}}^{\ast }$ is a global minimizer of the problem:

$$\begin{array}{c}\hfill \underset{\mathit{S}\in {\mathbb{S}}_{+}^p}{min}f({\mathit{L}}^{\ast },\mathit{S})+\sum _{(u,v)\in D^c}{\lambda }_{uv}^S{\left|S_{uv}\right|}^q,\end{array}$$

which can be equivalently rewritten as

$$\begin{array}{c}\hfill \underset{\mathit{S}\in {\mathbb{R}}^{p\times p}}{min}f({\mathit{L}}^{\ast },\mathit{S})+\sum _{(u,v)\in D^c}{\lambda }_{uv}^S{\left|S_{uv}\right|}^q+{\delta }_{{\mathbb{S}}_{+}^p}\left(\mathit{S}\right).\end{array}$$

Then, ${\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }}^{\ast }$ is a global minimizer of the problem:

$$\begin{array}{c}\hfill \underset{{\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }}}{min}f({\mathit{L}}^{\ast },{\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }},{\mathit{S}}_{{\left({\mathsf{\Gamma }}^{\ast }\right)}^c}^{\ast })+\sum _{(u,v)\in D^c\bigcap {\mathsf{\Gamma }}^{\ast }}{\lambda }_{uv}^S{\left|S_{uv}\right|}^q+{\delta }_{{\mathbb{S}}_{+}^p}\left([{\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }}:{\mathit{S}}_{{\left({\mathsf{\Gamma }}^{\ast }\right)}^c}^{\ast }]\right).\end{array}$$

(26)

It follows from the first-order optimality condition of Problem (26) and Exercise 10.10 in [32] that there exists a ${\mathit{Y}}_{{\mathsf{\Gamma }}^{\ast }}\in \partial h_{{\mathsf{\Gamma }}^{\ast }}\left({\mathit{S}}_{{\mathsf{\Gamma }}^{\ast }}\right)$ such that

$$\begin{array}{ccc}& & \hfill 0=L_{uv}^{\ast }+{\mathit{S}}_{uv}^{\ast }-{\tilde{\Sigma }}_{uv}+{\lambda }_{uv}^{S^{\ast }}\mathrm{sign}\left(S_{uv}^{\ast }\right){\left|S_{uv}^{\ast }\right|}^{q-1}+Y_{uv},\mathrm{for}\mathrm{all}(u,v)\in {\mathsf{\Gamma }}^{\ast }\bigcap D^c,\hfill \end{array}$$

(27)

$$\begin{array}{ccc}& & \hfill 0=L_{uv}^{\ast }+{\mathit{S}}_{uv}^{\ast }-{\tilde{\Sigma }}_{uv}+Y_{uv},\mathrm{for}\mathrm{all}(u,v)\in {\mathsf{\Gamma }}^{\ast }\bigcap D.\hfill \end{array}$$

(28)

Moreover, let $\mathit{Y}=[{\mathit{Y}}_{{\mathsf{\Gamma }}^{\ast }}:{\mathit{Y}}_{{\left({\mathsf{\Gamma }}^{\ast }\right)}^c}]$, where ${\mathit{Y}}_{{\left({\mathsf{\Gamma }}^{\ast }\right)}^c}$ is some finite array. For any $(u,v)\in {\left({\mathsf{\Gamma }}^{\ast }\right)}^c\bigcap D^c$, we have

$$\begin{array}{c}\hfill S_{uv}^{\ast }(L_{uv}^{\ast }+{\mathit{S}}_{uv}^{\ast }-{\tilde{\Sigma }}_{uv})+{\lambda }_{uv}^S{\left|S_{uv}\right|}^q+S_{uv}^{\ast }{Y}_{uv}=0,\end{array}$$

(29)

and for any $(u,v)\in {\left({\mathsf{\Gamma }}^{\ast }\right)}^c\bigcap D$, we have

$$\begin{array}{c}\hfill S_{uv}^{\ast }(L_{uv}^{\ast }+{\mathit{S}}_{uv}^{\ast }-{\tilde{\Sigma }}_{uv})+S_{uv}^{\ast }{Y}_{uv}=0,\end{array}$$

(30)

Combining (27)-(30) with the definition of $\mathit{Y}$, (10) holds. Thus, $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ is a first-order stationary point of Problem (7). □

Proof of Theorem 2. 

Since $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ is a first-order stationary point of Problem (7) and ${\mathit{S}}^{\ast }\in {\mathbb{S}}_{++}^p$, by Definition 1, it holds that

$$\begin{array}{c}\hfill \mathit{0}={\mathit{S}}^{\ast }\circ ({\mathit{L}}^{\ast }+{\mathit{S}}^{\ast }-\tilde{\Sigma })+{\mathit{Q}}^{\ast },\end{array}$$

(31)

where ${\mathit{Q}}^{\ast }:={\left(Q_{uv}^{\ast }\right)}_{u,v=1}^p$ and

$$\begin{array}{c}\hfill Q_{uv}^{\ast }=\left\{\begin{array}{cc}q{\lambda }_{uv}^S{\left|S_{uv}^{\ast }\right|}^q,\hfill & \mathrm{if}(u,v)\in D^c,\hfill \\ 0,\hfill & \mathrm{if}(u,v)\in D.\hfill \end{array}\right.\end{array}$$

For all $(u,v)\in D^c\bigcap {\mathsf{\Gamma }}^{\ast }$, by (31), we have

$$\begin{array}{c}\hfill 0=S_{uv}^{\ast }(L_{uv}^{\ast }+S_{uv}^{\ast }-{\tilde{\Sigma }}_{uv})+q{\lambda }_{uv}^S{\left|S_{uv}^{\ast }\right|}^q.\end{array}$$

(32)

Moreover, the statement $\left(iv\right)$ in Proposition 1 shows that ${\nabla }_{S}f(\mathit{L},\mathit{S})$ is Lipschitz continuous with a constant 1, and it follows that

$$f({\mathit{L}}^{\ast },\mathit{Y})\le f({\mathit{L}}^{\ast },\mathit{X})+\langle {\nabla }_{S}f({\mathit{L}}^{\ast },\mathit{X}),\mathit{Y}-\mathit{X}\rangle +\frac{1}{2}{\parallel \mathit{Y}-\mathit{X}\parallel }_F^2,\forall \mathit{X},\mathit{Y}\in {\mathbb{R}}^{p\times p}.$$

Taking $\mathit{X}={\mathit{S}}^{\ast }$ and $\mathit{Y}={\mathit{S}}^{\ast }-{\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$, we obtain that

$$f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }-{\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }))\le f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })-\frac{1}{2}{\parallel {\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })\parallel }_F^2.$$

Note that

$$f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }-{\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }))\ge 0,f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })\le F({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })\le F({\mathit{L}}^0,{\mathit{S}}^0)+ϵ.$$

It follows that

$$\parallel {\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }){\parallel }_F\le \sqrt{2[F({\mathit{L}}^0,{\mathit{S}}^0)+ϵ]}.$$

This, together with the relation (32), implies that

$$q{\lambda }_{uv}^S|S_{uv}^{\ast }{|}^{q-1}=|{\left({\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })\right)}_{uv}|\le \parallel {\nabla }_{S}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }){\parallel }_F\le \sqrt{2[F({\mathit{L}}^0,{\mathit{S}}^0)+ϵ]},\forall (u,v)\in D^c\bigcap {\mathsf{\Gamma }}^{\ast }.$$

Thus, we have

$$|S_{uv}^{\ast }|\ge {\left(\frac{qmin\left\{{\lambda }_{uv}^S\right\}}{\sqrt{2[F({\mathit{L}}^0,{\mathit{S}}^0)+ϵ]}}\right)}^{\frac{1}{1-q}}\forall (u,v)\in D^c\bigcap {\mathsf{\Gamma }}^{\ast },$$

which completes this proof. □

Proof of Theorem 3. 

$\left(i\right)$ We just need to prove that Relation (20) is a sufficient condition for (14) and the inequality (15) holds. Since $({\mathit{L}}_{{\mu }_k}^{\ast },{\mathit{S}}_{{\mu }_k}^{\ast })$ is a first-order stationary point of Problem (16), there exists a matrix ${\mathit{Y}}_{{\mu }_k}^{\ast }\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)$ such that

$$\mathit{0}={\nabla }_S{F}_{{\mu }_k}\left({\mathit{L}}_{{\mu }_k}^{\ast }\right)+{\mathit{Y}}_{{\mu }_k}^{\ast }.$$

Let $\mathcal{N}:=\{(u,v)\in D^c:{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}>\mu \}$. It holds that

$$\begin{array}{c}\hfill \mathit{0}={\left({\mathit{L}}_{{\mu }_k}^{\ast }\right)}_{\mathcal{N}\bigcup D}+{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{\mathcal{N}\bigcup D}-{\tilde{\Sigma }}_{\mathcal{N}\bigcup D}+{\left[\mathsf{\Psi }\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)\right]}_{\mathcal{N}\bigcup D}+{\left({\mathit{Y}}_{{\mu }_k}^{\ast }\right)}_{\mathcal{N}\bigcup D}.\end{array}$$

(33)

By (11), it is known that

$${\left({\mathit{Y}}_{{\mu }_k}^{\ast }\right)}_{\mathcal{N}\bigcup D}\in {\mathcal{L}}_{\mathcal{N}\bigcup D}^{\ast }\left(N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)\right)\subseteq \partial h_{\mathcal{N}\bigcup D}\left({\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{\mathcal{N}\bigcup D}\right).$$

Then,

$$\mathrm{dist}\left(\mathit{0},{\left({\mathit{L}}_{{\mu }_k}^{\ast }\right)}_{\mathcal{N}\bigcup D}+{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{\mathcal{N}\bigcup D}-{\tilde{\Sigma }}_{\mathcal{N}\bigcup D}+{\left[\mathsf{\Psi }\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)\right]}_{\mathcal{N}\bigcup D}+{\left({\mathit{Y}}_{{\mu }_k}\right)}_{\mathcal{N}\bigcup D}:{\mathit{Y}}_{\mathcal{N}\bigcup D}\in \partial h_{\mathcal{N}\bigcup D}\left({\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{\mathcal{N}\bigcup D}\right)\right)\le \epsilon .$$

Moreover, it follows from the condition that $\mu \le \epsilon /p^2$ that

$$\parallel {\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{{({\mathsf{\Gamma }}^{\ast }\bigcup D)}^c}{\parallel }_F\le \sum _{(u,v)\in {({\mathsf{\Gamma }}^{\ast }\bigcup D)}^c}\left|{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}\right|\le \mu p^2\le \epsilon .$$

Therefore, $({\mathit{L}}_{{\mu }_k}^{\ast },{\mathit{S}}_{{\mu }_k}^{\ast })$ is an $\epsilon$-approximate first-order stationary point of Problem (7).

$\left(ii\right)$ Let $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ be an accumulation point of $({\mathit{L}}_{{\mu }_k}^{\ast },{\mathit{S}}_{{\mu }_k}^{\ast })$. By the definition of the first-order stationary point of Problem (16), for any $(u,v)\in {\mathcal{N}}^c\bigcap D^c$, there exists a matrix ${\mathit{Y}}_{{\mu }_k}^{\ast }\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)$ such that

$$\begin{array}{c}\hfill 0={\left({\mathit{L}}_{{\mu }_k}^{\ast }\right)}_{uv}+{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}-{\tilde{\Sigma }}_{uv}+q{\lambda }_{uv}{\left(\frac{{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}^2}{2{\mu }_k}+\frac{{\mu }_k}{2}\right)}^{q-1}\frac{{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}}{{\mu }_k}+{\left(Y_{{\mu }_k}^{\ast }\right)}_{uv}.\end{array}$$

(34)

Note that

$$\begin{array}{c}\hfill 0\le q{\lambda }_{uv}{\left(\frac{{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}^2}{2{\mu }_k}+\frac{{\mu }_k}{2}\right)}^{q-1}\frac{{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}^2}{{\mu }_k}\le q{\lambda }_{uv}{\left(\frac{{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}^2}{{\mu }_k}\right)}^q\le q{\lambda }_{uv}^S{\left|{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}\right|}^q.\end{array}$$

(35)

If $|{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}|\le {\mu }_k$ for arbitrarily large k, then $S_{uv}^{\ast }={lim}_{k\to \infty }{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}=0$, and by (35), we have

$$\underset{k\to \infty }{lim}q{\lambda }_{uv}{\left(\frac{{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}^2}{2{\mu }_k}+\frac{{\mu }_k}{2}\right)}^{q-1}\frac{{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}^2}{{\mu }_k}=0\mathrm{and}q{\lambda }_{uv}^S{\left|S_{uv}^{\ast }\right|}^q=0.$$

Since ${\mathit{Y}}_{{\mu }_k}^{\ast }\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)$, ${lim}_{k\to \infty }{\mathit{Y}}_{{\mu }_k}^{\ast }\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)$.

This, together with the relation (34), implies that, for any $(u,v)\in {\mathcal{N}}^c\bigcap D^c$,

$$\begin{array}{c}\hfill 0\in \left\{S_{uv}^{\ast }(L_{uv}^{\ast }+S_{uv}^{\ast }-{\tilde{\Sigma }}_{uv})+q{\lambda }_{uv}^S{\left|S_{uv}^{\ast }\right|}^q+S_{uv}^{\ast }{Y}_{uv},\mathit{Y}\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)\right\}.\end{array}$$

(36)

Moreover, for any $(u,v)\in \mathcal{N}\bigcup D$, we have

$$\begin{array}{c}\hfill 0={\left({\mathit{L}}_{{\mu }_k}^{\ast }\right)}_{uv}+{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}-{\tilde{\Sigma }}_{uv}+Q_{uv}^{{\mu }_k}+{\left(Y_{{\mu }_k}^{\ast }\right)}_{uv},\end{array}$$

(37)

where

$$Q_{uv}^{{\mu }_k}=\left\{\begin{array}{cc}0,\hfill & (u,v)\in D,\hfill \\ q{\lambda }_{uv}^S\mathrm{sign}\left({\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}\right){\left|{\left({\mathit{S}}_{{\mu }_k}^{\ast }\right)}_{uv}\right|}^{q-1},\hfill & (u,v)\in (\mathcal{N}\bigcup D)/D.\hfill \end{array}\right.$$

For any $(u,v)\in \mathcal{N}\bigcup D$, it follows that

$$\begin{array}{c}\hfill 0\in \left\{S_{uv}^{\ast }(L_{uv}^{\ast }+S_{uv}^{\ast }-{\tilde{\Sigma }}_{uv})+S_{uv}^{\ast }{Q}_{uv}+S_{uv}^{\ast }{Y}_{uv},\mathit{Y}\in N_{{\mathbb{S}}_{+}^p}\left({\mathit{S}}^{\ast }\right)\right\},\end{array}$$

(38)

where

$$Q_{uv}=\left\{\begin{array}{cc}0,\hfill & (u,v)\in D,\hfill \\ q{\lambda }_{uv}^S\mathrm{sign}\left({\mathit{S}}_{uv}^{\ast }\right){\left|{\mathit{S}}_{uv}^{\ast }\right|}^{q-1},\hfill & (u,v)\in (\mathcal{N}\bigcup D)/D.\hfill \end{array}\right.$$

Since the $Y_{{\mu }_k}^{\ast }$ in (34) coincides with that in (37), by (36) and (38), we have that $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ is the first-order stationary point of Problem (16). □

Proof of Theorem 5. 

$\left(i\right)$ By $\left(iv\right)$ of Proposition 1, we have that

$$\frac{1}{2}\parallel {\mathit{L}}_1+\mathit{S}-\tilde{\Sigma }{\parallel }_F^2\le \frac{1}{2}\parallel {\mathit{L}}_2+\mathit{S}-\tilde{\Sigma }{\parallel }_F^2+\langle {\mathit{L}}_2+\mathit{S}-\tilde{\Sigma },{\mathit{L}}_1-{\mathit{L}}_2\rangle +\frac{1}{2}{\parallel {\mathit{L}}_1-{\mathit{L}}_2\parallel }_F^2,\forall {\mathit{L}}_1,{\mathit{L}}_2\in {\mathbb{R}}^{p\times p}.$$

It follows that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill F_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^m)& =\frac{1}{2}\parallel {\mathit{L}}^{m+1}+{\mathit{S}}^m-\tilde{\Sigma }{\parallel }_F^2+{\lambda }_L{\parallel {\mathit{L}}^{m+1}\parallel }_q^q+\sum _{(u,v)\in D^c}{\lambda }_{uv}^S{\psi }_{\mu }{\left(S_{uv}^m\right)}^q,\hfill \\ \hfill & \le \frac{1}{2}\parallel {\mathit{L}}^m+{\mathit{S}}^m-\tilde{\Sigma }{\parallel }_F^2+\langle {\mathit{L}}^m+{\mathit{S}}^m-\tilde{\Sigma },{\mathit{L}}^{m+1}-{\mathit{L}}^m\rangle +\frac{1}{2}\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m{\parallel }_F^2+{\lambda }_L{\parallel {\mathit{L}}^{m+1}\parallel }_q^q\hfill \\ \hfill & \hfill +\sum _{(u,v)\in D^c}{\lambda }_{uv}^S{\psi }_{\mu }{\left(S_{uv}^m\right)}^q\hfill \\ \hfill & \le \frac{1}{2}\parallel {\mathit{L}}^m+{\mathit{S}}^m-\tilde{\Sigma }{\parallel }_F^2+\langle {\mathit{L}}^m+{\mathit{S}}^m-\tilde{\Sigma },{\mathit{L}}^{m+1}-{\mathit{L}}^m\rangle +\frac{\ell }{2}\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m{\parallel }_F^2+{\lambda }_L{\parallel {\mathit{L}}^{m+1}\parallel }_q^q\hfill \\ \hfill & \hfill +\sum _{(u,v)\in D^c}{\lambda }_{uv}^S{\psi }_{\mu }{\left(S_{uv}^m\right)}^q-\frac{\ell -1}{2}{\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m\parallel }_F^2\hfill \\ \hfill & \le \frac{1}{2}\parallel {\mathit{L}}^m+{\mathit{S}}^m-\tilde{\Sigma }{\parallel }_F^2+{\lambda }_L\parallel {\mathit{L}}^m{\parallel }_q^q+\sum _{(u,v)\in D^c}{\lambda }_{uv}^S{\psi }_{\mu }{\left(S_{uv}^m\right)}^q-\frac{\ell -1}{2}{\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m\parallel }_F^2\hfill \\ \hfill & =F_{\mu }({\mathit{L}}^m,{\mathit{S}}^m)-\frac{\ell -1}{2}{\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m\parallel }_F^2,\hfill \end{array}\end{array}$$

where the third inequality is due to the iteration of $\mathit{L}$ in AU. Then,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill F_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^{m+1})& \le F_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^m)-\frac{c}{2}{\parallel {\mathit{S}}^{m+1}-{\mathit{S}}^m\parallel }_F^2\hfill \\ \hfill & \le F_{\mu }({\mathit{L}}^m,{\mathit{S}}^m)-\frac{c}{2}\parallel {\mathit{S}}^{m+1}-{\mathit{S}}^m{\parallel }_F^2-\frac{\ell -1}{2}{\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m\parallel }_F^2.\hfill \end{array}\end{array}$$

$\left(ii\right)$ It is easy to show that $F_{\mu }$ is coercive, which implies that $F_{\mu }$ is level bounded. By $\left(i\right)$ in this Theorem, we have $F_{\mu }({\mathit{L}}^m,{\mathit{S}}^m)\le F_{\mu }({\mathit{L}}^0,{\mathit{S}}^0)$. Then, $\left\{({\mathit{L}}^m,{\mathit{S}}^m)\right\}$ is bounded.

$\left(iii\right)$ It follows from the inequality in $\left(i\right)$ that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \sum _{m=0}^{\infty }\left(\frac{c}{2}\parallel {\mathit{S}}^{m+1}-{\mathit{S}}^m{\parallel }_F^2+\frac{\ell -1}{2}{\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m\parallel }_F^2\right)& \hfill \le \sum _{m=0}^{\infty }\left(F_{\mu }({\mathit{L}}^m,{\mathit{S}}^m)-F_{\mu }({\mathit{L}}^{m+1},{\mathit{S}}^{m+1})\right)\hfill \\ \hfill & \hfill \le F_{\mu }({\mathit{L}}^0,{\mathit{S}}^0)-\underset{m\to \infty }{lim}{F}_{\mu }({\mathit{L}}^m,{\mathit{S}}^m)<\infty ,\hfill \end{array}\end{array}$$

where the last inequality is due to $F_{\mu }$ being bounded from below. Hence, ${lim}_{m\to \infty }{\parallel {\mathit{L}}^{m+1}-{\mathit{L}}^m\parallel }_F=0$, ${lim}_{m\to \infty }{\parallel {\mathit{S}}^{m+1}-{\mathit{S}}^m\parallel }_F=0$.

$\left(iv\right)$ Let $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ be an accumulation point of $({\mathit{L}}^m,{\mathit{S}}^m)$. We assumed that ${lim}_{j\to \infty }{\mathit{L}}^{m_j}={\mathit{L}}^{\ast }$ and ${lim}_{j\to \infty }{\mathit{S}}^{m_j}={\mathit{S}}^{\ast }$, where $m_j\to \infty$ as $j\to \infty$. By the results in $\left(iii\right)$,

$$\underset{j\to \infty }{lim}{\mathit{L}}^{m_j+1}=\underset{j\to \infty }{lim}{\mathit{L}}^{m_j}+({\mathit{L}}^{m_j+1}-{\mathit{L}}^{m_j})={\mathit{L}}^{\ast },\underset{j\to \infty }{lim}{\mathit{S}}^{m_j+1}=\underset{j\to \infty }{lim}{\mathit{S}}^{m_j}+({\mathit{S}}^{m_j+1}-{\mathit{S}}^{m_j})={\mathit{S}}^{\ast }.$$

By the AU method in this paper, we have that

$$\begin{array}{ccc}& & \hfill {\mathit{L}}^{m_j+1}=\underset{\mathit{L}\in {\mathbb{S}}_{+}^p}{\mathrm{argmin}}\langle {\nabla }_{L}f({\mathit{L}}^{m_j},{\mathit{S}}^{m_j}),\mathit{L}-{\mathit{L}}^{\ast }\rangle +\frac{\ell }{2}\parallel \mathit{L}-{\mathit{L}}^{m_j}{\parallel }_F^2+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q,\hfill \\ & & \hfill \langle {\mathit{S}}^{m_j+1}-({\mathit{S}}^{m_j}-{\alpha }_{m_j}{\nabla }_S{F}_{\mu }({\mathit{L}}^{m_j},{\mathit{S}}^{m_j})),\mathit{S}-{\mathit{S}}^{m_j+1}\rangle \ge 0,\forall \mathit{S}\in {\mathbb{S}}_{+}^p.\hfill \end{array}$$

Taking the limits on both sides, it follows that

$$\begin{array}{ccc}& & {\mathit{L}}^{\ast }=\underset{\mathit{L}\in {\mathbb{S}}_{+}^p}{\mathrm{argmin}}\langle {\nabla }_{L}f({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }),\mathit{L}-{\mathit{L}}^{\ast }\rangle +\frac{\ell }{2}\parallel \mathit{L}-{\mathit{L}}^{m_j}{\parallel }_F^2+{\lambda }_L{\parallel \mathit{L}\parallel }_q^q,\hfill \\ & & \langle {\nabla }_S{F}_{\mu }({\mathit{L}}^{\ast },{\mathit{S}}^{\ast }),\mathit{S}-{\mathit{S}}^{\ast }\rangle \ge 0,\forall \mathit{S}\in {\mathbb{S}}_{+}^p,\hfill \end{array}$$

i.e., $({\mathit{L}}^{\ast },{\mathit{S}}^{\ast })$ is a first-order stationary point of Problem (16). □