Diagonally-Dominant Principal Component Analysis

TL;DR

This paper introduces Diagonally-Dominant PCA (DD-PCA), a novel matrix decomposition method that enhances covariance estimation and multiple testing procedures for strongly correlated data, with demonstrated practical benefits.

Contribution

The paper develops DD-PCA, a new decomposition technique, and applies it to improve covariance estimation and global null hypothesis testing in correlated data scenarios.

Findings

DD-PCA provides accurate covariance and precision matrix estimation.

New tests outperform existing methods in multiple testing scenarios.

Algorithms for efficient DD-PCA computation are proposed.

Abstract

We consider the problem of decomposing a large covariance matrix into the sum of a low-rank matrix and a diagonally dominant matrix, and we call this problem the "Diagonally-Dominant Principal Component Analysis (DD-PCA)". DD-PCA is an effective tool for designing statistical methods for strongly correlated data. We showcase the use of DD-PCA in two statistical problems: covariance matrix estimation, and global detection in multiple testing. Using the output of DD-PCA, we propose a new estimator for estimating a large covariance matrix with factor structure. Thanks to a nice property of diagonally dominant matrices, this estimator enjoys the advantage of simultaneous good estimation of the covariance matrix and the precision matrix (by a plain inversion). A plug-in of this estimator to linear discriminant analysis and portfolio optimization yields appealing performance in real data. We also propose two new tests for testing the global null hypothesis in multiple testing when the $z$-scores have a factor covariance structure. Both tests first use DD-PCA to adjust the individual $p$-values and then plug in the adjusted $p$-values to the Higher Criticism (HC) test. These new tests significantly improve over the HC test and compare favorably with other existing tests. For computation of DD-PCA, we propose an iterative projection algorithm and an ADMM algorithm.

Tables

Table 1. Table 1: Performance of the proposed ADMM in Experiment 1.

Dimension $p$	rank($𝑳$)	rank($\widehat{𝑳}$)	$\frac{{\Vert \widehat{𝑳}+\widehat{𝑨}-𝑺\Vert }_F}{{\Vert 𝑺\Vert }_F}$	$\frac{{\Vert \widehat{𝑳}-𝑳\Vert }_F}{{\Vert 𝑳\Vert }_F}$	$\frac{{\Vert \widehat{𝑨}-𝑨\Vert }_F}{{\Vert 𝑨\Vert }_F}$
500	25	25	0.008	0.011	0.045
1000	50	50	0.010	0.008	0.034
2000	100	100	0.013	0.006	0.026

Table 2. Table 2: Performance of Algorithm  1 in Experiment 2.

Dimension $p$	rank($𝑳$)	rank($\widehat{𝑳}$)	$\frac{{\Vert \widehat{𝑳}+\widehat{𝑨}-𝑺\Vert }_F}{{\Vert 𝑺\Vert }_F}$	$\frac{{\Vert \widehat{𝑳}-𝑳\Vert }_F}{{\Vert 𝑳\Vert }_F}$	$\frac{{\Vert \widehat{𝑨}-𝑨\Vert }_F}{{\Vert 𝑨\Vert }_F}$
500	25	25	0.264	0.166	0.340
1000	50	50	0.269	0.163	0.286
2000	100	100	0.274	0.160	0.243

Table 3. Table 3: Estimation errors of DD-POET and its robustness to a misspecified K 𝐾 K .

$(p,K)$	Target	Norm	$k$						POET∗ $\left(k=3\right)$
			1	2	3	4	5	6
(100,3)	${𝚺}_u$	Frobenius	48.26	27.52	3.28	3.46	3.63	3.83	3.24
		Spectral	22.51	16.80	0.80	1.00	1.08	1.14	0.85
	${𝚺}_u^{-1}$	Frobenius	9.10	7.50	3.02	3.17	3.34	3.56	3.65
		Spectral	1.34	1.34	0.61	0.72	0.83	0.93	0.64
(300,3)	${𝚺}_u$	Frobenius	95.00	56.96	6.22	6.23	6.26	6.32	6.04
		Spectral	32.52	26.04	0.82	0.86	0.90	0.93	0.90
	${𝚺}_u^{-1}$	Frobenius	41.50	17.33	5.68	5.66	5.66	5.68	6.37
		Spectral	27.75	7.16	0.66	0.66	0.65	0.64	0.64
(500,3)	${𝚺}_u$	Frobenius	126.50	75.74	8.38	8.35	8.35	8.35	7.99
		Spectral	38.10	30.86	0.84	0.87	0.89	0.91	0.95
	${𝚺}_u^{-1}$	Frobenius	25.87	18.19	7.66	7.61	7.57	7.55	8.26
		Spectral	9.23	1.76	0.69	0.68	0.68	0.67	0.64

Equations

$$\boldsymbol{\Sigma}=\boldsymbol{L}+\boldsymbol{A},\qquad\mbox{where}\quad\mathrm{rank}(\boldsymbol{L})=K\ll p,\quad\mbox{and $\boldsymbol{A}$ is a ``nice" matrix}.$$

$$\mbox{Each diagonal of $\boldsymbol{A}$ is large compared with other entries in the same row}.$$

$$\mathcal{SDD}_{c}^{+}=\Bigl{\{}\boldsymbol{A}=(a_{ij})_{p\times p}:\boldsymbol{A}^{T}=\boldsymbol{A},\;a_{jj}\geq c\sum_{i:i\not=j}|a_{ji}|\ \mbox{ for\ all\ }1\leq j\leq p\Bigr{\}}.$$

$$\underset{(\mathbf{L},\mathbf{A})}{\min }\Vert \mathbf{S}-\mathbf{L}-\mathbf{A}{\Vert }_F,\text{subject to}\mathrm{rank}(\mathbf{L})\le K,\mathbf{L}={\mathbf{L}}^T,\mathbf{A}\in {\mathcal{S}\mathcal{DD}}_c^{+},$$

$${\mathbf{\Sigma }}^{-1}={\mathbf{A}}^{-1}-{\mathbf{A}}^{-1}\mathbf{B}({\mathbf{I}}_K+{\mathbf{B}}^T{\mathbf{A}}^{-1}\mathbf{B}{)}^{-1}{\mathbf{B}}^T{\mathbf{A}}^{-1}.$$

$$\Vert {\mathbf{A}}^{-1}\Vert \le \frac{c}{c-1}\Vert [\mathrm{diag}(\mathbf{A}){]}^{-1}\Vert ,\text{mbox}forany\mathbf{A}\in {\mathcal{S}\mathcal{DD}}_c^{+},\text{mbox}wherec>1.$$

$$X(j)=\sum _{k=1}^K{b}_k(j)W_k+Z(j),1\le j\le p,$$

$$\mathbf{A}\equiv \mathrm{Cov}(Z)\in {\mathcal{S}\mathcal{DD}}^{+}.$$

$$\boldsymbol{\Sigma}=\boldsymbol{B}\mathrm{Cov}(W)\boldsymbol{B}^{T}+\boldsymbol{A},\qquad\mbox{where}\quad\mathrm{rank}\bigl{(}\boldsymbol{B}\mathrm{Cov}(W)\boldsymbol{B}^{T}\bigr{)}=K\quad\mbox{and}\quad\boldsymbol{A}\in\mathcal{SDD}^{+}.$$

$$\widehat{\boldsymbol{\Sigma}}_{ddpca}=\widehat{\boldsymbol{L}}+\widehat{\boldsymbol{A}},\qquad\mbox{where $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ is the output of one-step DD-PCA.}$$

$${\hat{\mathbf{\Sigma }}}_{poet}={\hat{\mathbf{L}}}_{\ast }+{\mathcal{T}}_a({\hat{\mathbf{A}}}_{\ast }),\text{mbox}where{\hat{\mathbf{L}}}_{\ast }=\sum _{k=1}^K{\lambda }_k{\xi }_k{\xi }_k^T,{\hat{\mathbf{A}}}_{\ast }=\sum _{k=K+1}^p{\lambda }_k{\xi }_k{\xi }_k^T.$$

$$\mathcal{T}_{a}(\widehat{\boldsymbol{A}}_{*})=\widehat{\boldsymbol{D}}^{\frac{1}{2}}H_{a}\Bigl{(}\widehat{\boldsymbol{D}}^{-\frac{1}{2}}\widehat{\boldsymbol{A}}_{*}\widehat{\boldsymbol{D}}^{\frac{1}{2}}\Bigr{)}\widehat{\boldsymbol{D}}^{\frac{1}{2}}\qquad\text{where}\quad\widehat{\boldsymbol{D}}=\mathrm{diag}(\widehat{\boldsymbol{A}}_{*}),$$

$${\mathbf{\Sigma }}^{-1}={\mathbf{A}}^{-1}-{\mathbf{A}}^{-1}\mathbf{B}({\mathbf{I}}_K+{\mathbf{B}}^T{\mathbf{A}}^{-1}\mathbf{B}{)}^{-1}{\mathbf{B}}^T{\mathbf{A}}^{-1}.$$

$${\hat{\mathbf{\Sigma }}}^{-1}={\hat{\mathbf{A}}}^{-1}-{\hat{\mathbf{A}}}^{-1}\hat{\mathbf{B}}({\mathbf{I}}_K+{\hat{\mathbf{B}}}^T{\hat{\mathbf{A}}}^{-1}\hat{\mathbf{B}}{)}^{-1}{\hat{\mathbf{B}}}^T{\hat{\mathbf{A}}}^{-1}.$$

$$\widehat{\boldsymbol{\Sigma}}^{-1}=\widehat{\boldsymbol{A}}^{-1}-\widehat{\boldsymbol{A}}^{-\frac{1}{2}}\biggl{(}\sum_{k=1}^{K}\frac{1}{\hat{\sigma}_{k}^{-2}+1}\hat{\eta}_{k}\hat{\eta}_{k}^{\prime}\biggr{)}\widehat{\boldsymbol{A}}^{-\frac{1}{2}}$$

$$\Vert {\hat{\mathbf{A}}}^{-1}-{\mathbf{A}}^{-1}\Vert \le \Vert {\hat{\mathbf{A}}}^{-1}\Vert \Vert {\mathbf{A}}^{-1}\Vert \Vert \hat{\mathbf{A}}-\mathbf{A}\Vert .$$

$$\lambda_{\min}(\widehat{\boldsymbol{A}})\geq\min_{1\leq j\leq p}\Bigl{\{}\widehat{a}_{jj}-\sum_{i:i\neq j}|\widehat{a}_{ji}|\Bigr{\}}\geq\min_{1\leq j\leq p}\Bigl{\{}\widehat{a}_{jj}-c^{-1}\widehat{a}_{jj}\Bigr{\}}\geq\frac{c-1}{c}\min_{1\leq j\leq p}\widehat{a}_{jj}.$$

$$\Vert {\hat{\mathbf{A}}}^{-1}\Vert \le \frac{c}{c-1}\Vert [\mathrm{diag}(\hat{\mathbf{A}}){]}^{-1}\Vert .$$

$$\boldsymbol{w}_{*}=\mathrm{argmin}_{\boldsymbol{w}^{T}\boldsymbol{1}=1}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w},\qquad\mbox{$\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$: asset covariance matrix}.$$

$${\hat{\mathbf{w}}}_{\ast }={\mathrm{argmin}}_{{\mathbf{w}}^{T}1=1}{\mathbf{w}}^T\hat{\mathbf{\Sigma }}\mathbf{w}.$$

$${\hat{\mathbf{w}}}_{\ast }=(1^T{\hat{\mathbf{\Sigma }}}^{-1}1{)}^{-1}({\hat{\mathbf{\Sigma }}}^{-1}1).$$

$$R({\hat{\mathbf{w}}}_{\ast })=\frac{1}{T}\sum _{t=1}^T({\mathbf{y}}_t^T{\hat{\mathbf{w}}}^{\ast }{)}^2,$$

$$X\sim {\mathcal{N}}_p(\mu ,\mathbf{\Sigma }),$$

$$H_{0}:\mu=0,\qquad v.s.\qquad H_{1}:\mbox{$\mu\neq 0$, and $\mu$ is sparse}.$$

