Sparse and Low-Rank Covariance Matrices Estimation

1 Introduction

Estimation of population covariance matrices from samples of multivariate data has draw many attentions in the last decade owing to its fundamental importance in multivariate analysis. With dramatic advances in technology in recent years, various research fields, such as genetic data, brain imaging, spectroscopic imaging, climate data and so on, have been used to deal with massive high-dimensional data sets, whose sample sizes can be very small relative to dimension. In such settings, the standard and the most usual sample covariance matrices often performs poorly

[1, 2, 11].

Fortunately, regularization as a class of new methods to estimate covariance matrices has recently emerged to overcome those shortages of using traditional sample covariance matrices. These methods encompass several specified forms, banding [1, 6, 17], tapering [4, 10] and thresholding [2, 5, 8, 16] for instance.

Moreover, there are many cases where the model is known to be structured in several ways at the same time. In recent years, one of research contents is to estimate a covariance matrix possessing both sparsity and positive definiteness. For instance, Rothman [15] gave the following model:

{^Σ}_{+} = argmin Σ ≻ 0 \frac{1}{2} {∥ Σ - Σ_{n} ∥}_{n}^{2} - τ log det (Σ) + λ {∥ ∥ Σ^{-} ∥ ∥}_{1}^{-},

(1)

where $Σ_{n}$ is the sample covariance matrix, $∥ \cdot ∥_{F}$ is the Frobenious norm, $∥ \cdot ∥_{1}$ is the element-wise $ℓ_{1}$ -norm and $Σ−:=Σ−\textmdDiag(Σ)$ . From the optimization viewpoint, (1) is similar to the graphical lasso criterion [9] which also has a log-determinant part and the element-wise $ℓ_{1}$ -penalty. Rothman [15] derived an iterative procedure to solve (1). While Xue, Ma and Zou [18] omitted the log-determinant part and considered the positive definite constraint ${Σ⪰ε\emph{I}}$ for some arbitrarily small $ε > 0$ :

^Σε=argminΣ⪰ε\emph{I} 12∥Σ−Σn∥2F+λ∥∥Σ−∥∥1,

(2)

They utilized an efficient alternating direction method (ADM) to solve the challenging problem (2) and established its convergence properties.

Most of the literatures, e.g.,[12, 15, 18], required the population covariance matrices being positive definite, and thus there is no essence of pursuing the low-rank of the estimator. By contrast, newly appeared research topic is to consider simultaneously the sparsity and low-rank of a structured model, which implies that the population covariance matrices are no longer restricted to the positive definite matrix cone and can be relaxed to the positive semidefinite cone. In addition, the models with structure of being simultaneously the sparsity and low-rank are widely applied into practice, such as sparse signal recovery from quadratic measurements and sparse phase retrieval, see [13] for example. Moreover, Richard et al. [14]

showed that both sparse and low-rank model can be derived in covariance matrix when the random variables are highly correlated in groups, which means this covariance matrix has a block diagonal structure.

With stimulations of those ideas, we construct the following convex model encompassing the $ℓ_{1}$ -norm and nuclear norm for estimating the covariance matrix:

^Σ = argmin Σ ⪰ 0 \frac{1}{2} {∥ Σ - Σ_{n} ∥}_{n}^{2} + λ {∥ Σ ∥}_{1} + τ {∥ Σ ∥}_{*},

(3)

where $λ \geq 0, τ \geq 0$ are tuning parameters. The $ℓ_{1}$ -norm penalty ${∥ Σ ∥}_{1} = \sum_{i, j} | σ_{i j} |$ is also called lasso-type penalty and is used to encourage sparse solutions. The nuclear norm , ${∥ Σ ∥}_{*} = \sum_{i} | λ_{i} (Σ) |$ with $λ_{i} (Σ)$

being the eigenvalue of

Σ

, is the trace norm when

Σ ⪰ 0

and ensures low-rank solutions of (3). Here we inroduce the approximate rank to interpret the low-rank, which is defined as

a r (A) := a r \in {1, 2, \dots, p}

being the smallest number such that

\frac{σ_{a r + 1} (A)}{σ_{1} (A)} \leq γ,

(4)

where $σ_{i} (A) (i = 1, 2, \dots, p)$

are the singular value of

A

with

σ_{1} (A) \geq σ_{2} (A) \geq \dots \geq σ_{p} (A)

, and

γ > 0

could be chosen based on the needs, throughout our paper we fix

γ = 0.001

for simplicity.