$$X=\mu+\sum_{k=1}^{K}w_{k}\eta_{k}+z,\quad w_{k}\sim\mathcal{N}(0,\nu_{k}),\;z\sim\mathcal{N}_{p}(0,\boldsymbol{A}),\;\mbox{$w_{1},\ldots,w_{K},z$ are independent}.$$

$$H{C}_p^{\ast }=\underset{1\le j\le p/2}{\max }H{C}_{p,j},\text{mbox}whereH{C}_{p,j}=\frac{\sqrt{p}[(j/p)-{\pi }_{(j)}]}{\sqrt{{\pi }_{(j)}(1-{\pi }_{(j)})}},1\le j\le p.$$

$$X-\sum _{k=1}^K{w}_k{\eta }_k\sim {\mathcal{N}}_p(\mu ,\mathbf{A}).$$

$$X=\mu +\sum _{k=1}^K{w}_k{\hat{\eta }}_k+z.$$

$${\mathcal{DD}}_c^{+}=\{\mathbf{A}=(a_{ij}{)}_{p\times p}:a_{jj}\ge c\sum _{i:i\ne j}|a_{ji}|\text{mbox}forallj\}$$

$$\underset{(\mathbf{L},\mathbf{A})}{\min }\frac{1}{2}\Vert \mathbf{S}-\mathbf{L}-\mathbf{A}{\Vert }_F^2+\lambda \Vert \mathbf{L}{\Vert }_{\ast }\text{subject to}\mathbf{A}\in {\mathcal{S}\mathcal{DD}}_c^{+}.$$

Full text

Keywords: ADMM, Alternating Projection, Approximate Factor Model, Covariance Matrix Estimation, De-Correlation, Higher Criticism, Global Testing, POET

1 Introduction

The approximate low-rankness is a popular structural assumption on covariance matrices. It assumes that a $p\times p$ covariance matrix $\boldsymbol{\Sigma}$ decomposes into

[TABLE]

Equivalently, it introduces a latent factor model on any random vector $X$ whose covariance matrix is $\boldsymbol{\Sigma}$, where $\boldsymbol{A}$ is the “residual covariance matrix” after the effects of latent variables are removed. Such a decomposition is not unique and varies with the meaning of a “nice” $\boldsymbol{A}$. One can impose different requirements on $\boldsymbol{A}$ to facilitate different applications. In the classical factor models for econometrics and finance, $\boldsymbol{A}$ is assumed a diagonal matrix (Fama and French, 1993) or a sparse matrix (Chamberlain and Rothschild, 1983), to enforce that the idiosyncratic noise accounts for little cross-sectional risk. In large-scale multiple testing, it is often assumed that the covariance matrix of test statistics has the above decomposition with $\boldsymbol{A}$ being a diagonal matrix (Leek and Storey, 2008) or having a small Frobenius norm (Fan et al., 2012). The motivation there is development of factor-adjusted multiple testing procedures, to make it legitimate to use conventional multiple testing methods on the post-factor-removal data. In image processing, a similar decomposition on image matrices was proposed (Candès et al., 2011), where $\boldsymbol{A}$ is assumed sparse, for the purpose of capturing details of images. In this paper, we explore a new type of approximate low-rankness where

[TABLE]

Translated to the latent variable representation, it means, after the effects of latent variables are removed, the correlation matrix of “residual” variables have uniformly small off-diagonal entries. One motivation of imposing this condition is to take into account the varying scale of the diagonal elements of $\boldsymbol{A}$. Most aforementioned approximate low-rank decompositions first perform PCA on $\boldsymbol{\Sigma}$ (or an empirical version of it) to remove the first a few principal components, and then conduct operations on the remaining matrix. It is often observed that the diagonal elements of the remaining matrix has considerable variations in magnitude. To deal with it requires careful adjustment on the post-PCA operations, such as adaptive thresholding (Cai and Liu, 2011). On the contrary, we impose the requirement (2) directly in the decomposition (1), in hopes of improving the PCA factors and easing the post-PCA operations. Another motivation of adopting the assumption (2) is to guarantee that $\boldsymbol{A}^{-1}$ is well-behaved. In many applications such as portfolio management and linear discriminant analysis, $\boldsymbol{A}^{-1}$ plays a key role (see Section 2). In the decomposition (1), forcing $\boldsymbol{A}$ to be a strictly diagonal matrix can ensure both $\boldsymbol{A}$ and $\boldsymbol{A}^{-1}$ are well-behaved, but this requirement is often too restrictive, and (2) is a natural relaxation. We note that imposing the common sparsity assumption on $\boldsymbol{A}$ does not even guarantee positive definiteness. Despite of remedies such as increasing the threshold or projection to the positive semi-definite cone (Fan et al., 2016), these approaches still don’t guarantee that $\boldsymbol{A}^{-1}$ is a “nice” matrix.

To formulate (2) mathematically, we define the set of “symmetric $c$-diagonally-dominant” matrices, for any $c>0$:

[TABLE]

For $c=1$, it reduces to the usual definition of diagonally-dominant matrices, and we omit the subscript and write $\mathcal{SDD}_{1}^{+}=\mathcal{SDD}^{+}$. Given a $p\times p$ positive semi-definite matrix $\boldsymbol{S}$, we introduce an optimization problem:

[TABLE]

where $\|\cdot\|_{F}$ is the matrix Frobenius norm. We call it Diagonally-Dominant Principal Component Analysis (DD-PCA). In this paper, we are primarily interested in $c=1$; discussions of $c\neq 1$ are deferred to Section 4.

The definition of DD-PCA is a nonconvex optimization with a rank constraint. Similar to solving other rank constrained optimizations in matrix completion, one can either solve a convex relaxation of (4) or develop an iterative algorithm that converges to a local minimum of (4). These ideas generate several variants of DD-PCA, as detailed in Section 4. Among those variants, one is of particular interest, which we call One-step DD-PCA:

• •

PCA: Obtain the $K$ leading eigenvalues and eigenvectors of $\boldsymbol{S}$, denoted as $\lambda_{1}\geq\ldots\geq\lambda_{K}\geq 0$ and $\xi_{1},\ldots,\xi_{K}\in\mathbb{R}^{p}$. Let $\boldsymbol{L}=\sum_{k=1}^{K}\lambda_{k}\xi_{k}\xi_{k}^{T}$.

• •

Projection to $\mathcal{SDD}^{+}$: Initialize $\boldsymbol{A}^{(0)}=\boldsymbol{S}-\boldsymbol{L}$ and $\boldsymbol{J}^{(0)}=\boldsymbol{0}$. For $t=1,2,\ldots,$

• –

Run the MRT algorithm (Mendoza et al., 1998)111The MRT algorithm computes the unique projection of a $p\times p$ matrix to the convex polyhedral cone consisting of all diagonally dominant matrices. It has a complexity of $O(p^{2}\log(p))$. See Section 4. to project $[\boldsymbol{A}^{(t-1)}-\boldsymbol{J}^{(t-1)}]$ into the diagonally-dominant cone. Let $\boldsymbol{G}^{(t)}$ be the projected matrix.

• –

Update $\boldsymbol{A}^{(t)}=\frac{1}{2}[\boldsymbol{G}^{(t-1)}+(\boldsymbol{G}^{(t-1)})^{T}]$ and $\boldsymbol{J}^{(t)}=\boldsymbol{J}^{(t-1)}+(\boldsymbol{G}^{(t)}-\boldsymbol{A}^{(t-1)})$.

• –

If $\|\boldsymbol{J}^{(t)}-\boldsymbol{J}^{(t-1)}\|_{F}\leq\epsilon$, stop and output $\boldsymbol{A}=\boldsymbol{A}^{(t)}$.

This method is obtained by running one outer-loop iteration in the iterative algorithm to be introduced in Section 4, explaining the name of one-step DD-PCA. It has the same philosophy as the one-step Huber estimator (Bickel, 1975) and one-step LLA implementation of non-convex penalized linear regressions (Zou and Li, 2008; Fan et al., 2014). It provides an approximate solution to (4), which is much faster to compute than solving (4) exactly.

Exploring the approximate low-rank structures is a powerful strategy for big data analysis. The classical PCA has motivated many statistical methods. Similarly, DD-PCA and one-step DD-PCA can also serve as building blocks for statistical methodology development. We exemplify it in two statistical problems: the first is estimating a large covariance matrix, and the second is testing of the global null hypothesis in multiple testing.

Estimation of large covariance matrices is a popular topic in statistical literatures (Fan et al., 2016). At the heart of it is two fundamental questions: (a) What structural assumption is appropriate? (b) How to evaluate the methods in real applications?

We adopt the structural assumption that the true covariance matrix $\boldsymbol{\Sigma}$ has an approximate low-rank decomposition with $\boldsymbol{A}\in\mathcal{SDD}_{c}^{+}$. This is a special type of factor covariance structures that are commonly used in econometrics (Fan and Fan, 2008), finance (Fama and French, 1993), genetics (Price et al., 2006) and many other fields. Our work is unique in the diagonal dominance assumption on $\boldsymbol{A}$. Intuitively, it is a natural relaxation of assuming $\boldsymbol{A}$ is diagonal, and it implies that, after the effects of latent factors are removed, the “residual variables” are almost uncorrelated. Compared with existing covariance matrix estimators that assume $\boldsymbol{A}$ is sparse (e.g., Fan et al. (2011, 2013a)), this diagonal dominance structure benefits simultaneous estimation of $\boldsymbol{\Sigma}$ and $\boldsymbol{\Sigma}^{-1}$: In factor covariance structures, the singular values of the low-rank matrix are much larger than $\|\boldsymbol{A}\|$, so the error of estimating $\boldsymbol{\Sigma}$ is dominated by the error of recovering the low-rank part. If our goal is merely to estimate $\boldsymbol{\Sigma}$, we do not gain much from exploring the diagonal dominance structure of $\boldsymbol{A}$. However, if we are also interested in estimating $\boldsymbol{\Sigma}^{-1}$, the error of estimating $\boldsymbol{A}^{-1}$ will play a key role. Note that there always exists a matrix $\boldsymbol{B}\in\mathbb{R}^{n\times K}$ such that $\boldsymbol{L}=\boldsymbol{B}\boldsymbol{B}^{T}$. It follows from the matrix inverse formula (Horn and Johnson, 2012) that

[TABLE]

Suppose we have obtained a good estimator $\widehat{\boldsymbol{\Sigma}}=\widehat{\boldsymbol{B}}\widehat{\boldsymbol{B}}^{T}+\widehat{\boldsymbol{A}}$ by fitting some factor covariance structure on the data. Even though $\|\widehat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}\|$ is small, it is still possible that $\|\widehat{\boldsymbol{A}}^{-1}-\boldsymbol{A}^{-1}\|$ is large so that $\widehat{\boldsymbol{\Sigma}}^{-1}$ is far from being a good estimator of $\boldsymbol{\Sigma}^{-1}$. Fortunately, exploring the diagonal dominance structure largely mitigates this issue, thanks to an appealing feature of the diagonally-dominant cone $\mathcal{SDD}_{c}^{+}$ (Horn and Johnson, 2012):

[TABLE]

Therefore, if we enforce $\widehat{\boldsymbol{A}}\in\mathcal{SDD}_{c}^{+}$ in fitting the factor covariance structure, for a constant $c>1$, then $\|\widehat{\boldsymbol{A}}^{-1}\|$ won’t explode, preventing ill behavior of $\widehat{\boldsymbol{\Sigma}}^{-1}$. To this end, we propose a new covariance matrix estimator $\widehat{\boldsymbol{\Sigma}}_{ddpca}$ using the solution of DD-PCA or one-step DD-PCA. We demonstrate in numerical studies: $\widehat{\boldsymbol{\Sigma}}_{ddpca}$ has comparable performance with state-of-art methods, but the new estimator is tuning free once $K$ is specified, so is more convenient to use. Moreover, when it comes to estimating $\boldsymbol{\Sigma}^{-1}$ by $\widehat{\boldsymbol{\Sigma}}^{-1}_{ddpca}$, the new estimator is significantly better than inverting other factor-based covariance matrix estimators.