The contributions of this paper mainly center on two aspects. For one thing, being different from [13, 14], we establish the theoretical statistical theory under different assumptions rather than giving the generalized error bound of the estimation. Especially, we acquire the estimation rate $O (\sqrt{(s log r) / n})$ under the Frobenious norm error, which improves the optimal rate $O (\sqrt{(s log p) / n})$ where the low-rank property of the estimator does not be considered [7, 15, 18] and $p$ is the samples’ dimension with $p > max {n, r}$ . For another, we take advantage of the alternating direction method of multipliers (ADMM), also can be seen in [18, 19], to combat our problem (3).

The organization of this paper is as follows. In Section 2 we will present some theoretical properties of the estimator derived by the proposed model (3). After that the alternating direction method of multipliers (ADMM) is going to be introduced to combat the problem, and numerical experiments are projected to show the performance of this method in Sections 3 and 4 respectively. We make a conclusion in the last section.

2 A Sparse and Low-Rank Covariance Estimator

Before the main part, we hereafter introduce some notations. $E (X)$ and $P (A)$ denote the expectation of $X$ and the probability of the incident $A$ occurring respectively. $\textmdCard(S)$ is the number of entries of the set $S$

. Normal distribution with mean

μ

and covariance

Σ

is written as

N (μ, Σ)

. Say

Y_{n} = O_{P} (1)

if for every

ε > 0

, there is a

C > 0

such that

P {| Y_{n} | > C} < ε

for all

n \geq n_{0} (ε)

, and say

Y_{n} = O_{P} (a_{n})

Y_{n} / a_{n} = O_{P} (1)

. If there are two constants

C_{1} \leq C_{2}

such that

C_{1} \leq X_{n} / Y_{n} \leq C_{2}

, we write as

X_{n} ≍ Y_{n}

For given observed independently and identically distributed (i.i.d. for short) $p$ -variate random variables $X_{1}, \dots,$ $X_{n}$ with covariance matrix $Σ_{0}$ and $p > n$ , the goal is to estimate the unknown matrix $Σ_{0}$ based on the sample ${X_{l} : l = 1, \dots, n}$ . This problem is called covariance matrix estimation which is of fundamental importance in multivariate analysis.

Given a random sample ${X_{1}, \dots, X_{n}}$ from $E (X) = 0$ (without loss of generality) and a population covariance matrix $Σ_{0} = {(σ_{0 i j})}_{1 \leq i, j \leq p} = E (X X^{⊤})$ , the sample covariance matrix is

Σ_{n} = {(σ_{n i j})}_{1 \leq i, j \leq p} = \frac{1}{n - 1} n \sum l = 1 (X_{l} - ¯ X) {(X_{l} - ¯ X)}_{l}^{⊤},

where $¯ X = \frac{1}{n} \sum_{l = 1}^{n} X_{l}$ . Denote $S$ the support set of the population covariance matrix $Σ_{0}$ as

S={(i,j):σ0ij≠0}, s=\textmdCard(S), r=\textmdrank(Σ0).

Assumption 2.1

For all $p$ , $0≤λ\emphmin(Σ0)≤λ\emphmax(Σ0)≤¯λ<∞$ , where $¯ λ$ is a constant.

Assumption 2.1 is a common used condition in covariance matrix estimation, on which a useful lemma based is recalled here for the sequel analysis. One can also refer it in [1].

Lemma 2.2

Let $X_{l} \in R^{p}, l = 1, \dots, n$ be i.i.d. $N (0, Σ_{0})$ and Assumption 2.1 holds, (i.e., for all $p$ , $λ\emphmax(Σ0)≤¯λ<∞$ ). Then, if $Σ_{0} = {(σ_{0 i j})}_{1 \leq i, j \leq p}$ ,

P{∣∣ ∣∣n∑l=1(XliXlj−σ0ij)∣∣ ∣∣≥nν}≤C1\emphexp{−C2nν2}, for |ν|≤δ

(5)

where constants $C_{1}, C_{2}$ and $δ$ depend on $¯ λ$ only.

Another lemma which plays an important role in our main results is stated below.

Lemma 2.3

Suppose that Assumption 2.1 holds, $λ \leq \frac{ε}{16 \sqrt{s}}, τ \leq \frac{ε}{8 \sqrt{r}}$ . Then for $ε > 0$ sufficiently small,

\emph{max}1≤i,j≤p∣∣σ0ij−σnij∣∣≤λ  implies  ∥^Σ−Σ0∥F≤ε,

(6)

where $^Σ$ is defined as $(???)$ .

Proof First make the eigenvalue decomposition of $Σ_{0}$ ( $r=\textmdrank(Σ0)$ ) as

Σ0=UΛΣ0U⊤≡U(\textmdDiag(λ(Σ0))0)U⊤,

where $U \in R^{p \times p}$ with $U U^{⊤} = U^{⊤} U = I_{p}$

is the matrix composed of eigenvectors,

\textmdDiag(λ(Σ0))∈Rr×r

is a diagonal matrix generated by eigenvalues with

λ (Σ_{0}) = (λ_{1} (Σ_{0}), \dots, λ_{r} (Σ_{0}))^{⊤}

and

λ_{j} (Σ_{0}) > 0, j = 1, \dots, r

.
By denoting

Σ := U Δ U^{⊤} + Σ_{0}

with

Δ = (\begin{matrix} Δ_{1} & 0 0 & 0 \end{matrix})

and

Δ_{1} \in R^{r \times r}

, which implies

Δ = U^{⊤} Σ U - Λ_{Σ_{0}}

, we consider the model

	$^Δ$	$:=$	$\textmdargminΔ+ΛΣ0⪰0 F(Δ)=\textmdargminΔ1+\textmdDiag(λ(Σ0))⪰0 F(Δ)$
		$\equiv$	$\textmdargminΔ+ΛΣ0⪰0 12∥∥UΔU⊤+Σ0−Σn∥∥2F+λ∥∥UΔU⊤+Σ0∥∥1+τ∥∥UΔU⊤+Σ0∥∥∗.$

Clearly, from (3) we have $^Σ = U^Δ U^{⊤} + Σ_{0}$ which implies $^Δ = U^{⊤}^Σ U - Λ_{Σ_{0}}$ . For a given $ε > 0$ sufficiently small and any $∥ Δ ∥_{F} = ε$ (i.e., $∥ Δ_{1} ∥_{F} = ε$ ), we compute

$F (Δ) - F (0)$	$=$	$\frac{1}{2} {∥ ∥ U Δ U^{⊤} + Σ_{0} - Σ_{n} ∥ ∥}_{0}^{⊤} + λ {∥ ∥ U Δ U^{⊤} + Σ_{0} ∥ ∥}_{1}^{⊤} + τ {∥ ∥ U Δ U^{⊤} + Σ_{0} ∥ ∥}_{*}^{⊤}$
		$- \frac{1}{2} {∥ Σ_{0} - Σ_{n} ∥}_{0}^{2} + λ {∥ Σ_{0} ∥}_{1} + τ {∥ Σ_{0} ∥}_{*}$
	$=$	$\frac{1}{2} ∥ Δ ∥_{F}^{2} + ⟨ U Δ U^{⊤}, Σ_{0} - Σ_{n} ⟩ + λ (∥ U Δ U^{⊤} + Σ_{0} ∥_{1} - ∥ Σ_{0} ∥_{1})$
		$+ τ (∥ U Δ U^{⊤} + Σ_{0} ∥_{} - ∥ Σ_{0} ∥_{})$
	$\equiv$	$12∥Δ∥2F+\textmdI+\textmdII+\textmdIII.$

For convenience we denote $¯ ¯¯¯ ¯ Δ = U Δ U^{⊤}$ . Then for I, it holds

	$\textmdI=⟨¯¯¯¯¯Δ,Σ0−Σn⟩$	$\geq$	$−∣∣\textmdtr(¯¯¯¯¯Δ(Σn−Σ0))∣∣$
		$=$	$- ∣ ∣ \sum i, j {¯ ¯¯¯ ¯ Δ}_{i j} (σ_{0 i j} - σ_{n i j}) ∣ ∣ \geq - max 1 \leq i, j \leq p ∣ ∣ σ_{0 i j} - σ_{n i j} ∣ ∣ ∥ ¯ ¯¯¯ ¯ Δ ∥_{1},$

For II, we obtain by noting $S = {(i, j) : σ_{0 i j} \neq 0}$ and $s=\textmdCard(S)$ that,

II	$=$	$λ (∥ Σ_{0} + ¯ ¯¯¯ ¯ Δ ∥_{1} - {∥ Σ_{0} ∥}_{1})$
	$=$	$λ (∥ Σ_{0 S} + {¯ ¯¯¯ ¯ Δ}_{S} ∥_{1} + ∥ {¯ ¯¯¯ ¯ Δ}_{S^{C}} ∥_{1} - {∥ Σ_{0 S} ∥}_{1})$
	$\geq$	$λ (∥ Σ_{0 S} + {¯ ¯¯¯ ¯ Δ}_{S} ∥_{1} + ∥ {¯ ¯¯¯ ¯ Δ}_{S^{C}} ∥_{1} - (∥ Σ_{0 S} + {¯ ¯¯¯ ¯ Δ}_{S} ∥_{1} + ∥ {¯ ¯¯¯ ¯ Δ}_{S} ∥_{1}))$
	$=$	$λ (∥ {¯ ¯¯¯ ¯ Δ}_{S^{C}} ∥_{1} - ∥ {¯ ¯¯¯ ¯ Δ}_{S} ∥_{1}) .$