In real applications, estimating the covariance matrix is rarely the ultimate goal. Often, it serves as an intermediate step for downstream tasks. We demonstrate the usefulness of our covariance estimator by evaluating its performance in two downstream tasks, portfolio management and linear discriminant analysis. In the former, an estimate of the covariance matrix is needed to obtain Markowitz’s optimal portfolio weights; in the latter, it is used to compute Fisher’s LDA classifier. Note that what is actually plugged into these downstream tasks is the inverse of estimated covariance matrix. As we have argued, the main advantage of our method is on estimating $\boldsymbol{\Sigma}^{-1}$ by $\widehat{\boldsymbol{\Sigma}}_{ddpca}^{-1}$, a perfect match to these applications. This is supported by encouraging real data results. It is worthwhile mentioning that our approach is different from the approach of plugging in an existing precision matrix estimator (e.g., the graphical lasso (Friedman et al., 2008)). These methods assume $\boldsymbol{\Sigma}^{-1}$ is sparse, while we assume a factor structure on $\boldsymbol{\Sigma}$. For portfolio data, adopting a factor covariance structure is the common practice. For classification problems, there are also many real data sets for which it is appropriate to assume a factor structure (see Section 2).

The second application of DD-PCA is for testing of the global null hypothesis in multiple testing. We are primarily interested in the setting where the mean vector in the alternative hypothesis is sparse. The Higher Criticism (HC) test (Donoho and Jin, 2004) is known to enjoy theoretical optimality and has gained increasing popularity in applications (Donoho and Jin, 2015; Wu et al., 2011). However, the orthodox HC test assumes that the individual $z$-scores are mutually independent. There is limited understanding of how to extend the HC test to the case where $z$-scores share common latent factors. We model that the covariance matrix of $z$-scores, $\boldsymbol{\Sigma}$, has a low rank plus diagonal dominance decomposition. We propose two variants of HC for this setting, both utilizing DD-PCA as a module. The first test is a modification of the Innovated HC test (Hall and Jin, 2010) by plugging in the estimator of $\boldsymbol{\Sigma}^{-1}$ from DD-PCA. The second test uses the low-rank matrix $\widehat{\boldsymbol{L}}$ from DD-PCA to remove the effects of latent factors and then applies the orthodox HC test. Both tests significantly improve over the orthodox HC test and the Innovated HC test in simulations. The rationale of these testing ideas is to “transform” and “decorrelate” the marginal $z$-scores. The first test extends the Innovated Transformation of $z$-scores (Hall and Jin, 2010), $X\mapsto\boldsymbol{\Sigma}^{-1}X$, from the case where $\boldsymbol{\Sigma}^{-1}$ is sparse to the case where $\boldsymbol{\Sigma}$ has a factor structure. The second test aims to estimate and remove the latent factors from $z$-scores so that the ‘residuals’ are almost uncorrelated (i.e., their covariance matrix is close to being diagonal). A similar idea was briefly mentioned in Fan et al. (2012) (called factor-adjusted $z$-scores), but it has never been used in the global testing problem. In both new tests, we can replace DD-PCA by the classical PCA, but the numerical performance will deteriorate. This indicates that exploring the diagonal dominance structure is beneficial.

The remaining of this paper is organized as follows. In Section 2, we introduce a new covariance matrix estimator powered by DD-PCA, and discuss its applications in portfolio optimization (Section 2.1) and linear discriminant analysis (Section 2.2). In Section 3, we propose two new test statistics, using DD-PCA as a building block, for testing of the global null hypothesis in multiple testing. In Section 4, we address the computation of DD-PCA, by introducing an ADMM algorithm and an iterative projection algorithm for conducting the decomposition (1)-(2) for any given covariance matrix. Section 5 contains simulations, and Section 6 contains concluding remarks. The appendix contains some algorithm details omitted in Section 4.

2 Estimating large covariance matrices by DD-PCA

Let $X\in\mathbb{R}^{p}$ be a multivariate random vector with a covariance matrix $\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$, where $p$ is presumably much larger than $n$. We adopt a factor model:

[TABLE]

where $W_{1},\ldots,W_{K}$ are unobserved random variables (factors), $b_{k}\in\mathbb{R}^{p}$ is a nonrandom vector containing the loadings of the $k$-th factor, and $Z\in\mathbb{R}^{p}$ is a random vector independent of the factors such that

[TABLE]

Given $iid$ data $X_{1},\ldots,X_{n}\in\mathbb{R}^{p}$, we are interested in estimating $\boldsymbol{\Sigma}$ and $\boldsymbol{\Sigma}^{-1}$.

By model (5)-(6), the covariance matrix of $X$ has a decomposition

[TABLE]

It has the low-rank plus diagonal dominance structure. We propose the following estimator: Let $\boldsymbol{S}=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\bar{X})(X_{i}-\bar{X})^{T}$ be the sample covariance matrix. Take $\boldsymbol{S}$ as the input to the one-step DD-PCA algorithm in Section 1 and let $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ be the output. We estimate $\boldsymbol{\Sigma}$ by

[TABLE]

We then estimate $\boldsymbol{\Sigma}^{-1}$ by the inverse of $\widehat{\boldsymbol{\Sigma}}_{ddpca}$. Here, $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ can be replaced by the output of other variants of DD-PCA (see Section 4). They give similar numerical performance, so we stick to one-step DD-PCA for computational convenience.

Different from existing covariance estimation methods under factor structures, our approach imposes the diagonal dominance constraint on $\boldsymbol{A}$. We now compare it with methods that impose the sparsity constraint on $\boldsymbol{A}$. One popular method is POET (Fan et al., 2013a). Let $\boldsymbol{S}=\sum_{k=1}^{p}\lambda_{k}\xi_{k}\xi_{k}^{T}$ be the eigen-decomposition of $\boldsymbol{S}$, where $\lambda_{k}$ and $\xi_{k}$ are the $k$-th eigenvalue and eigenvector, respectively. The POET estimator is

[TABLE]

Here, $\mathcal{T}_{a}(\cdot)$ can be any entry-wise adaptive thresholding operator (Rothman et al., 2009; Cai and Liu, 2011; Xue et al., 2012). Fan et al. (2013a) suggest using the hard-thresholding operator applied to a “correlation matrix” associated with $\widehat{\boldsymbol{A}}_{*}$, i.e.,

[TABLE]

and $H_{a}(\cdot)$ is the entry-wise hard-thresholding at the threshold $a>0$. Then, an estimate of $\boldsymbol{\Sigma}^{-1}$ is obtained by $\widehat{\boldsymbol{\Sigma}}^{-1}_{poet}$.

Figure 1 gives a simulation example. Fix $(p,n,K)=(2000,200,3)$. We generate data from the model (5), where the factors $\{W_{k}(i):1\leq k\leq K,1\leq i\leq n\}$ are drawn $iid$ from $N(0,1)$, the factor loadings $\{b_{k}(j):1\leq k\leq K,1\leq j\leq p\}$ are generated $iid$ from ${\cal N}(0,1)$, and the noise vectors $Z_{1},\ldots,Z_{n}$ are drawn $iid$ from a multivariate normal ${\cal N}_{p}(\boldsymbol{0},\boldsymbol{A})$, with $A(i,j)=0.5^{|i-j|+1}$ for $i\neq j$ and $1$ otherwise. For both methods, $K$ is unknown and treated as a tuning integer. POET has an additional tuning threshold $a$, which is selected by cross-validation (default procedure in the poet package).222This default procedure guarantees that $\widehat{\boldsymbol{\Sigma}}_{poet}$ is invertible.

In the top four panels of Figure 1, we show the average estimation errors on $\boldsymbol{\Sigma}$ and $\boldsymbol{\Sigma}^{-1}$ over $100$ repetitions. Since $K$ is unknown, we implement both methods for the true $K=3$ and misspecified $K\in\{4,5,\ldots,8\}$.333We don’t include the results of $K\in\{1,2\}$, as the errors are much larger. For estimating $\boldsymbol{\Sigma}$, the two methods give very similar performance. This is not surprising. Since the eigenvalues of the low-rank part are much larger than $\|\boldsymbol{A}\|$, the error of estimating $\boldsymbol{\Sigma}$ is dominated by the error of recovering the low-rank part. Our method and POET has the same low-rank part (the $\widehat{\boldsymbol{L}}$ from one-step DD-PCA and the $\widehat{\boldsymbol{L}}_{*}$ in (8) are indeed the same), so they have similar errors on estimating $\boldsymbol{\Sigma}$. From the bottom left two panels of Figure 1, we can see that our method does a better job on estimating $\boldsymbol{A}$; especially, the spectral norm error is $10$-$20$% smaller. However, this improvement is almost negligible compared with the errors on recovering the low-rank part. We conclude that our method and POET have similar performance on estimating $\boldsymbol{\Sigma}$. Still, our method has an advantage: It has no tuning threshold and is more convenient to use.

How about the performance on estimating $\boldsymbol{\Sigma}^{-1}$? The top right two panels of Figure 1 clearly suggest that our method has a significant advantage. When $K=3$, the spectral norm error of our method is only one half of the error of POET. Interestingly, the performance of POET improves with an overshooting $K$; but even for $K=8$, its spectral norm error is still $10\%$ larger than the error of our method. For the Frobenius norm error, our estimator also outperforms POET for all choices of $K$. This phenomenon is due to that $\widehat{\boldsymbol{A}}$ plays a dominating role when we compute the inverse of $\widehat{\boldsymbol{\Sigma}}$, while the low-rank part has a negligible effect, so the advantage of our method on recovering $\boldsymbol{A}$ becomes prominent. This is illustrated in the bottom right two panels of Figure 1. Recall that $\widehat{\boldsymbol{A}}$ is from one-step DD-PCA and $\widehat{\boldsymbol{A}}_{*}$ is as in (8). The Frobenius/spectral norm of $(\widehat{\boldsymbol{A}}^{-1}-\boldsymbol{A}^{-1})$ is significantly smaller than the Frobenius/spectral norm of $(\widehat{\boldsymbol{A}}_{*}^{-1}-\boldsymbol{A}^{-1})$. Additionally, by comparing the top right two panels with the bottom right two panels, we can see that the error of estimating $\boldsymbol{\Sigma}^{-1}$ is almost determined by the error of estimating $\boldsymbol{A}^{-1}$.

This numerical example delivers two messages: First, compared with competitive factor-based methods, the major advantage of our method is on estimating $\boldsymbol{\Sigma}^{-1}$ by $\widehat{\boldsymbol{\Sigma}}^{-1}$. Second, such an advantage is driven by the better accuracy on recovering $\boldsymbol{A}^{-1}$. Below, we explain them using linear algebra.

Without loss of generality, in model (5), we assume the covariance matrix of $W$ equals to the identify matrix. Then, $\boldsymbol{\Sigma}=\boldsymbol{B}\boldsymbol{B}^{T}+\boldsymbol{A}$. By matrix inverse formula,

[TABLE]

Suppose we construct an estimator $\widehat{\boldsymbol{\Sigma}}=\widehat{\boldsymbol{B}}\widehat{\boldsymbol{B}}^{T}+\widehat{\boldsymbol{A}}$ from fitting a factor-type covariance structure. Then, $\widehat{\boldsymbol{\Sigma}}^{-1}$ (if it exists) has a similar decomposition:

[TABLE]

By some basic linear algebra, we can derive the following proposition:

Proposition 1**.**

Let $\widehat{\boldsymbol{A}}^{-\frac{1}{2}}\widehat{\boldsymbol{B}}=\sum_{k=1}^{K}\hat{\sigma}_{k}\hat{\eta}_{k}\hat{h}_{k}^{\prime}$ be the singular value decomposition of $\widehat{\boldsymbol{A}}^{-\frac{1}{2}}\widehat{\boldsymbol{B}}$, where $\hat{\sigma}_{k}>0$ is the $k$-th singular value and $\hat{\eta}_{k}\in\mathbb{R}^{p}$ and $\hat{h}_{k}\in\mathbb{R}^{K}$ are the corresponding left and right singular vectors. Then,

[TABLE]

By (10), the error of recovering the low-rank part only affects the matrix in the brackets. For $(\widehat{\boldsymbol{A}},\widehat{\boldsymbol{B}})$ obtained in factor-based methods, nonzero eigenvalues of $\widehat{\boldsymbol{B}}\widehat{\boldsymbol{B}}^{T}$ are much larger than $\|\widehat{\boldsymbol{A}}\|$, so $\widehat{\sigma}_{k}$’s are all very large. Then, the matrix in the brackets can hardly bring in a large error in $\widehat{\boldsymbol{\Sigma}}^{-1}$. The error in $\widehat{\boldsymbol{\Sigma}}^{-1}$ mainly comes from the error in $\widehat{\boldsymbol{A}}^{-1}$.