From the H $¨ o$ lder Inequality, one can prove that

	$∥ ¯ ¯¯¯ ¯ Δ ∥_{} = ∥ Δ ∥_{} = ∥ Δ_{1} ∥_{*} = r \sum i = 1 \| λ_{i} (Δ_{1}) \| \leq \sqrt{r} ∥ Δ_{1} ∥_{F} = \sqrt{r} ∥ Δ ∥_{F},$		(8)
	$∥ {¯ ¯¯¯ ¯ Δ}_{S} ∥_{1} \leq \sqrt{s} ∥ {¯ ¯¯¯ ¯ Δ}_{S} ∥_{F} .$		(9)

For III, combining with (8) we get that

Since ${max}_{i, j} ∣ ∣ σ_{0 i j} - σ_{n i j} ∣ ∣ \leq λ$ and (9),

$G (Δ)$	$\equiv$	$F (Δ) - F (0)$
	$=$	$12∥Δ∥2F+\textmdI+\textmdII+\textmdIII$
	$\geq$	$\frac{1}{2} {∥ Δ ∥}_{F}^{2} - λ ∥ ¯ ¯¯¯ ¯ Δ ∥_{1} + λ (∥ {¯ ¯¯¯ ¯ Δ}_{S^{C}} ∥_{1} - ∥ {¯ ¯¯¯ ¯ Δ}_{S} ∥_{1}) - τ \sqrt{r} ∥ Δ ∥_{F}$
	$=$	$\frac{1}{2} {∥ Δ ∥}_{F}^{2} - λ (∥ ¯ Δ_{S} ∥_{1} + ∥ {¯ ¯¯¯ ¯ Δ}_{S^{C}} ∥_{1} - (∥ {¯ ¯¯¯ ¯ Δ}_{S^{C}} ∥_{1} - ∥ {¯ ¯¯¯ ¯ Δ}_{S} ∥_{1})) - τ \sqrt{r} ∥ Δ ∥_{F}$
	$=$	$\frac{1}{2} {∥ Δ ∥}_{F}^{2} - 2 λ ∥ {¯ ¯¯¯ ¯ Δ}_{S} ∥_{1} - τ \sqrt{r} ∥ Δ ∥_{F}$
	$\geq$	$\frac{1}{2} {∥ Δ ∥}_{F}^{2} - 2 λ \sqrt{s} ∥ {¯ ¯¯¯ ¯ Δ}_{S} ∥_{F} - τ \sqrt{r} ∥ Δ ∥_{F}$
	$\geq$	$\frac{1}{2} {∥ Δ ∥}_{F}^{2} - 2 λ \sqrt{s} ∥ ¯ ¯¯¯ ¯ Δ ∥_{F} - τ \sqrt{r} ∥ Δ ∥_{F}$
	$=$	$\frac{1}{2} {∥ Δ ∥}_{F}^{2} - 2 λ \sqrt{s} ∥ Δ ∥_{F} - τ \sqrt{r} ∥ Δ ∥_{F} .$

Therefore, by $λ \leq \frac{ε}{16 \sqrt{s}}, τ \leq \frac{ε}{8 \sqrt{r}}$

G (Δ) \geq \frac{1}{2} {∥ Δ ∥}_{F}^{2} - 2 λ \sqrt{s} ∥ Δ_{S} ∥_{F} - τ \sqrt{r} {∥ Δ ∥}_{F} \geq \frac{ε^{2}}{2} - \frac{ε^{2}}{8} - \frac{ε^{2}}{8} = \frac{ε^{2}}{4} > 0.

Hence we prove that if $λ \leq \frac{ε}{16 \sqrt{s}}, τ \leq \frac{ε}{8 \sqrt{r}}$ and ${max}_{i, j} ∣ ∣ σ_{0 i j} - σ_{n i j} ∣ ∣ \leq λ$ , for any $∥ Δ ∥_{F} = ε$ , it holds

G (Δ) \equiv F (Δ) - F (0) > 0.

In addition, from (2), we have

^Δ=\textmdargminΔ+\textmdDiag(λΣ0)⪰0  F(Δ)−F(0) =\textmdargminΔ+\textmdDiag(λΣ0)⪰0  G(Δ),

(10)

which implies that $∥^Δ ∥_{F} \leq ε$ . Otherwise, we suppose $∥^Δ ∥_{F} > ε$ , then $∥ ε^Δ / ∥^Δ ∥_{F} ∥_{F} = ε$ . Since for any $∥ Δ ∥_{F} = ε$ , it follows $G (Δ) > 0 = G (0)$ which is contradicted with the fact $G (\cdot)$ is a convex function and $G (^Δ) \leq G (0) = 0$ , because

	$0 < G (\frac{ε}{∥^Δ ∥_{F}}^Δ)$	$=$	$G (\frac{ε}{∥^Δ ∥_{F}}^Δ + (1 - \frac{ε}{∥^Δ ∥_{F}}) 0)$
		$\leq$	$\frac{ε}{∥^Δ ∥_{F}} G (^Δ) + (1 - \frac{ε}{∥^Δ ∥_{F}}) G (0) \leq 0.$

Finally $∥^Δ ∥_{F} \leq ε$ indicates that $∥^Σ - Σ_{0} ∥_{F} \leq ε$ due to $∥^Δ ∥_{F} = ∥ U^Δ U^{⊤} ∥_{F} = ∥^Σ - Σ_{0} ∥_{F}$ . Hence the desired result is obtained.. $□$

Then in order to acquiring the rate of the estimation, the two following commonly used assumptions are needed to introduced, and also can be seen [1, 15]. Assumption 2.4 holds, for example, if $X_{l 1}, \dots, X_{l p}$ are Gaussian.

Assumption 2.4

$E{\emphexp(tX2lj)}≤C1<∞$ hold for all $j = 1, \dots, p$ and $0 < | t | < t_{0}$ , where $t_{0}$ and $C_{1}$ are two constants.

Assumption 2.5

$E {{∣ ∣ X_{l j} ∣ ∣}_{l j}^{2 α}} \leq C_{2} < \infty$ hold for all $j = 1, \dots, p$ , where some $α \geq 2$ and $C_{2}$ is a constant.

Built on the two assumptions, we give our main results with regard to rates of the estimator of (3).

Theorem 2.6

For some $δ \in [1, 2)$ , let $K_{1}$ be a sufficiently large constant and suppose Assumptions 2.1 and 2.4 hold. If $λ = K_{1} \sqrt{\frac{log r}{n} + \frac{log p}{n K_{1}^{δ}}}, τ = O_{P} (\sqrt{\frac{s log r}{n r} + \frac{s log p}{n r K_{1}^{δ}}})$ and $s log r + (s log p) / K_{1}^{δ} = o (n)$ , then

∥^Σ - Σ_{0} ∥_{F} = O_{P} ⎛ ⎝ \sqrt{\frac{s log r}{n} + \frac{s log p}{n K_{1}^{δ}}} ⎞ ⎠ .

(11)

Proof Since Assumptions 2.1 and 2.4 hold, a fact employed by Rothman et al. [16] is that for $ν > 0$ sufficiently small,

P {max 1 \leq i, j \leq n ∣ ∣ σ_{0 i j} - σ_{n i j} ∣ ∣ \geq ν} \leq C_{3} p^{2} exp {- C_{4} n ν^{2}},

where $C_{3}$ and $C_{4}$ are some constants. We then apply the bound $s log r + (s log p) / K_{1}^{δ} = o (n)$ and Lemma 2.3 with

ε = 16 K_{1} \sqrt{\frac{s log r}{n} + \frac{s log p}{n K_{1}^{δ}}}, λ = K_{1} \sqrt{\frac{log r}{n} + \frac{log p}{n K_{1}^{δ}}}

to obtain that

P{∥^Σ−Σ0∥F≤ε}≥P{max1≤i,j≤n∣∣σ0ij−σnij∣∣≤λ}≥1−C3p2−C4K2−δ1r−C4K21.

Evidently, $1 - C_{3} p^{2 - C_{4} K_{1}^{2 - δ}} r^{- C_{4} K_{1}^{2}}$ can be arbitrarily close to one by choosing $K_{1}$ sufficiently large.. $□$

Corollary 2.7

For some $δ \in [1, 2)$ , let $K_{1}$ be a sufficiently large constant such that $K_{1}^{δ} log r ≍ log p$ , and suppose Assumptions 2.1, 2.4 hold. If $λ = K_{1} \sqrt{\frac{log r}{n}}, τ = O_{P} (\sqrt{\frac{s log r}{n r}})$ and $s log r = o (n)$ , then