We further investigate the error in $\widehat{\boldsymbol{A}}^{-1}$. Note that

[TABLE]

To achieve a small $\|\widehat{\boldsymbol{A}}-\boldsymbol{A}\|$ by imposing structural assumptions on $\boldsymbol{A}$ is not too difficult. However, it typically does not prevent $\|\widehat{\boldsymbol{A}}^{-1}\|$ from exploding. For example, if $\widehat{\boldsymbol{A}}$ is obtained from entry-wise thresholding, we need a comparably large threshold to control $\|\widehat{\boldsymbol{A}}^{-1}\|$, but unfortunately we cannot let the threshold be too large as it significantly increases $\|\widehat{\boldsymbol{A}}-\boldsymbol{A}\|$. It turns out that, if we restrict $\widehat{\boldsymbol{A}}\in\mathcal{SDD}_{c}^{+}$ for a constant $c>1$, then it is automatically guaranteed that $\|\widehat{\boldsymbol{A}}^{-1}\|$ has a nice bound. As a property of diagonally-dominant matrices (Horn and Johnson, 2012), for $c>1$,

[TABLE]

It follows that

[TABLE]

This explains why the constraint of $\widehat{\boldsymbol{A}}\in\mathcal{SDD}_{c}^{+}$ helps significantly reduce the errors in $\widehat{\boldsymbol{A}}^{-1}$ and (ultimately) the errors in $\widehat{\boldsymbol{\Sigma}}^{-1}$.

The above argument applies to $c>1$. In our method, $c=1$. Sometimes, we may even have to use $c<1$, so that the assumption $\boldsymbol{A}\in\mathcal{SDD}_{c}^{+}$ is not too restrictive (see Section 4). For $c\leq 1$, we do not have a solid argument as (12), but a similar phenomenon is observed in numerical studies.

Below, we use two real applications to further demonstrate that exploring the diagonal dominance factor structures is a useful strategy.

2.1 Application to portfolio management

Given a collection of $p$ assets, portfolio management aims to determine the weights allocated to each asset. It is often desirable to construct the minimum risk portfolio, where the asset weights $\boldsymbol{w}^{*}=(w_{1}^{*},\ldots,w_{p}^{*})$ are determined by

[TABLE]

In practice, $\boldsymbol{\Sigma}$ is unknown. We first obtain an estimate $\widehat{\boldsymbol{\Sigma}}$ using asset returns $\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{n}\in\mathbb{R}^{p}$ during a period of $n$ days, then we estimate the weights by

[TABLE]

This optimization has an explicit solution:

[TABLE]

Since what we actually need is $\widehat{\boldsymbol{\Sigma}}^{-1}$, exploring the low-rank plus diagonal dominance structure is a potentially useful strategy.

We compare our method with POET on real data. We collected the daily returns of stocks in S&P 100 index from January 1st 2006 to December 31st 2016. After removing companies that were listed after 2006, there are $80$ stocks in total. On the first trading day of each month, we created two portfolios from (13), where $\widehat{\boldsymbol{\Sigma}}$ is estimated using daily returns for the proceeding 12 months ($n=252$) by our method and by POET, respectively. We set $K=3$ for both methods. The threshold in POET is chose by cross-validation (we use the default cross-validation procedure in poet package). On the last trading day of the same month, we measure the actual risk of each portfolio by

[TABLE]

where $T$ is the number of trading days in this month ($T=21$ for most months) and $\boldsymbol{y}_{t}\in\mathbb{R}^{80}$ contains the stock returns on day $t$ of the month. Define $r=(R_{poet}-R_{ddpca})/R_{ddpca}$; note that a positive $r$ indicates that the portfolio created using our method is superior to that of POET. Figure 2 displays the histogram of $r$ over $120$ months in our data range. It suggests that our method improves POET by $9.5\%$ on average and $14.7\%$ in the median.

2.2 Application to linear discriminant analysis

In binary classification, given feature vectors $X_{1},\ldots,X_{n}\in\mathbb{R}^{p}$ and training labels $\ell_{1},\ldots,\ell_{n}\in\{1,2\}$, we aim to construct a linear classifier. In the classical regime where $p$ is fixed as the training sample size grows, Fisher’s LDA is an effective linear classifier. In the modern high dimensional settings where $p\gg n$, it has been well understood that feature screening is necessary before one applies Fisher’s LDA (Fan and Fan, 2008; Donoho and Jin, 2008), and that it is desirable to plug in a good estimate of the inverse covariance matrix that explores structural assumptions (Cai et al., 2011). Recently, Fan et al. (2013b) proposed a linear classifier that uses an estimate of inverse covariance matrix in both the screening step and LDA step, and they showed that this classifier is rate-optimal under a multivariate normal model with even extremely weak signal strength. This classifier was later applied to several large real classification problems with superior results (Huang et al., 2016). We shall combine our covariance matrix estimator with this classifier to see whether exploring the low-rank plus diagonal-dominance structure is helpful.

Given an estimate $\widehat{\boldsymbol{\Omega}}$ of the inverse covariance matrix and a threshold $t>0$, the classifier has four steps (Huang et al., 2016):

Calculate the feature-wise $t$-score: For $1\leq j\leq p$, let $Z(j)=[\bar{X}_{1}(j)-\bar{X}_{2}(j)]/(n\cdot s_{j})$, where $\bar{X}_{1}(j)$ and $\bar{X}_{2}(j)$ are the within-class sample means of feature $j$ and $s_{j}>0$ is the pooled standard deviation of feature $j$. Write $Z=(Z(1),\ldots,Z(p))^{T}$. 2. 2.

Apply the Innovated Transformation (Fan et al., 2013b) to get $\tilde{Z}=\widehat{\boldsymbol{\Omega}}Z$. 3. 3.

Feature-wise thresholding: For $1\leq j\leq p$, let $w(j)=\text{sgn}(\tilde{Z}(j))\cdot 1\{|\tilde{Z}(j)|\geq t\}$. Write $w=(w(1),w(2),\dots,w(p))^{T}$. 4. 4.

Classification by LDA. Given a test feature vector $\widetilde{X}\in\mathbb{R}^{p}$, for $1\leq j\leq p$, normalize $\widetilde{X}(j)$ to $\widetilde{X}^{\ast}(j)=[\tilde{X}(j)-\frac{1}{2}(\bar{X}_{1}(j)+\bar{X}_{2}(j))]/s_{j}$, where $(\bar{X}_{1}(j),\bar{X}_{2}(j),s_{j})$ are the same as in Step 1. Write $\widetilde{X}^{\ast}=(\widetilde{X}^{\ast}(1),\ldots,\widetilde{X}^{\ast}(p))^{T}$. We classify the test sample to class 1 if $w^{T}\widehat{\boldsymbol{\Omega}}\widetilde{X}^{\ast}>0$ and to class 2 otherwise.

In this classifier, the matrix $\widehat{\boldsymbol{\Omega}}$ plays two roles: First, it is used in the Innovated Transformation, so different $\widehat{\boldsymbol{\Omega}}$ leads to different feature rankings. Second, it is used in the LDA step, so $\widehat{\boldsymbol{\Omega}}$ also affects the classification boundary.

We compare the classification performance of plugging in three versions of $\widehat{\boldsymbol{\Omega}}$: The first is $\widehat{\boldsymbol{\Sigma}}_{ddpca}^{-1}$, the second is $\widehat{\boldsymbol{\Sigma}}^{-1}_{poet}$, and the last is $[\mathrm{diag}(\boldsymbol{S})]^{-1}$, where $\boldsymbol{S}$ is the sample covariance matrix. We note that the last approach is indeed the method FAIR (Fan and Fan, 2008). The above classifier also requires a threshold $t>0$. To minimize the effects of selecting $t$, for each $1\leq k\leq p$, we set the threshold such that $k$ features are retained and record the classification error. This generates an error curve for each method as $k$ ranges from $1$ to $p$.

We consider two datasets: the lung cancer dataset (Gordon et al., 2002) and the breast cancer dataset (Wang et al., 2005). They were downloaded from http://blog.nus.edu.sg/staww/softwarecode/. For both datasets, we conducted a pre-processing by ranking all features by the feature-wise $t$-score and retaining $p_{0}$ top-ranked features, where $p_{0}$ is a number that is for sure larger than the true number of useful features (but $p_{0}\ll p$).

[TABLE]

The lung cancer dataset was analyzed in various papers (Tibshirani et al., 2002; Fan and Fan, 2008). The estimated number of useful features by these methods is around $30$, so we confidently set $p_{0}=100$. The breast cancer dataset is a more difficult one and requires a lot more retained features. Jin and Wang (2016) analyzed the dataset under a clustering framework and suggested that the number of useful features is 728, so we set $p_{0}=1000$. We also tried other choices of $p_{0}$ (e.g., $p_{0}=200$ for lung cancer data and $p_{0}=2000$ for breast cancer data), and the results are similar.

We evaluate the classification performance by a 5-fold cross-validation procedure with stratified sampling. In detail, we randomly divide samples from class 1 into five folds and do the same to samples from class 2; we then re-combine them to five folds, such that the fraction of class 1 is the same across all folds. Next, we successively leave out each fold, train the classifier on remaining samples, and compute the test error on leave-out samples. The misclassification error reported is the average over 5 folds.

Figure 3 displays the results on lung cancer dataset. POET and DD-PCA have a tuning integer $K$, and we tried $K\in\{1,2,3\}$. The results suggest that, as long as more than $10$ features are retained, the classifier powered by DD-PCA uniformly outperforms the other two. Especially, for $K\in\{2,3\}$, the error keeps as low as $1/181$ once the number of retained features exceeds $60$. The performance of POET is slightly worse than FAIR for $K\in\{1,2\}$, and slightly better for $K=3$. We emphasize that the estimated inverse covariance matrix affects both the feature ranking and the LDA; therefore, even when the number of retained features is the same, the actual retained features are different across different methods.

Figure 4 displays the results on breast cancer dataset. For both POET and DD-PCA, $K\in\{4,5\}$ is favored to $K=3$. When $K=4$, as the number of retained features is in the interval of $[500,700]$, DD-PCA achieves the smallest error of $114/276$. In all three panels, the lowest attainable error of DD-PCA is smaller than those of POET and FAIR.

3 Detecting sparse mean effects by DD-PCA

The global detection is a problem of great interest in multiple testing (Simes, 1986; Donoho and Jin, 2004; Wu et al., 2011). Let $X_{1},\ldots,X_{p}\in\mathbb{R}$ be the $z$-scores of $p$ tests, where $p$ is presumably large. We assume

[TABLE]

where $\mu$ contains the true effects of these tests and $\boldsymbol{\Sigma}$ captures the dependence among the $z$-scores. We are interested in testing

[TABLE]

When $\boldsymbol{\Sigma}$ is a diagonal matrix, this problem has been studied extensively in the literature. Various tests were proposed, such as the $\chi^{2}$ test (or Hotelling’s $T^{2}$ test), maximum entry test (or minimum $p$-value test), Higher-Criticism test (Donoho and Jin, 2004), Berk-Jones test (Jager and Wellner, 2007), etc.. The Higher-Criticism (HC) test achieves the theoretically optimal detection boundary when the nonzero effects in $\mu$ are rare and weak (Donoho and Jin, 2004) and has gained increasing popularity in real applications (Donoho and Jin, 2015). However, there is limited understanding of how to use the HC test beyond a diagonal $\boldsymbol{\Sigma}$. When $\boldsymbol{\Sigma}$ is heavily non-sparse, a brute-forth application of the orthodox HC test leads to suboptimal performance (Hall and Jin, 2008). A satisfactory answer is only available when $\boldsymbol{\Sigma}$ is row-wise sparse. Hall and Jin (2010) introduced the “Innovated Transformation” on data, $X\mapsto\widehat{\boldsymbol{\Omega}}X$, where $\widehat{\boldsymbol{\Omega}}$ is an estimator of $\boldsymbol{\Sigma}^{-1}$. They showed that an application of the HC test on the post-transformation data lead to theoretically optimal testing performance.

Motivated by the popularity of adopting factor covariance structures in multiple testing (Leek and Storey, 2008; Fan et al., 2012), we consider the global testing problem (14)-(15) by assuming that $\boldsymbol{\Sigma}$ has a low rank plus diagonal dominance structure as in (1)-(2). Denote by $\sum_{k=1}^{K}\nu_{k}\eta_{k}\eta_{k}^{T}$ the eigen-decomposition of $\boldsymbol{L}$. Equivalently, we model that

[TABLE]

Here, $w_{1},\ldots,w_{K}$ are latent variables that account for most of the heavy dependence among test statistics. Under this model, $\boldsymbol{\Sigma}$ is heavily non-sparse, so the orthodox HC test performs unsatisfactorily (Hall and Jin, 2008). At the same time, $\boldsymbol{\Sigma}^{-1}$ may not be row-wise sparse, so the Innovated HC test (Hall and Jin, 2010) is not necessarily a good choice either. How to adapt the HC test to factor covariance structures is still largely unclear. We use DD-PCA (and its variants such as one-step DD-PCA) to develop two modifications of the HC test. Both tests significantly outperform the orthodox HC and Innovated HC, when the factor covariance structure holds.

Without loss of generality, we assume an estimate of $\boldsymbol{\Sigma}$ is available, denoted as $\widehat{\boldsymbol{\Sigma}}$. Note that it is common that a $z$-score $X_{j}$ is computed from a number of repeated observations. Suppose we observe $iid$ samples $X_{1},X_{2},\ldots,X_{n}$ and obtain the $z$-scores as $X=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}X_{i}$, then the sample covariance matrix of $X_{i}$’s can be used as $\widehat{\boldsymbol{\Sigma}}$.

We begin with describing the orthodox HC. Let $\pi_{1},\pi_{2},\ldots,\pi_{p}$ be the marginal $P$-values computed from $X_{j}\sim{\cal N}(0,\widehat{\Sigma}_{jj})$, where $\widehat{\Sigma}_{jj}$ is the $j$-th diagonal of $\widehat{\boldsymbol{\Sigma}}$. Sort the $P$-values in the descending order and denote by $\pi_{(k)}$ the $k$-th smallest $P$-value, $1\leq k\leq p$. The OHC test statistic is

[TABLE]

The null distribution of $HC_{p}^{*}$ is often approximated by a Gumbel distribution or by simulating the null data (Donoho and Jin, 2015). Although OHC was proposed for the case of a diagonal $\boldsymbol{\Sigma}$, we can treat it as a blackbox procedure (we actually know what is happening in the ‘blackbox’, but we have no intention to interfere): It takes as input the $P$-value for each individual test, and outputs a test statistic for the global null hypothesis which is an aggregation of all individual $P$-values. Now, if we feed this ‘blackbox’ with a different set of individual $P$-values, it will output a different test statistic. Following this strategy, we modify OHC by constructing individual $P$-values using $\widehat{\boldsymbol{\Sigma}}$, in hopes of borrowing strength from each other.

The first test we propose is IHC-DD test, where ‘IHC’ stands for Innovated HC and ‘DD’ stands for DD-PCA. As the name has suggested, this test is built on top of the Innovated HC test—we replace $\widehat{\boldsymbol{\Omega}}$ in IHC by the DD-PCA estimator of $\boldsymbol{\Sigma}^{-1}$ introduced in Section 2.

• •

Take $\widehat{\boldsymbol{\Sigma}}$ as the input to the one-step DD-PCA algorithm in Section 1 and let $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ be the output. Let $\widehat{\boldsymbol{\Sigma}}_{ddpca}=\widehat{\boldsymbol{L}}+\widehat{\boldsymbol{A}}$.

• •

Obtain $\widetilde{X}=\widehat{\boldsymbol{\Sigma}}_{ddpca}^{-1}X$. Compute the indivisual $P$-values $\widetilde{\pi}_{j}$ from the null distribution of $\widetilde{X}_{j}\sim{\cal N}(0,\widehat{\Omega}_{jj})$, where $\widehat{\Omega}_{jj}$ is the $j$-th diagonal of $\widehat{\boldsymbol{\Sigma}}_{ddpca}^{-1}$.

• •

Input the $P$-values $\widetilde{\pi}_{1},\ldots,\widetilde{\pi}_{p}$ to the OHC procedure (17) to get a test statistic.

The individual $P$-values fed to OHC are different from before: Each $\widetilde{\pi}_{j}$ borrows information from $z$-scores of other tests, taking advantage of the dependence between tests. The IHC-DD test can be viewed as another application of the DD-PCA covariance estimator proposed in Section 2, where we simply plug the estimated $\boldsymbol{\Sigma}^{-1}$ into the existing IHC test.

The second test we propose is more customized to the factor covariance structure. We call it DD-HC test, to differentiate it from the test above. By (16),

[TABLE]

Provided with estimates of $w_{k}$’s, $\eta_{k}$’s, and $\boldsymbol{A}$, we can use $X^{*}=X-\sum_{k=1}^{K}\widehat{w}_{k}\widehat{\eta}_{k}$ as the new $z$-scores for all individual tests, and we can compute the individual $P$-values from the null distribution of $X_{j}^{*}\sim{\cal N}(0,\widehat{A}_{jj})$, where $\widehat{A}_{jj}$ is the $j$-th diagonal of $\widehat{\boldsymbol{A}}$. Let’s discuss how to estimate $(w_{1},\ldots,w_{K},\eta_{1},\ldots,\eta_{K},\boldsymbol{A})$ from DD-PCA. Let $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ be the output of DD-PCA with $\widehat{\boldsymbol{\Sigma}}$ as the input (we will discuss below which DD-PCA algorithm to use). Recall that in model (16), $\eta_{k}$’s are eigenvectors of $\boldsymbol{L}$. This motivates us to estimate $\eta_{k}$’s by the leading eigenvectors of $\widehat{\boldsymbol{L}}$. Once we have $\widehat{\eta}_{1},\ldots,\widehat{\eta}_{K}$, we can approximate model (16) by

[TABLE]

This is indeed a linear regression model with $K$ covariates and $p$ observations. The regression coefficients are $w_{1},\ldots,w_{K}$, and each observation has an individual intercept $\mu_{j}$. Since $\mu$ is a sparse mean vector, the majority of $\mu_{j}$’s are zero; we thus estimate $w_{1},\ldots,w_{K}$ by a robust regression via minimizing $\|X-\sum_{k=1}^{K}w_{k}\widehat{\eta}_{k}\|_{1}$, with respect to $w_{k}$’s. Last, we discuss the estimation of $\boldsymbol{A}$. A straightforward idea is to use $\widehat{\boldsymbol{A}}$ output by DD-PCA. However, we don’t recommend this approach. The estimator of $\boldsymbol{A}$ shall be used to approximate the null distribution of $X_{j}^{*}=X_{j}-\sum_{k=1}^{K}\widehat{w}_{k}\widehat{\eta}_{k}(j)$, in order to compute individual $P$-values. So, our goal is not to estimate $\boldsymbol{A}$ accurately but to provide an accurate estimate of the variance of $X_{j}^{*}$. We have to take into account the plug-in effect of $\widehat{w}_{k}$ and $\widehat{\eta}_{k}$. Since $X_{j}^{*}$ is the residual of fitting a linear regression, a more reasonable choice is to use the diagonals of $\widehat{\boldsymbol{\Sigma}}-\widehat{\boldsymbol{L}}$ as estimates of the variances of $X_{j}^{*}$’s. We summarize the procedure as follows:

• •

Take $\widehat{\boldsymbol{\Sigma}}$ as the input to a DD-PCA algorithm and and let $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ be the output.

• •

Conduct PCA on $\widehat{\boldsymbol{L}}$, and denote by $\widehat{\eta}_{1},\ldots,\widehat{\eta}_{K}$ the first $K$ eigenvectors of $\widehat{\boldsymbol{L}}$.

• •

Regress $X$ on $(\widehat{\eta}_{1},\ldots,\widehat{\eta}_{K})$ using a robust regression: $\min_{w_{1},\ldots,w_{K}}\|X-\sum_{k=1}^{K}w_{k}\widehat{\eta}_{k}\|_{1}$. Let $(\widehat{w}_{1},\ldots,\widehat{w}_{K})$ be the estimated coefficients.

• •

Obtain $X^{*}=X-\sum_{k=1}^{K}\widehat{w}_{k}\widehat{\eta}_{k}$. Compute the individual $P$-values $\pi^{*}_{j}$ from the null distribution of $X_{j}^{*}\sim{\cal N}(0,\widehat{R}_{jj})$, where $\widehat{R}_{jj}$ is the $j$-th diagonal of $\widehat{\boldsymbol{R}}\equiv\widehat{\boldsymbol{\Sigma}}-\widehat{\boldsymbol{L}}$.

• •

Input the $P$-values $\pi^{*}_{1},\ldots,\pi^{*}_{p}$ to the OHC procedure (17) to get a test statistic.

Since the variance of $X_{j}^{*}$ is smaller than the variance of $X_{j}$, we expect that feeding to OHC with these new $P$-values $\pi_{j}^{*}$ will lead to more testing power. Now, let’s consider the choice of the DD-PCA algorithm. We shall use the iterative projection algorithm to be introduced in Section 4.2. The one-step DD-PCA in Section 1 is equivalent to running this iterative algorithm with only one iteration. For the sake of constructing the DD-HC test statistic, we need more iterations. By design of the iterative algorithm, as the number of iterations increases, $\|\widehat{\boldsymbol{\Sigma}}-\widehat{\boldsymbol{L}}-\widehat{\boldsymbol{A}}\|_{F}$ continues to decrease. It indicates that $\widehat{\boldsymbol{R}}$ and $\widehat{\boldsymbol{A}}$ are closer to each other, and so $\widehat{\boldsymbol{R}}$ becomes more diagonal dominant (note that $\widehat{\boldsymbol{A}}$ is forced to be in the diagonal dominant cone). The matrix $\widehat{\boldsymbol{R}}$ approximately captures the dependence structure in $X^{*}$. When $\widehat{\boldsymbol{R}}$ is close to being diagonal, it means the dependence among original $z$-scores has been fully utilized and there is minimal loss by ignoring the dependence among entries of $X_{j}^{*}$.

Remark. We are not the first to consider using a factor covariance structure to adjust $P$-values in multiple testing. Fan et al. (2012) proposed a similar idea, but their $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ are from the classical PCA. Our innovation is two fold. First, we are the first to incorporate a factor covariance structure in testing against the global null hypothesis. The main focus of Fan et al. (2012) is on estimating the false discovery proportion (FDP). Although they briefly mentioned the idea of computing factor-adjusted $P$-values, they didn’t use it in their method (their FDP estimator is still based on the original $P$-values). Second, our method is equipped with the fresh DD-PCA algorithm, in contrast with the classical PCA used in Fan et al. (2012).

We investigate the performance of new tests on extensive simulations. Given $(n,p,s,\tau)$, let $\mu_{j}=\tau\cdot 1\{1\leq j\leq s\}$ for $1\leq j\leq p$. Generate the matrix $\boldsymbol{\Sigma}=\boldsymbol{F}\boldsymbol{F}^{T}+\boldsymbol{A}$, where $\boldsymbol{F}$ is a $p\times 2$ matrix whose entries are $i.i.d$ drawn from $\mathcal{N}(0,1/2)$ and $A_{i,j}=0.5^{|i-j|}$ for $1\leq i,j\leq p$. In the null and alternative hypothesis, we generate $X_{1},X_{2},\ldots,X_{n}$ $iid$ from ${\cal N}_{p}(0,\boldsymbol{\Sigma})$ and ${\cal N}_{p}(\mu,\boldsymbol{\Sigma})$, respectively. The vector of $z$-scores is $X=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}X_{i}$, and $\widehat{\boldsymbol{\Sigma}}$ is chosen as the sample covariance matrix of $X_{i}$’s. For each test, we record the Ideal Testing error, defined as the sum of type I and type II errors with the optimal cut-off value for this test (computed via $1000$ repetitions). We fix $n=50$ and let $(p,s,\tau)$ take different values, where $s$ controls the sparsity level and $\tau$ controls the signal strength.