∥^Σ - Σ_{0} ∥_{F} = O_{P} (\sqrt{\frac{s log r}{n}}) .

(12)

Remark 2.8

Clearly, if the $λ_{min} (Σ_{0}) > 0$ in Assumption 2.1, the better rate $O_{P} (\sqrt{\frac{s log r}{n}})$ would reduce to $O_{P} (\sqrt{\frac{s log p}{n}})$ . It is worth mentioning that under the the Assumption 2.1, the minimax optimal rate of convergence under the Frobenius norm in Theorem 4 of [7] is $O_{P} (\sqrt{\frac{s log p}{n}})$ which also has been obtained by [18]. However, to attain the same rate in the presence of the log-determinant barrier term (1), Rothman [15] instead would require that $λ_{min}$ , the minimal eigenvalue of the true covariance matrix, should be bounded away from zero by some positive constant, and also that the barrier parameter should be bounded by some positive quantity. [18] illustrated this theory requiring a lower bound on $λ_{min}$ is not very appealing.

Theorem 2.9

Let $K_{2}$ be a sufficiently large constant and suppose that Assumptions 2.1 and 2.5 hold. If $λ = K_{2} \sqrt{\frac{p^{4 / α}}{n}}, τ = O_{P} (\sqrt{\frac{s p^{4 / α}}{n r}})$ and $s p^{4 / α} = o (n)$ , then

∥^Σ - Σ_{0} ∥_{F} = O_{P} ⎛ ⎝ \sqrt{\frac{s p^{4 / α}}{n}} ⎞ ⎠ .

(13)

Proof Since Assumptions 2.1 and 2.5 hold, one can modify a result of Bickel & Levina (2008a) and show that for $ν > 0$ sufficiently small

P {max 1 \leq i, j \leq n ∣ ∣ σ_{0 i j} - σ_{n i j} ∣ ∣ \geq ν} \leq p^{2} C_{5} n^{- α / 2} ν^{- α},

where $C_{5}$ is a constant. We then apply the bound $s p^{4 / α} = o (n)$ and Lemma 2.3 with $ε = 16 K_{2} \sqrt{\frac{s p^{4 / α}}{n}}, λ = K_{2} \sqrt{\frac{p^{4 / α}}{n}}$ to get

P{∥∥^Σ−Σ0∥∥F≤ε}≥P{max% 1≤i,j≤n∣∣σ0ij−σnij∣∣≤λ}≥1−C6K−α2.

Apparently, the bound $1 - C_{6} K_{2}^{- α}$ can be arbitrarily close to one by taking $K_{2}$ sufficiently large.. $□$

3 Alternating Direction Method of multipliers

In this section, we will construct the alternating direction method of multipliers (ADMM) to solve problem (3). By introducing an auxiliary variable $Γ$ , problem (3) can be rewritten as

	$\textmdmin 12∥Σ−Σn∥2F+λ∥Σ∥1+τ∥Γ∥∗,$		(14)
	$\textmds.t. Γ⪰0, Σ−Γ=0.$

The constraint $Γ ⪰ 0$ can be put into the objective function by using an indicator function:

	$I (Γ ⪰ 0) =$	$0, Γ ⪰ 0$
	$I (Γ ⪰ 0) =$	$+∞, \textmdotherwise.$

This leads to the following equivalent reformulation of (14):

	$\textmdmin 12∥Σ−Σn∥2F+λ∥Σ∥1+τ∥Γ∥∗+I(Γ⪰0),$		(15)
	$\textmds.t. Σ−Γ=0.$

Recently, the alternating direction method of multipliers (ADMM) has been studied extensively for solving (15). A typical iteration of ADMM for solving (15) can be described as

		$Γk+1:=\textmdargminΓ Lμ(Σk,Γ,Λk)$			(16)
		$Σk+1:=\textmdargminΣ Lμ(Σ,Γk+1,Λk)$			(17)
		$Λ^{k + 1} := Λ^{k} - \frac{1}{μ} (Γ^{k + 1} - Σ^{k + 1}),$			(18)

where the augmented Lagrangian function $L_{μ} (Σ, Γ, Λ)$ is defined as

	$L_{μ} (Σ, Γ, Λ)$	$:=$			(19)
			$+ I (Γ ⪰ 0) - ⟨ Λ, Γ - Σ ⟩ + \frac{1}{2 μ} ∥ Σ - Γ ∥_{F}^{2},$		(19)