We compare our test with the $\chi^{2}$-test (test statistic: $\|X\|^{2}$), maximum test (test statistic: $\max_{1\leq j\leq p}|X_{j}|$), the HC test, and the Innovated HC test ($\widehat{\boldsymbol{\Omega}}$ is taken as the generalized inverse of sample covariance matrix). The results are displayed in Figure 6. We have several observations. First, the two proposed tests, IHC-DD and DD-HC, significantly outperform the other tests. Especially, the DD-HC test yields the lowest error in almost all settings. A possible reason is that DD-HC is customized to the factor covariance structure. Second, the IHC-DD test is much better than the IHC test. Since IHC-DD is a variant of IHC by plugging in $\widehat{\boldsymbol{\Sigma}}_{ddpca}^{-1}$, for a fair comparison, we include two other variants of IHC by plugging $\widehat{\boldsymbol{\Sigma}}_{poet}^{-1}$ and $\widehat{\boldsymbol{\Sigma}}_{Ind}^{-1}$, respectively, where $\widehat{\boldsymbol{\Sigma}}_{poet}$ is the POET estimator and $\widehat{\boldsymbol{\Sigma}}_{Ind}$ is a special case of POET by keeping only the diagonal entries after PCA. Interestingly, these two variants behave similarly as IHC, without much improvement. These observations suggest that the IHC idea is still quite powerful for factor covariance structure, except that we need to plug in a good estimate of $\boldsymbol{\Sigma}^{-1}$; they also reconfirm that the DD-PCA covariance estimator does a good job in estimating $\boldsymbol{\Sigma}^{-1}$. Last, the $\chi^{2}$-test, orthodox HC, and maximum test perform unsatisfactorily. Since $\boldsymbol{\Sigma}$ is heavily non-sparse, these tests lose power due to not exploring the covariance structure; comparably, the maximum test is less affected.

4 Optimization for DD-PCA

This section studies the optimization problem for DD-PCA. Section 4.1 proposes a three-block ADMM with provable theoretical guarantees to solve the convex relaxation of DD-PCA. Section 4.2 proposes an iterative projection algorithm to directly solve the nonconvex optimization of DD-PCA. Section 4.3 gives a comparison between two approaches.

Before proceeding, we introduce the efficient projection onto $\mathcal{SDD}_{c}^{+}=\mathcal{S}\cap\mathcal{DD}_{c}^{+}$, the set of “symmetric $c$-diagonally-dominant” matrices, where $\mathcal{S}$ is the set of symmetric matrices and $\mathcal{DD}_{c}^{+}$ is the set of $c$-diagonally-dominant matrices with nonnegative diagonal entries:

[TABLE]

It is not difficult to see that both $\mathcal{DD}_{c}^{+}$ and $\mathcal{SDD}_{c}^{+}$ are closed and convex polyhedral cones. To solve DD-PCA, we shall obtain the (Euclidean) projection of a matrix $\boldsymbol{A}$ onto the convex cone $\mathcal{SDD}_{c}^{+}$ or $\mathcal{DD}_{c}^{+}$, denoted by $\mathcal{P}_{\mathcal{SDD}_{c}^{+}}(\boldsymbol{A})$ or $\mathcal{P}_{\mathcal{DD}_{c}^{+}}(\boldsymbol{A})$. As summarized in Appendix A, the Mendoza-Raydan-Tarazaga (MRT) algorithm computes the efficient projection $\mathcal{P}_{\mathcal{DD}_{c}^{+}}(\boldsymbol{A})$. Following Theorem 2.1 of Mendoza, Raydan, and Tarazaga (1998), we have the convergence guarantee that $\boldsymbol{X}$ obtained by MRT Algorithm is the unique projection of $\boldsymbol{A}$ onto $\mathcal{DD}_{c}^{+}$. The computational complexity of MRT Algorithm is $O(p^{2}\log(p))$.

4.1 Convex relaxation and ADMM

This subsection solves the convex relaxation of (4) by replacing nonconvex rank constraints with convex nuclear norm constraints. To be specific, we consider the convex optimization:

[TABLE]

where $\lambda$ is a tuning parameter to strike a balance between the approximation error and the low-rank. A large $\lambda$ would encourage the solution $\hat{\boldsymbol{L}}$ to be low-rank, whereas a smaller $\lambda$ would lead to relatively smaller approximation error but higher rank in $\hat{\boldsymbol{L}}$.

We introduce a new variable $\boldsymbol{E}$ and rewrite the optimization problem as follows:

[TABLE]

The objective function would be separable in three blocks, subject to an equality constraint. Now, we define the following augmented Lagrange function:

[TABLE]

where $\boldsymbol{\Lambda}$ is the Lagrange multiplier associated with the equality constraint, and $\rho$ is a given penalty parameter. The proposed three-block ADMM proceeds as follows till convergence:

[TABLE]

Each subproblem can be efficiently solved. In the $\boldsymbol{L}$ step, we solve $\boldsymbol{L}^{(t)}$ from

[TABLE]

We have $\boldsymbol{L}^{(t)}=\mathcal{D}_{\rho^{-1}\lambda}\left(\boldsymbol{S}-\boldsymbol{A}^{(t-1)}-\boldsymbol{E}^{(t-1)}-\rho^{-1}\boldsymbol{\Lambda}^{(t-1)}\right)$, where $\mathcal{D}_{\tau}$ is the singular value thresholding operator $\mathcal{D}_{\tau}(\boldsymbol{\Omega})=\boldsymbol{U}s_{\tau}(\boldsymbol{D})\boldsymbol{V}^{T}$ for any singular value decomposition $\boldsymbol{\Omega}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^{T}$, and $s_{\tau}$ denotes the soft-thresholding operator $s_{\tau}(x)=\mbox{sgn}(x)\max(|x|-\tau,0)$.

In the $\boldsymbol{A}$ step, we need to obtain the projection on $\mathcal{SDD}_{c}^{+}$:

[TABLE]

To this end, we follow Mendoza, Raydan, and Tarazaga (1998) to use Dykstra’s alternating projection algorithm between $\mathcal{DD}_{c}^{+}$ and $\mathcal{S}$. The details of this alternating projection are included in Appendix A. Alternatively, we may follow the proximal-gradient-based ADMM (Ma et al., 2013) to solve the $\boldsymbol{A}$ step. See Section 4 of Ma et al. (2013) for more details.

In the $\boldsymbol{E}$ step, it is straightforward to solve

[TABLE]

Hence, the proposed three-block ADMM can be summarized in Algorithm 1.

Algorithm 1**.**

ADMM for Solving the Convex Relaxation of DD-PCA

Given a sample covariance matrix $\boldsymbol{S}$, do:

• •

Let $\boldsymbol{A}^{(0)}=\boldsymbol{E}^{(0)}=\boldsymbol{\Lambda}^{(0)}=\boldsymbol{0}$.

• •

For $t=1,2,\dots$

• –

$\boldsymbol{L}^{(t)}=\mathcal{D}_{\rho^{-1}\lambda}\left(\boldsymbol{S}-\boldsymbol{A}^{(t-1)}-\boldsymbol{E}^{(t-1)}-\rho^{-1}\boldsymbol{\Lambda}^{(t-1)}\right)$* where $\mathcal{D}_{\tau}(\boldsymbol{\Omega})$ is the singular value thresholding operator given by $\mathcal{D}_{\tau}(\boldsymbol{\Omega})=\boldsymbol{U}s_{\tau}(\boldsymbol{D})\boldsymbol{V}^{T}$ for any singular value decomposition $\boldsymbol{\Omega}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^{T}$, and $s_{\tau}$ denotes the soft-thresholding operator given by $s_{\tau}(x)=\mbox{sgn}(x)\max(|x|-\tau,0)$.*

• –

$\boldsymbol{A}^{(t)}=\mathcal{P}_{\mathcal{SDD}_{c}^{+}}\left(\boldsymbol{S}-\boldsymbol{L}^{(t)}-\boldsymbol{E}^{(t-1)}-\rho^{-1}\boldsymbol{\Lambda}^{(t-1)}\right)$.

• –

$\boldsymbol{E}^{(t)}=\frac{\rho}{\rho+1}\left(\boldsymbol{S}-\boldsymbol{A}^{(t)}-\boldsymbol{L}^{(t)}-\rho^{-1}\boldsymbol{\Lambda}^{(t-1)}\right)$.

• –

$\boldsymbol{\Lambda}^{(t)}=\boldsymbol{\Lambda}^{(t-1)}+\rho\left(\boldsymbol{A}^{(t)}+\boldsymbol{L}^{(t)}+\boldsymbol{E}^{(t)}-\boldsymbol{S}\right)$.

• •

Stop if the convergence criterion is met.

Although three-block ADMM does not necessarily converge in general (Chen et al., 2016), DD-PCA belongs to a class of regularized least squares decomposition problem. For this class of regularized problems, the global convergence of the proposed three-block ADMM is always guaranteed such that any cluster point of the iterated solutions is an optimal primal and dual pair of DD-PCA (See Theorem 3.2 of Lin et al. (2018)).

4.2 An iterative projection algorithm

In the sequel, we introduce an iterative projection algorithm that directly tackles the nonconvex optimization in DD-PCA. The key observation is that we attempt to find a matrix $\boldsymbol{L}^{\ast}$ in the set $\mathcal{L}_{K}=\{\boldsymbol{L}:\text{rank}(\boldsymbol{L})=K\}$ that is closest to the set $\mathcal{M}_{\boldsymbol{S}}=\{\boldsymbol{S}-\boldsymbol{A}:\boldsymbol{A}\in\mathcal{SDD}_{c}^{+}\}$. Inspired by Netrapalli et al. (2014), a natural approach would be to iteratively project $(\boldsymbol{S}-\boldsymbol{L})$ onto $\mathcal{SDD}_{c}^{+}$ to update $\boldsymbol{A}$ and then to project $(\boldsymbol{S}-\boldsymbol{A})$ onto $\mathcal{L}_{K}$ to update $\boldsymbol{L}$. To reduce the computational cost, we replace the projection onto $\mathcal{SDD}_{c}^{+}$ by the projection onto $\mathcal{DD}_{c}^{+}$, followed by symmetrization. Algorithm 2 summarizes the details.

Algorithm 2**.**

Iterative Projection Algorithm for Solving the DD-PCA

Given a sample covariance matrix $\boldsymbol{S}$ and integer $k$, do:

• •

Let $\boldsymbol{A}^{(0)}=\mathbf{0}$.

• •

For $t=1,2,\dots$

• –

$\boldsymbol{L}^{(t)}=\mathcal{P}_{\mathcal{L}_{K}}(\boldsymbol{S}-\boldsymbol{A}^{(t-1)})$**

• –

$\tilde{\boldsymbol{A}}^{(t)}=\mathcal{P}_{\mathcal{DD}_{c}^{+}}(\boldsymbol{S}-\boldsymbol{L}^{(t)})$.

• –

$\boldsymbol{A}^{(t)}=\left(\tilde{\boldsymbol{A}}^{(t)}+(\tilde{\boldsymbol{A}}^{(t)})^{T}\right)/2$.

• •

Stop if the convergence criterion is met.

In Algorithm 2, we need to calculate $\mathcal{P}_{\mathcal{L}_{K}}$ and $\mathcal{P}_{\mathcal{DD}_{c}^{+}}$. The calculation of $\mathcal{P}_{\mathcal{DD}_{c}^{+}}$ is given in Appendix A. The calculation of $\mathcal{P}_{\mathcal{L}_{K}}$ is given as follows: for any symmetric matrix $\boldsymbol{A}$, we write its eigenvalue decomposition as $\boldsymbol{A}=\boldsymbol{Q}\boldsymbol{\Lambda}\boldsymbol{Q}^{T}$ where $\boldsymbol{\Lambda}=\mathrm{diag}\{\lambda_{1},\lambda_{2},\dots,\lambda_{p}\}$ with $|\lambda_{1}|\geq|\lambda_{2}|\geq\dots\geq|\lambda_{p}|$. Hence, the best rank-$K$ approximation is given by $\mathcal{P}_{\mathcal{L}_{K}}(\boldsymbol{A})=\boldsymbol{Q}_{K}\boldsymbol{\Lambda}_{K}\boldsymbol{Q}_{K}^{T}$ where $\boldsymbol{Q}_{K}$ contains the first $K$ columns of $\boldsymbol{Q}$ and $\boldsymbol{\Lambda}_{K}=\mathrm{diag}\{\lambda_{1},\lambda_{2},\dots,\lambda_{K}\}$.