in which $Λ$ is the Lagrange multiplier and $μ > 0$ is a penalty parameter. Note that ADMM (16-18) can be written explicitly as

		$Γk+1:=\textmdargminΓ⪰0 τ∥Γ∥∗+12μ∥Γ−(Σk+μΛk)∥2F$			(20)
		$Σk+1:=\textmdargminΣ λ∥Σ∥1+μ+12μ∥Σ−μμ+1(Σn+1μΓk+1−Λk)∥2F$			(21)
		$Λ^{k + 1} := Λ^{k} - (Γ^{k + 1} - Σ^{k + 1}) / μ .$			(22)

We now show that the two subproblems (20) and (21) can be easily solved. For the subproblem (20),

$Γ^{k + 1}$	$=$	$\textmdargminΓ⪰0 τ∥Γ∥∗+12μ∥Γ−(Σk+μΛk)∥2F$	(23)
	$=$	$\textmdargminΓ⪰0 τ⟨Γ,I⟩+12μ∥Γ−(Σk+μΛk)∥2F$
	$=$	$\textmdargminΓ⪰0 ∥Γ−(Σk+μΛk−μτI)∥2F$
	$=$	$(Σ^{k} + μ Λ^{k} - μ τ I)_{+},$

where $(X)_{+}$ denote the projection of a matrix $X$ onto the convex positive semidefinite cone $S_{+}^{n}$ . Namely $(X)+=U\textmdDiag(max{λ(X),0})U⊤$ , where $X=U\textmdDiag(λ(X))U⊤$ and $λ (X) = (λ_{1} (X), λ_{2} (X), \dots, λ_{p} (X))^{⊤} .$

The solution of the second subproblem (21) is given by the $ℓ_{1}$ -shrinkage operation

	$Σ^{k + 1}$	$=$	$μμ+1\textmdShrink(Σn+1μΓk+1−Λk,λ)$		(24)
		$=$	$μμ+1(\textmdmax {∣∣Pij∣∣−λ,0}\textmdsign(Pij))1≤i,j≤p,$		(24)

where $P = {(P_{i j})}_{1 \leq i, j \leq p} := Σ_{n} + \frac{1}{μ} Γ^{k + 1} - Λ^{k}$ and $\textmdsign(⋅)$ is a sign function.

Therefore, combining with (22)-(24), the whole algorithm is written as follows

ADMM: Alternating Direction Method of Multipliers
Initialize $μ, λ, τ, Σ^{0}, Λ^{0};$
Repeat
Compute $Γ^{k + 1} = (Σ^{k} + μ Λ^{k} - μ τ I)_{+}$ ;
Compute $Σk+1=μμ+1\textmdShrink(Σn+1μΓk+1−Λk,λ);$
Compute $Λ^{k + 1} := Λ^{k} - \frac{1}{μ} (Γ^{k + 1} - Σ^{k + 1}) .$
Untill Convergence
Return

Table 1: The framework of the ADMM.

To end this section, we prove that the sequence $(Γ^{k}, Σ^{k}, Λ^{k})$ produced by the alternating direction method of multipliers (Table 1) converges to $(^Γ,^Σ,^Λ)$ , where $(^Γ,^Σ)$ is an optimal solution of (15) and $^Λ$ is the optimal dual variable. Now we label some necessary notations for the ease of presentation. Let $H$ be a $2 p \times 2 p$ matrix defined as

H=(1μ\emphIp×p00μ\emphIp×p),

the weighted norm ${∥ \cdot ∥}_{H}^{2}$ stands for ${∥ V ∥}_{H}^{2} := ⟨ V, H V ⟩$ and the corresponding inner product $⟨ \cdot, \cdot ⟩_{H}$ is $⟨ U, V ⟩_{H} := ⟨ U, H V ⟩$ . Before presenting the main theorem with regard to the global convergence of ADMM, we introduce the following lemma.

Lemma 3.1

Assume that $(^Γ,^Σ)$ is an optimal solution of $(???)$ and $^Λ$ is the corresponding optimal dual variable associated with the equality constraint $Σ - Γ = 0$ . Then the sequence ${(Γ^{k}, Σ^{k}, Λ^{k})}$ produced by ADMM satisfies

∥^V - V^{k} ∥_{H}^{2} - ∥^V - V^{k + 1} ∥_{H}^{2} \geq ∥ V^{k + 1} - V^{k} ∥_{H}^{2},

(25)

where $V^{k} = (Λ^{k}, Σ^{k})^{⊤}$ and $^V = (^Λ,^Σ)^{⊤}$ .

Based on the lemma above, the convergent theorem can be derived immediately.

Theorem 3.2

The sequence ${(Γ^{k}, Σ^{k}, Λ^{k})}$ generated by Algorithm 1 from any starting point converges to an optimal solution of $(???)$ .

4 Numerical Simulations

In this section we will exploit the proposed method ADMM to tackle two examples, one of which possessed the block structured population covariance matrix, and another utilized the banded population covariance matrix. Actually as the constraint $Σ ⪰ 0$ , our proposed model (3) is equivalent to

(26)

So similar to the method in [18], one can solve the soft-thresholding estimator

Σ_{s t} := argmin Σ \frac{1}{2} {∥ Σ - (Σ_{n} - τ I) ∥}_{F}^{2} + λ {∥ Σ ∥}_{1}

to initialize the $Σ^{0}$ . If the derived $Σ_{s t} ⪰ 0$ then the recovered sparse and low-rank semidefinite estimator $^Σ = Σ_{s t}$ . In our stimulation, we uniformly initialize $Λ^{0}$ as the matrix with all entries being 1, $Σ^{0}$

as zero matrix and

Σ_{s t}

respectively. Unlike

λ

and

τ

μ

does not change the final covariance estimator, thus we fixed

μ = 1

just for simplicity and the stop criteria is set as

max {∥ Γ^{k + 1} - Γ^{k} ∥_{F}, ∥ Σ^{k + 1} - Σ^{k} ∥_{F}} \leq 5 \times 10^{- 4} .

For the sample dimensions, we always take $n = 50$ and $p = 100, 200, 500, 1000$ .

4.1 Example I: Block Structure

Analogous to the model, modified slightly here, emerged in [14] who synthesized $n$ samples $X_{l} \sim N (0, Σ_{0})$ for a block diagonal population covariance matrix $Σ_{0} \in R^{p \times p}$ , we will use $K (= 5, 10, 20)$ blocks of random sizes, and each block is generated by $v v^{⊤}$ where the entries of $v$

are drawn i.i.d. from the uniform distribution on

[- 1, 1]

. Evidently, the rank of

Σ_{0}

produced in the way is

K

. Corresponding MATLAB code of generating

X_{l}

X_{l} = Σ_{0}^{1 / 2} r a n d n (p, 1), l = 1, 2, \dots, n

, thereby deriving the sample covariance matrix

Σ_{n} = \frac{1}{n - 1} n \sum l = 1 (X_{l} - ¯ X) {(X_{l} - ¯ X)}_{l}^{⊤} .

What is worth mentioning is that if $Σ_{0}$ is a positive definite matrix, the solution our method obtains would be a positive definite matrix with full rank. But fortunately, compared to some largest singular values of $Σ_{0}$ , the left are relatively small

Sparse and Low-Rank Covariance Matrices Estimation

A finite sample estimator for large covariance matrices

Square-root nuclear norm penalized estimator for panel data models with approximately low-rank unobserved heterogeneity

Convex Optimization without Projection Steps

Low-Rank Covariance Function Estimation for Multidimensional Functional Data

A Sparsity Inducing Nuclear-Norm Estimator (SpINNEr) for Matrix-Variate Regression in Brain Connectivity Analysis

Simultaneously sparse and low-rank abundance matrix estimation for hyperspectral image unmixing

Low rank prior and l0 norm to remove impulse noise in images

1 Introduction

2 A Sparse and Low-Rank Covariance Estimator

Assumption 2.1

Lemma 2.2

Lemma 2.3

Assumption 2.4

Assumption 2.5

Theorem 2.6

Corollary 2.7

Remark 2.8

Theorem 2.9

3 Alternating Direction Method of multipliers

Lemma 3.1

Theorem 3.2

4 Numerical Simulations

4.1 Example I: Block Structure

Sparse and Low-Rank Covariance Matrices Estimation

Related Research

A finite sample estimator for large covariance matrices

Square-root nuclear norm penalized estimator for panel data models with approximately low-rank unobserved heterogeneity

Convex Optimization without Projection Steps

Low-Rank Covariance Function Estimation for Multidimensional Functional Data

A Sparsity Inducing Nuclear-Norm Estimator (SpINNEr) for Matrix-Variate Regression in Brain Connectivity Analysis

Simultaneously sparse and low-rank abundance matrix estimation for hyperspectral image unmixing

Low rank prior and l0 norm to remove impulse noise in images

1 Introduction

2 A Sparse and Low-Rank Covariance Estimator

Assumption 2.1

Lemma 2.2

Lemma 2.3

Assumption 2.4

Assumption 2.5

Theorem 2.6

Corollary 2.7

Remark 2.8

Theorem 2.9

3 Alternating Direction Method of multipliers

Lemma 3.1

Theorem 3.2

4 Numerical Simulations

4.1 Example I: Block Structure