To use Algorithm 2, we need to estimate the rank $K$ if it is unknown. A simple estimate of $K$ is to look at the eigenvalues of $\boldsymbol{S}$ to pick the $K$ such that there is a significant gap in magnitude between the first $K$ eigenvalues and the remaining ones. In Section 5, we investigate the robustness of the iterative projection algorithm to the estimation of $K$.

4.3 Comparing convex and nonconvex approaches

The convex approaches do not require the knowledge of rank $K$ of the low rank matrix $\boldsymbol{L}$, and the global convergence of the proposed ADMM is guaranteed. However, its convergence rate could be slow. The nonconvex approaches, on the other hand, can be faster in terms of convergence. The per-iteration cost for Algorithm 2 is $O(p^{2}\max\{\log(p),K\})$, compared to $O(p^{3})$ for Algorithm 1. But the convergence guarantee of Algorithm 2 is an open question.

The rigorous convergence analysis of the iterative projection algorithm is difficult due to the non-convexity of the set $\mathcal{L}_{K}$. The existing result (e.g., Drusvyatskiy et al. (2015)) proves the local linear convergence of the alternating projections for two closed sets if the two sets intersect transversally at the converging point. We conjecture that such condition would hold for most cases in our setting, therefore the convergence would be guaranteed. In practice, our algorithms are stable and always converge to a valid solution in simulations.

5 Simulation studies

This section investigates several numerical properties of DD-PCA, including the estimation performance, necessity and robustness, and application to covariance matrix estimation.

Experiment 1: Exact DD-PCA.

We first examine the numerical performance of two-block ADMM (see Appendix B) and the iterative projection algorithm in Algorithm 2 to solve the exact DD-PCA. Fixing $(p,K)$, we first generate a rank-$K$ matrix $\boldsymbol{L}=XX^{T}$ where $X$ is a $p\times k$ matrix whose entries are i.i.d drawn from $\mathcal{N}(0,1/p)$. We then generate a matrix $\boldsymbol{A}_{0}$ with entries sampled i.i.d from $\mathcal{N}(0,1/p^{2})$ and set $\boldsymbol{A}=\boldsymbol{A}_{0}+\boldsymbol{A}_{0}^{T}+\boldsymbol{D}$, where $\boldsymbol{D}$ is a diagonal matrix whose $j$-th diagonal is equal to $\sum_{i:i\neq j}|A_{0}(j,i)+A_{0}(i,j)|-2A_{0}(j,j)$ for $1\leq j\leq p$; it follows that $\boldsymbol{A}$ is a diagonally dominant matrix. We then let $\boldsymbol{S}=\boldsymbol{L}+\boldsymbol{A}$. We consider $p=500,1000,2000$ and fix $K=0.05\cdot p$ for each choice of $p$.

First, for two-block ADMM, we use the solution $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ after $20$ iterations. Table 1 displays the comparison between $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ and the true $(\boldsymbol{L},\boldsymbol{A})$ based on $20$ repetitions. It suggests that $\widehat{\boldsymbol{L}}$ always has the same rank as that of $\boldsymbol{L}$ and that $\widehat{\boldsymbol{L}}$ and $\widehat{\boldsymbol{A}}$ are reasonably close to their respective counterparts. Since the ADMM algorithm is an iterative algorithm, $\widehat{\boldsymbol{L}}+\widehat{\boldsymbol{A}}$ are not exactly equal to $\boldsymbol{S}$, but the two matrices are reasonably close after $20$ iterations. The results also suggest that, for a large $p$, more iterations are needed for the convergence of ADMM.

Next, we look at Algorithm 2. Instead of fixing the maximum number of iterations, we investigate how the solution $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ evolves over iterations. Note that $\widehat{\boldsymbol{L}}$ is guaranteed to have rank $K$ in all iterations, and $\widehat{\boldsymbol{A}}$ is equal to the projection of $(\boldsymbol{S}-\widehat{\boldsymbol{L}})$ into ${\cal DD}^{+}$ (with symmetrization). We introduce a quantity to measure the diagonal dominance of $(\boldsymbol{S}-\widehat{\boldsymbol{L}})$. For any matrix $p\times p$ matrix $\boldsymbol{B}$, define

[TABLE]

It measures how close a matrix is to the diagonally dominant cone. If $\zeta(\boldsymbol{\Sigma}-\widehat{\boldsymbol{L}})$ continues to increase and eventually gets close to zero, then the algorithm converges. The left panel of Figure 7 shows the evolution of $\zeta(\boldsymbol{S}-\widehat{\boldsymbol{L}})$ over iterations, and it suggests that the algorithm converges quickly. The right panel of Figure 7 displays the evolution of the relative approximation error $\|\widehat{\boldsymbol{L}}+\widehat{\boldsymbol{A}}-\boldsymbol{S}\|_{F}/\|\boldsymbol{S}\|_{F}$, which also decreases quickly over iterations.

Experiment 2: Approximate DD-PCA.

We investigate the performance of Algorithm 1 (an ADMM algorithm) and Algorithm 2 (an iterative projection algorithm) for approximate DD-PCA. Fixing $(p,K)$ and $\sigma>0$, we generate a rank $K$ matrix $\boldsymbol{L}$ and a diagonally dominant matrix $\boldsymbol{A}$ in the same way as in Experiment 1. We then generate a $p\times p$ symmetric matrix $\boldsymbol{E}$ whose upper triangular entries are sampled i.i.d from $\mathcal{N}(0,\sigma^{2}/p)$. Last, let $\boldsymbol{S}=\boldsymbol{L}+\boldsymbol{A}+\boldsymbol{E}$.

First, we study Algorithm 1, which is an ADMM algorithm. Fix $\sigma=1$. We consider $p=500,1000,2000$, and set $K=0.05\cdot p$. The tuning parameter in the algorithm is set as $\lambda=3$, and we look at the solution $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ after 50 iterations. The results are displayed in Table 2. For all three settings, the algorithm exactly recovers the true rank of $\boldsymbol{L}$, however, the convergence of $(\widehat{\boldsymbol{L}},\widehat{\boldsymbol{A}})$ is relatively slow. As we shall see below, the performance of Algorithm 1 is not as good as the iterative projection algorithm—Algorithm 2, but Algorithm 1 is theoretically more tractable.

Next, we study Algorithm 2, the iterative projection algorithm. In Experiment 1, we have investigated its performance when $\boldsymbol{S}$ has an exact decomposition to the sum of a low-rank matrix and a diagonally dominant matrix. In this experiment, we apply the same algorithm to $\boldsymbol{S}$ which does not have such an exact decomposition. We run the algorithm for $20$ iterations and measure the relative approximation error $\|\widehat{\boldsymbol{L}}+\widehat{\boldsymbol{A}}-\boldsymbol{S}\|_{F}/\|\boldsymbol{S}\|_{F}$. The results are shown in Figure 8. In the left panel, $K/p$ is fixed as $0.05$ and the noise level $\sigma$ varies from $0.5$ to $5$. In the right panel, $\sigma$ is fixed to be $1$ and $K/p$ varies from $0.01$ to $0.1$. For each value of $p$, the relative approximation error increases, as both $\sigma$ and $K$ increase. For the same values of $\sigma$ and $K/p$, a larger $p$ comes with a smaller relative approximation error. Furthermore, if we compare the results with those in Table 2, Algorithm 2 has a better practical performance than Algorithm 1.

Experiment 3: Necessity of DD-PCA.

If $\boldsymbol{\Sigma}$ truly satisfies the assumption of “low-rank plus diagonal dominance”, it is a natural question to know whether one can simply apply PCA and robust PCA (Candès et al., 2011) to get a diagonally dominant $\boldsymbol{A}$. Unfortunately, this is often not the case. Let us consider applying PCA to a $\boldsymbol{\Sigma}$ which has the decomposition $\boldsymbol{\Sigma}=\boldsymbol{L}_{0}+\boldsymbol{A}_{0}$ such that $\mathrm{rank}(\boldsymbol{L}_{0})=K$ and $\boldsymbol{A}_{0}$ is diagonally dominant. Let $\lambda_{k}$ and $\xi_{k}$ be the $k$-th eigenvalue and eigenvector, respectively, $1\leq k\leq p$. We construct $\boldsymbol{L}=\sum_{k=1}^{K}\lambda_{k}\xi_{k}\xi_{k}^{T}$ and $\boldsymbol{A}=\boldsymbol{\Sigma}-\boldsymbol{L}$. We can only hope $\boldsymbol{A}$ is diagonally dominant when $\boldsymbol{A}$ and $\boldsymbol{A}_{0}$ are entrywise close to each other, or equivalently, when $\|\boldsymbol{L}-\boldsymbol{L}_{0}\|_{\max}$ is small ($\|\cdot\|_{\max}$ is the entrywise max norm). However, from the literatures on perturbation analysis of PCA, it requires strong conditions to guarantee that $\|\boldsymbol{L}-\boldsymbol{L}_{0}\|_{\max}$ is small (Ke and Yang, 2017). In particular, when $K$ is moderately large, these conditions may be violated. Similarly, robust PCA cannot produce a diagonally dominant $\boldsymbol{A}$ in general. Therefore, it is necessary to develop new algorithms that are specifically designed for DD-PCA. In Figure 9, we present a numerical example, where the output $\boldsymbol{A}$ from our DD-PCA algorithm is much more “diagonally dominant” than the $\boldsymbol{A}$ constructed from PCA.

Experiment 4: Robustness to the misspecification of

$K$.

We use the same setup as in Experiment 1 and investigate the performance of Algorithm 2 with a misspecified $K$. Consider two settings where $(p,K)=(500,25)$ and $(p,K)=(2000,100)$, respectively. For each setting, we plug $K=k$ in the algorithm and take the solution after ten iterations. Figure 10 shows the relative difference between $\widehat{\boldsymbol{L}}$ and $\boldsymbol{L}$ for various choices of $k$. It suggests that as long as $k\geq K$, the performance of the algorithm is very stable. Hence, in practice, we recommend that the users pick a relatively large $k$ when the true $K$ is hard to estimate.

Experiment 5: Application to covariance matrix estimation

We expand the numerical study in Section 2 and investigate the performance of DD-POET on more simulation settings. Given $K=3$ and $p\in\{100,300,500\}$, we generate data in the same way as in the numerical example of Section 2. First, we compare the performance of DD-POET and POET. For both methods, we use the true $K=3$. POET has an additional threshold, which we set as the ideal one that minimizes the estimation error (the ideal threshold varies as we change the error measure). The results are contained in Column 6 and Column 10 of Table 3, where, in all settings, DD-POET has a comparable performance as POET with an ideal threshold, and in some settings, DD-POET is even better. The ideal threshold for POET is not practically feasible, and it is unclear how to set the threshold in a data-driven fashion; however, DD-POET is tuning free once $K$ is given. Second, we investigate the performance of DD-POET when we plug in $K=k$ with $k\in\{1,2,\ldots,6\}$; see Table 3. If $k$ is misspecified but $k\geq K$, the estimation errors remain relatively stable; if $k<K$, the performance deteriorates. It suggests that an overshooting of $K$ is better than an undershooting. This is consistent with the observations made by Fan et al. (2013a).

6 Discussion

The diagonally dominant matrices have been well studied in linear algebra (Feingold and Varga, 1962; Barker and Carlson, 1975) and optimization (Barlow and Demmel, 1990; Mendoza et al., 1998), motivated by the appealing properties of these matrices for computation. Our work has a very different motivation: We recognize that diagonally dominant (covariance) matrices also have appealing statistical properties, and propose exploring the “low-rank plus diagonal dominance” covariance structure in data analysis. We demonstrate the benefit of exploring such structure in two statistical problems. For covariance matrix estimation, we propose DD-PCA as a new estimator. For testing of the global null hypothesis in multiple testing, we propose IHC-DD and DD-HC as two new tests. These new methods have shown encouraging numerical performance in simulations and real applications, especially when the data have a factor covariance structure.

The above methods rely on the availability of algorithms to decompose any given covariance matrix (approximately) into the sum of a low-rank matrix and a diagonally dominant matrix. To obtain such decomposition is a non-convex optimization. We propose two algorithms — an ADMM algorithm that solves a convex relaxation, and an iterative projection algorithm that solves the noncovex problem directly. In comparison, the ADMM algorithm is theoretically more tractable, and the iterative projection algorithm shows very appealing numerical performance.

The study here motivates several interesting future directions, such as the uniqueness of the low-rank plus diagonal dominance decomposition, the statistical error of recovering $\boldsymbol{L}$ and $\boldsymbol{A}$, as well as the convergence of proposed algorithms. We leave to future works.

Our work is related to the literatures of factor models and the literatures of “low-rank plus sparse” matrix decomposition. These two lines of works have found wide applications in many areas. Similarly, the use of DD-PCA is not limited to covariance matrix estimation and multiple testing. We expect that DD-PCA will have applications in classification (Zhu and Hastie, 2005; Tong et al., 2016), clustering (Wang and Zhu, 2008; Clarke et al., 2009), dimension reduction (Ma and Zhu, 2013; Zou and Xue, 2018) and forecasting (Fan et al., 2017).

Appendix A Efficient projection onto $\mathcal{SDD}_{c}^{+}$

Recall that $\mathcal{S}$ is the set of symmetric matrices and $\mathcal{DD}_{c}^{+}$ is the set of $c$-diagonally-dominant matrices with nonnegative diagonal entries. Now, we present the (Euclidean) projection of a matrix $\boldsymbol{A}$ onto the convex cone $\mathcal{SDD}_{c}^{+}$ or $\mathcal{DD}_{c}^{+}$, denoted by $\mathcal{P}_{\mathcal{SDD}_{c}^{+}}(\boldsymbol{A})$ or $\mathcal{P}_{\mathcal{DD}_{c}^{+}}(\boldsymbol{A})$.

Algorithm 3**.**

Mendoza-Raydan-Tarazaga (MRT) Algorithm

Given a $p\times p$ matrix $\boldsymbol{A}$, where the $j$th row of $\boldsymbol{A}$ is denoted by $\mathbf{a}_{j}$. For $1\leq j\leq p$, the $j$th row of the projection $\boldsymbol{X}$, denoted by $\mathbf{x}_{j}$, is given by

• •

If $a_{jj}\geq\sum_{l:l\not=j}|a_{jl}|$, then $\mathbf{x}_{j}=\mathbf{a}_{j}$.

• •

If $-\sum_{l:l\not=j}|a_{jl}|\leq a_{jj}<0$ and $|a_{jj}|>|a_{jl}|$ for all $l\not=j$, or $a_{jj}<-\sum_{l:l\not=j}|a_{jl}|$, then $\mathbf{x}_{j}=\mathbf{0}$.

• •

If $-\sum_{l:l\not=j}|a_{jl}|\leq a_{jj}<0$ and $|a_{jj}|\leq|a_{jl}|$ for some $l\not=j$, or $0\leq a_{jj}<\sum_{l:l\not=j}|a_{jl}|$, then $\mathbf{x}_{j}$ is generated as follows:

Sort $|\mathbf{a}_{j}|$, excluding $a_{jj}$, in the ascending order, and denote the reordered vector as $e$. Note that $e_{j}=a_{jj}$ and $|e_{i}|\leq|e_{l}|$ for all $i<l,i\not=j,l\not=j$. 2. 2.

For $m\not=j$, compute $d_{m}=\sum_{l=m}^{p}|e_{l}|\cdot I_{\{j\neq l\}}-e_{j}$ and $\bar{d}_{m}=d_{m}/(p-m+1)\cdot I_{\{m<j\}}+d_{m}/(p-m+2)\cdot I_{\{m>j\}}$ 3. 3.

Solve $m^{\star}$ as the smallest integer among $m=1,\ldots,p$ such that $m\neq j$, $|e_{m}|>0$ and $|e_{m}|\geq\bar{d}_{m}$ 4. 4.

Solve $\mathbf{x}_{j}=(x_{j1},\dots,x_{jp})$ such that $x_{jj}=a_{jj}+\bar{d}_{m^{\ast}}$; $x_{ji}=(a_{ji}-\bar{d}_{m^{\ast}})^{+}$ if $a_{ji}\geq 0$ for $i\not=j$; $x_{ji}=-(a_{ji}+\bar{d}_{m^{\ast}})^{-}$ if $a_{ji}<0$ for $i\not=j$, where $(z)^{+}=\max\{z,0\}$ and $(z)^{-}=-\min\{z,0\}$.

Mendoza et al. (1998) applied Dykstra’s alternating projection algorithm between $\mathcal{DD}^{+}$ and $\mathcal{S}$ to obtain the projection on $\mathcal{SDD}^{+}$. The algorithm is summarized in Algorithm 4.

Algorithm 4**.**

Efficient Projection onto $\mathcal{SDD}^{+}$

Given a $p\times p$ matrix $\boldsymbol{A}$,

• •

Let $\boldsymbol{G}^{(0)}=\boldsymbol{A}$ and $\boldsymbol{I}^{(0)}=\mathbf{0}$

• •

For $t=1,2,\dots$

• –

$\boldsymbol{G}^{(t)}=\mathcal{P}_{\mathcal{DD}^{+}}\left(\frac{1}{2}(\boldsymbol{G}^{(t-1)}+(\boldsymbol{G}^{(t-1)})^{T})-\boldsymbol{I}^{(t-1)}\right)$**

• –

$\boldsymbol{I}^{(t)}=\boldsymbol{G}^{(t)}-\left(\frac{1}{2}(\boldsymbol{G}^{(t-1)}+(\boldsymbol{G}^{(t-1)})^{T})-\boldsymbol{I}^{(t-1)}\right)$**

• •

Stop if the convergence criterion is met.

When $c=1$, the convergence result of Algorithm 4 can be similarly established as in Boyle and Dykstra (1986) such that the iterated solutions converge in the Frobenius norm to the unique solution of the projection on $\mathcal{SDD}^{+}$. More details can be found in Mendoza et al. (1998). When $c\not=1$, MRT algorithm can’t be directly used. In this case, we obtain $\mathcal{P}_{\mathcal{DD}_{c}^{+}}(\boldsymbol{A})$ through Quadratic Programming (QP). The key observation is that the problem can be separated as $p$ independent row-wise projection. For each $1\leq j\leq p$, the $j$th row projection can be written as

[TABLE]

and the solution $(v_{1},\dots,v_{p})$ would be the $j$th row of $\mathcal{P}_{\mathcal{DD}_{c}^{+}}(\boldsymbol{A})$. We can reformulate (20) as

[TABLE]

It’s easy to see that for $i\not=j$, we should let $\text{sign}(\delta_{i})=\text{sign}(a_{ji})$ and $|\delta_{i}\leq a_{ji}|$, and hence $|a_{ji}-\delta_{i}|=|a_{ji}|-|\delta_{i}|$. Without loss of generality, we assume $a_{ji}\geq 0$ for all $i\not=j$ so we can restrict $\delta_{i}\geq 0$ for all $i\not=j$. Then (21) becomes

[TABLE]

which is a QP problem and can be solved using standard solver.

Appendix B Convex relaxation and ADMM for Exact DD-PCA

The exact DD-PCA is difficult to solve due to the nonconvex rank minimization. Consider the following convex relaxation of the exact DD-PCA:

[TABLE]

where $\|\cdot\|_{*}$ is the matrix nuclear norm.

Given the efficient projection onto $\mathcal{DD}^{+}$ in Algorithm 3, we introduce a new variable $\boldsymbol{B}$, satisfying the equality that $\boldsymbol{A}=\boldsymbol{B}$, to separate the symmetric and diagonally-dominant constraints as follows:

[TABLE]

where $\mathcal{I}_{C}$ is the indicator function which equals to 0 if condition $C$ is satisfied, and equals to infinity otherwise (Boyd and Vandenberghe, 2004).

We define the following augmented Lagrange function:

[TABLE]

where $\boldsymbol{\Lambda}_{1}$ and $\boldsymbol{\Lambda}_{2}$ are the Lagrangian multipliers associated with the equality constraints, and $\rho$ is a given penalty parameter. We propose an efficient ADMM to solve the exact DD-PCA from $\mathcal{L}_{\rho}(\boldsymbol{L},\boldsymbol{A},\boldsymbol{B},\boldsymbol{\Lambda}_{1},\boldsymbol{\Lambda}_{2})$, which proceeds as follows till convergence:

[TABLE]

Our proposed ADMM is a two-block ADMM with two blocks $\{\boldsymbol{L},\boldsymbol{B}\}$ and $\boldsymbol{A}$, and its global convergence is always guaranteed (Boyd et al., 2011). In what follows, we explicitly show how to obtain closed-form solutions for each subproblem. In the $\boldsymbol{L}$ step, we have

[TABLE]

It’s easy to show that the solution is given by $\boldsymbol{L}^{(t)}=\mathcal{D}_{\rho^{-1}}\left(\boldsymbol{S}-\boldsymbol{A}^{(t-1)}-\rho^{-1}\boldsymbol{\Lambda}_{1}^{(t-1)}\right)$ where $\mathcal{D}_{\tau}(\boldsymbol{\Omega})$ is the singular value thresholding operator given by $\mathcal{D}_{\tau}(\boldsymbol{\Omega})=\boldsymbol{U}s_{\tau}(\boldsymbol{D})\boldsymbol{V}^{T}$ for any singular value decomposition $\boldsymbol{\Omega}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^{T}$, and $s_{\tau}$ denotes the soft-thresholding operator given by $s_{\tau}(x)=\mbox{sgn}(x)\max(|x|-\tau,0)$.

In the $\boldsymbol{B}$ step, we also have the following closed-form solution:

[TABLE]

Finally in the $\boldsymbol{A}$ step, we have

[TABLE]

We summarize our proposed two-block ADMM in Algorithm 5.

Algorithm 5**.**

Two-Block ADMM for Solving the Exact DD-PCA

Given the sample covariance matrix $\boldsymbol{S}$, do

• •

Let $\boldsymbol{A}^{(0)}=\boldsymbol{\Lambda}_{1}^{(0)}=\boldsymbol{\Lambda}_{2}^{(0)}=\boldsymbol{0}$

• •

For $t=1,2,\dots$

• –

$\boldsymbol{L}^{(t)}=\mathcal{D}_{\rho^{-1}}\left(\boldsymbol{S}-\boldsymbol{A}^{(t-1)}-\rho^{-1}\boldsymbol{\Lambda}_{1}^{(t-1)}\right)$.

• –

$\boldsymbol{B}^{(t)}=\frac{1}{2}\left[\left(\boldsymbol{A}^{(t-1)}+\rho^{-1}\boldsymbol{\Lambda}_{2}^{(t-1)}\right)+\left(\boldsymbol{A}^{(t-1)}+\rho^{-1}\boldsymbol{\Lambda}_{2}^{(t-1)}\right)^{T}\right]$**

• –

$\boldsymbol{A}^{(t)}=\mathcal{P}_{\mathcal{DD}^{+}}\left(\frac{1}{2}(\boldsymbol{S}-\boldsymbol{L}^{(t)}+\boldsymbol{B}^{(t)}-\rho^{-1}\boldsymbol{\Lambda}_{1}^{(t-1)}-\rho^{-1}\boldsymbol{\Lambda}_{2}^{(t-1)})\right)$**

• –

$\boldsymbol{\Lambda}_{1}^{(t)}=\boldsymbol{\Lambda}_{1}^{(t-1)}+\rho(\boldsymbol{A}^{(t)}+\boldsymbol{L}^{(t)}-\boldsymbol{S})$**

• –

$\boldsymbol{\Lambda}_{2}^{(t)}=\boldsymbol{\Lambda}_{2}^{(t-1)}+\rho(\boldsymbol{A}^{(t)}-\boldsymbol{B}^{(t)})$**

• •

Stop if the convergence criterion is met.