Compressed Covariance Estimation With Automated Dimension Learning

TL;DR

This paper introduces a scalable method for covariance matrix estimation that combines low-rank and diagonal structures, utilizing data compression and an automated dimension learning framework based on Stein's Unbiased Risk Estimation.

Contribution

It presents a novel covariance estimation technique that does not require the low-rank component to be sparse, with an automated dimension selection method.

Findings

Effective in high-dimensional settings

Scalable and computationally efficient

Outperforms existing methods in simulations

Abstract

We propose a method for estimating a covariance matrix that can be represented as a sum of a low-rank matrix and a diagonal matrix. The proposed method compresses high-dimensional data, computes the sample covariance in the compressed space, and lifts it back to the ambient space via a decompression operation. A salient feature of our approach relative to existing literature on combining sparsity and low-rank structures in covariance matrix estimation is that we do not require the low-rank component to be sparse. A principled framework for estimating the compressed dimension using Stein's Unbiased Risk Estimation theory is demonstrated. Experimental simulation results demonstrate the efficacy and scalability of our proposed approach.

Tables

Table 1. Table 1 . Simulation Setting 1 with s = 0.5 𝑠 0.5 s=0.5 . Top panel compares k ^ sure superscript ^ 𝑘 sure \widehat{k}^{\text{sure}} with k opt superscript 𝑘 opt k^{\text{opt}} . Bottom two panels compare the proposed approach with POET and AT in terms of ‖ Σ ^ − Σ 0 ‖ 2 / p subscript norm ^ Σ subscript Σ 0 2 𝑝 \|{\widehat{\Sigma}}-{\Sigma_{0}}\|_{2}/p and ‖ Σ ^ − Σ 0 ‖ F / p subscript norm ^ Σ subscript Σ 0 𝐹 𝑝 \|{\widehat{\Sigma}}-{\Sigma_{0}}\|_{F}/p respectively. Standard errors are in parenthesis.

ktr	10			50
p	250	500	1000	250	500	1000
$k^{\text{opt}}$	240	480	950	210	420	830
${\hat{k}}^{\text{sure}}$	240	480	950	210	410	820
CD	$\underset{\left(0\cdot 04\right)}{0\cdot 25}$	$\underset{\left(0\cdot 03\right)}{0\cdot 28}$	$\underset{\left(0\cdot 03\right)}{0\cdot 27}$	$\underset{\left(0\cdot 03\right)}{0\cdot 63}$	$\underset{\left(0\cdot 02\right)}{0\cdot 57}$	$\underset{\left(0\cdot 02\right)}{0\cdot 54}$
AT	$\underset{\left(0\cdot 05\right)}{0\cdot 28}$	$\underset{\left(0\cdot 05\right)}{0\cdot 30}$	$\underset{\left(0\cdot 05\right)}{0\cdot 30}$	$\underset{\left(0\cdot 09\right)}{0\cdot 88}$	$\underset{\left(0\cdot 08\right)}{0\cdot 77}$	$\underset{\left(0\cdot 07\right)}{0\cdot 71}$
POET	$\underset{\left(0\cdot 05\right)}{0\cdot 28}$	$\underset{\left(0\cdot 05\right)}{0\cdot 30}$	$\underset{\left(0\cdot 05\right)}{0\cdot 30}$	$\underset{\left(0\cdot 09\right)}{0\cdot 93}$	$\underset{\left(0\cdot 08\right)}{0\cdot 90}$	$\underset{\left(0\cdot 07\right)}{0\cdot 92}$
CD	$\underset{\left(0\cdot 04\right)}{0\cdot 56}$	$\underset{\left(0\cdot 03\right)}{0\cdot 57}$	$\underset{\left(0\cdot 03\right)}{0\cdot 54}$	$\underset{\left(0\cdot 03\right)}{2\cdot 14}$	$\underset{\left(0\cdot 03\right)}{2\cdot 14}$	$\underset{\left(0\cdot 03\right)}{2\cdot 12}$
AT	$\underset{\left(0\cdot 06\right)}{0\cdot 63}$	$\underset{\left(0\cdot 05\right)}{0\cdot 63}$	$\underset{\left(0\cdot 05\right)}{0\cdot 59}$	$\underset{\left(0\cdot 08\right)}{3\cdot 50}$	$\underset{\left(0\cdot 08\right)}{3\cdot 56}$	$\underset{\left(0\cdot 08\right)}{3\cdot 57}$
POET	$\underset{\left(0\cdot 05\right)}{0\cdot 60}$	$\underset{\left(0\cdot 05\right)}{0\cdot 61}$	$\underset{\left(0\cdot 05\right)}{0\cdot 58}$	$\underset{\left(0\cdot 08\right)}{2\cdot 60}$	$\underset{\left(0\cdot 07\right)}{2\cdot 62}$	$\underset{\left(0\cdot 08\right)}{2\cdot 60}$

Table 2. Table 2 . Simulation Setting 1 with s = 0.1 𝑠 0.1 s=0.1 . Top panel compares k ^ sure superscript ^ 𝑘 sure \widehat{k}^{\text{sure}} with k opt superscript 𝑘 opt k^{\text{opt}} . Bottom two panels compare the proposed approach with POET and AT in terms of ‖ Σ ^ − Σ 0 ‖ 2 / p subscript norm ^ Σ subscript Σ 0 2 𝑝 \|{\widehat{\Sigma}}-{\Sigma_{0}}\|_{2}/p and ‖ Σ ^ − Σ 0 ‖ F / p subscript norm ^ Σ subscript Σ 0 𝐹 𝑝 \|{\widehat{\Sigma}}-{\Sigma_{0}}\|_{F}/p respectively. Standard errors are in parenthesis.

ktr	10			50
p	250	500	1000	250	500	1000
$k^{\text{opt}}$	240	470	950	210	420	820
${\hat{k}}^{\text{sure}}$	240	470	940	210	410	810
CD	$\underset{\left(0\cdot 06\right)}{0\cdot 49}$	$\underset{\left(0\cdot 06\right)}{0\cdot 52}$	$\underset{\left(0\cdot 06\right)}{0\cdot 50}$	$\underset{\left(0\cdot 05\right)}{1\cdot 16}$	$\underset{\left(0\cdot 04\right)}{1\cdot 03}$	$\underset{\left(0\cdot 02\right)}{0\cdot 95}$
AT	$\underset{\left(0\cdot 09\right)}{0\cdot 52}$	$\underset{\left(0\cdot 10\right)}{0\cdot 56}$	$\underset{\left(0\cdot 09\right)}{0\cdot 56}$	$\underset{\left(0\cdot 12\right)}{1\cdot 67}$	$\underset{\left(0\cdot 10\right)}{1\cdot 39}$	$\underset{\left(0\cdot 06\right)}{1\cdot 29}$
POET	$\underset{\left(0\cdot 09\right)}{0\cdot 52}$	$\underset{\left(0\cdot 09\right)}{0\cdot 56}$	$\underset{\left(0\cdot 09\right)}{0\cdot 56}$	$\underset{\left(0\cdot 18\right)}{1\cdot 74}$	$\underset{\left(0\cdot 12\right)}{1\cdot 67}$	$\underset{\left(0\cdot 14\right)}{1\cdot 61}$
CD	$\underset{\left(0\cdot 07\right)}{0\cdot 92}$	$\underset{\left(0\cdot 07\right)}{0\cdot 98}$	$\underset{\left(0\cdot 07\right)}{0\cdot 96}$	$\underset{\left(0\cdot 05\right)}{3\cdot 88}$	$\underset{\left(0\cdot 06\right)}{3\cdot 80}$	$\underset{\left(0\cdot 06\right)}{3\cdot 76}$
AT	$\underset{\left(0\cdot 10\right)}{1\cdot 01}$	$\underset{\left(0\cdot 10\right)}{1\cdot 07}$	$\underset{\left(0\cdot 10\right)}{1\cdot 05}$	$\underset{\left(0\cdot 42\right)}{6\cdot 36}$	$\underset{\left(0\cdot 36\right)}{6\cdot 30}$	$\underset{\left(0\cdot 17\right)}{6\cdot 37}$
POET	$\underset{\left(0\cdot 08\right)}{0\cdot 96}$	$\underset{\left(0\cdot 09\right)}{1\cdot 03}$	$\underset{\left(0\cdot 09\right)}{1\cdot 01}$	$\underset{\left(0\cdot 16\right)}{4\cdot 68}$	$\underset{\left(0\cdot 14\right)}{4\cdot 63}$	$\underset{\left(0\cdot 15\right)}{4\cdot 59}$

Table 3. Table 3 . Simulation Setting 2 with s = 0.5 𝑠 0.5 s=0.5 . Top panel compares k ^ sure superscript ^ 𝑘 sure \widehat{k}^{\text{sure}} with k opt superscript 𝑘 opt k^{\text{opt}} . Bottom two panels compare the proposed approach with POET and AT in terms of ‖ Σ ^ − Σ 0 ‖ 2 / p subscript norm ^ Σ subscript Σ 0 2 𝑝 \|{\widehat{\Sigma}}-{\Sigma_{0}}\|_{2}/p and ‖ Σ ^ − Σ 0 ‖ F / p subscript norm ^ Σ subscript Σ 0 𝐹 𝑝 \|{\widehat{\Sigma}}-{\Sigma_{0}}\|_{F}/p respectively. Standard errors are in parenthesis.

ktr	10			50
p	250	500	1000	250	500	1000
$k^{\text{opt}}$	240	480	950	210	420	830
${\hat{k}}^{\text{sure}}$	240	480	940	210	410	810
CD	$\underset{\left(0\cdot 04\right)}{0\cdot 31}$	$\underset{\left(0\cdot 03\right)}{0\cdot 26}$	$\underset{\left(0\cdot 03\right)}{0\cdot 28}$	$\underset{\left(0\cdot 03\right)}{0\cdot 63}$	$\underset{\left(0\cdot 02\right)}{0\cdot 57}$	$\underset{\left(0\cdot 02\right)}{0\cdot 53}$
AT	$\underset{\left(0\cdot 06\right)}{0\cdot 33}$	$\underset{\left(0\cdot 05\right)}{0\cdot 29}$	$\underset{\left(0\cdot 05\right)}{0\cdot 29}$	$\underset{\left(0\cdot 07\right)}{0\cdot 88}$	$\underset{\left(0\cdot 05\right)}{0\cdot 77}$	$\underset{\left(0\cdot 03\right)}{0\cdot 71}$
POET	$\underset{\left(0\cdot 06\right)}{0\cdot 33}$	$\underset{\left(0\cdot 05\right)}{0\cdot 29}$	$\underset{\left(0\cdot 05\right)}{0\cdot 30}$	$\underset{\left(0\cdot 08\right)}{0\cdot 93}$	$\underset{\left(0\cdot 08\right)}{0\cdot 92}$	$\underset{\left(0\cdot 07\right)}{0\cdot 89}$
CD	$\underset{\left(0\cdot 05\right)}{0\cdot 57}$	$\underset{\left(0\cdot 04\right)}{0\cdot 51}$	$\underset{\left(0\cdot 04\right)}{0\cdot 53}$	$\underset{\left(0\cdot 03\right)}{2\cdot 09}$	$\underset{\left(0\cdot 03\right)}{2\cdot 09}$	$\underset{\left(0\cdot 03\right)}{2\cdot 09}$
AT	$\underset{\left(0\cdot 06\right)}{0\cdot 62}$	$\underset{\left(0\cdot 06\right)}{0\cdot 56}$	$\underset{\left(0\cdot 06\right)}{0\cdot 58}$	$\underset{\left(0\cdot 17\right)}{3\cdot 39}$	$\underset{\left(0\cdot 20\right)}{3\cdot 45}$	$\underset{\left(0\cdot 16\right)}{3\cdot 49}$
POET	$\underset{\left(0\cdot 06\right)}{0\cdot 60}$	$\underset{\left(0\cdot 05\right)}{0\cdot 53}$	$\underset{\left(0\cdot 05\right)}{0\cdot 55}$	$\underset{\left(0\cdot 05\right)}{2\cdot 52}$	$\underset{\left(0\cdot 08\right)}{2\cdot 54}$	$\underset{\left(0\cdot 07\right)}{2\cdot 55}$

Table 4. Table 4 . Simulation Setting 2 with s = 0.1 𝑠 0.1 s=0.1 . Top panel compares k ^ sure superscript ^ 𝑘 sure \widehat{k}^{\text{sure}} with k opt superscript 𝑘 opt k^{\text{opt}} . Bottom two panels compare the proposed approach with POET and AT in terms of ‖ Σ ^ − Σ 0 ‖ 2 / p subscript norm ^ Σ subscript Σ 0 2 𝑝 \|{\widehat{\Sigma}}-{\Sigma_{0}}\|_{2}/p and ‖ Σ ^ − Σ 0 ‖ F / p subscript norm ^ Σ subscript Σ 0 𝐹 𝑝 \|{\widehat{\Sigma}}-{\Sigma_{0}}\|_{F}/p respectively. Standard errors are in parenthesis.

ktr	10			50
p	250	500	1000	250	500	1000
$k^{\text{opt}}$	240	480	950	210	420	830
${\hat{k}}^{\text{sure}}$	240	470	940	210	410	810
CD	$\underset{\left(0\cdot 06\right)}{0\cdot 49}$	$\underset{\left(0\cdot 06\right)}{0\cdot 50}$	$\underset{\left(0\cdot 05\right)}{0\cdot 50}$	$\underset{\left(0\cdot 05\right)}{1\cdot 17}$	$\underset{\left(0\cdot 04\right)}{1\cdot 04}$	$\underset{\left(0\cdot 02\right)}{0\cdot 95}$
AT	$\underset{\left(0\cdot 09\right)}{0\cdot 54}$	$\underset{\left(0\cdot 10\right)}{0\cdot 53}$	$\underset{\left(0\cdot 09\right)}{0\cdot 56}$	$\underset{\left(0\cdot 17\right)}{1\cdot 67}$	$\underset{\left(0\cdot 10\right)}{1\cdot 43}$	$\underset{\left(0\cdot 09\right)}{1\cdot 26}$
POET	$\underset{\left(0\cdot 09\right)}{0\cdot 54}$	$\underset{\left(0\cdot 10\right)}{0\cdot 53}$	$\underset{\left(0\cdot 09\right)}{0\cdot 55}$	$\underset{\left(0\cdot 18\right)}{1\cdot 75}$	$\underset{\left(0\cdot 13\right)}{1\cdot 66}$	$\underset{\left(0\cdot 12\right)}{1\cdot 62}$
CD	$\underset{\left(0\cdot 08\right)}{0\cdot 90}$	$\underset{\left(0\cdot 09\right)}{0\cdot 92}$	$\underset{\left(0\cdot 08\right)}{0\cdot 94}$	$\underset{\left(0\cdot 03\right)}{3\cdot 85}$	$\underset{\left(0\cdot 06\right)}{3\cdot 81}$	$\underset{\left(0\cdot 05\right)}{3\cdot 75}$
AT	$\underset{\left(0\cdot 11\right)}{0\cdot 98}$	$\underset{\left(0\cdot 12\right)}{0\cdot 99}$	$\underset{\left(0\cdot 11\right)}{1\cdot 03}$	$\underset{\left(0\cdot 59\right)}{6\cdot 18}$	$\underset{\left(0\cdot 36\right)}{6\cdot 36}$	$\underset{\left(0\cdot 28\right)}{6\cdot 33}$
POET	$\underset{\left(0\cdot 09\right)}{0\cdot 94}$	$\underset{\left(0\cdot 11\right)}{0\cdot 95}$	$\underset{\left(0\cdot 10\right)}{0\cdot 98}$	$\underset{\left(0\cdot 17\right)}{4\cdot 64}$	$\underset{\left(0\cdot 13\right)}{4\cdot 64}$	$\underset{\left(0\cdot 13\right)}{4\cdot 57}$

Equations

$$\hat{\Sigma }(\varphi ;k)={\varphi }^{\ast }{\Sigma }_w\varphi =\frac{1}{n}{\varphi }^{\ast }(W{W}^{\ast })\varphi ={\varphi }^{\ast }(\varphi {\hat{\Sigma }}_x{\varphi }^{\ast })\varphi .$$

$$\displaystyle\widehat{\Sigma}_{{\mathrm{\scriptscriptstyle CD}}}(k)=\mathbb{E}_{\phi}\big{[}\widehat{\Sigma}_{{\mathrm{\scriptscriptstyle CD}}}(\phi;k)\big{]}.$$

$$\widehat{\Sigma}_{{\mathrm{\scriptscriptstyle CD}}}(k)=\frac{k}{(p^{2}-1)p}\bigg{[}(pk-1){{\widehat{\Sigma}}_{x}}+(p-k)\mbox{Tr}({\widehat{\Sigma}}_{x})I_{p}\bigg{]},$$

$$R(k)=\mathbb{E}\Vert {\hat{\Sigma }}_{\mathrm{CD}}(k)-{\hat{\Sigma }}_x{\Vert }_F^2-\sum _{i,j}^p\text{var}({\hat{\sigma }}_{ij})+2\sum _{i,j}^p\text{cov}({\hat{\sigma }}_{ij}^{(k)},{\hat{\sigma }}_{ij})$$

$$\widehat{\text{mbox}var}({\hat{\sigma }}_{ij})$$

$$\widehat{\text{var}}({\hat{\sigma }}_{ii})$$

$$\widehat{\text{mbox}cov}({\hat{\sigma }}_{ll},{\hat{\sigma }}_{ii})$$

$$\mbox{SURE}(k)=(\eta-1)^{2}\|{\widehat{\Sigma}}\circ{\widehat{\Sigma}}\|_{F}+p{\gamma^{2}}+2{\gamma}(\eta-1){\mbox{Tr}}({\widehat{\Sigma}})\\ +2\bigg{\{}(a_{n}\eta+d_{n}\gamma)\bigg{(}\|{\widetilde{\Sigma}}\circ{\widetilde{\Sigma}}\|_{F}\;-\;\mbox{Tr}({\widetilde{\Sigma}}\circ{\widetilde{\Sigma}})\bigg{)}\\ +(b_{n}\eta+c_{n}\gamma)\bigg{(}{{\mbox{Tr}}({\widetilde{\Sigma}})}^{2}-{\mbox{Tr}}({\widetilde{\Sigma}}\circ{\widetilde{\Sigma}})\bigg{)}+c_{n}(\eta+\gamma){\mbox{Tr}}\bigg{(}{\widetilde{\Sigma}}\circ{\widetilde{\Sigma}}\bigg{)}\bigg{\}}$$

$${\hat{k}}^{\text{sure}}={\text{mbox}argmin}_k\text{mbox}SURE(k).$$

$$\mathbb{E}(({\tilde{\sigma }}_{ij}^s{)}^2)=\frac{n}{n-1}{\sigma }_{ij}^2+\frac{{\sigma }_{ii}{\sigma }_{jj}}{n-1}$$

$$\mathbb{E}({\tilde{\sigma }}_{ii}^s{\tilde{\sigma }}_{jj}^s)=\frac{n+1}{n-1}{\sigma }_{ii}{\sigma }_{jj}+\frac{2(n+2)}{n(n-1)}{\sigma }_{ij}^2.$$

$$\displaystyle{\mathbb{E}}\bigg{[}\frac{n(n^{2}-1)}{n^{3}+n^{2}-2n-4}({\widetilde{\sigma}}_{ij}^{s})^{2}-\frac{n(n-1)}{n^{3}+n^{2}-2n-4}{\widetilde{\sigma}}_{ii}^{s}{\widetilde{\sigma}}_{jj}^{s}\bigg{]}={\sigma_{ij}^{2}}.$$

$$\displaystyle{\mathbb{E}}\bigg{[}\frac{2(n-1)(n+2)}{2n+4-n^{3}-n^{2}}({\widetilde{\sigma}}_{ij}^{s})^{2}-\frac{n^{2}(n-1)}{2n+4-n^{3}-n^{2}}{\widetilde{\sigma}}_{ii}^{s}{\widetilde{\sigma}}_{jj}^{s}\bigg{]}={\sigma_{ii}}{\sigma_{jj}}.$$

$$\widehat{\text{mbox}Var}({\tilde{\sigma }}_{ij}^s)=\frac{{\hat{\sigma }}_{ij}^2+\widehat{{\sigma }_{ii}{\sigma }_{jj}}}{n-1}.$$

$$\text{mbox}Cov({\tilde{\sigma }}_{jj}^s,{\tilde{\sigma }}_{ii}^s)=\frac{2(n+2)}{n(n-1)}{\sigma }_{ij}^2+\frac{2}{n-1}{\sigma }_{ii}{\sigma }_{jj}.$$

$$\Vert {\hat{\Sigma }}_{\mathrm{CD}}(k)-\hat{\Sigma }{\Vert }_F^2$$

$$\displaystyle={\displaystyle\sum\limits_{i\neq j}^{p}}{(\eta-1)^{2}{\widehat{\sigma}}_{ij}^{2}}+\displaystyle\sum\limits_{i=j}^{p}\bigg{[}(\eta-1){\widehat{\sigma}}_{ii}+\gamma\bigg{]}^{2}$$

$$=(\eta -1{)}^2\Vert \hat{\Sigma }\circ \hat{\Sigma }{\Vert }_F+p{\gamma }^2+2\gamma (\eta -1)\text{mbox}Tr(\hat{\Sigma }),$$

$$\displaystyle={\displaystyle\sum\limits_{i\neq j}^{p}}{\eta}{\mbox{Var}}({\widehat{\sigma}}_{ij})+{\displaystyle\sum\limits_{i=1}^{p}}\bigg{\{}\eta{\mbox{Var}}({\widehat{\sigma}}_{ii})+\gamma\;{\mbox{cov}}\bigg{(}{\displaystyle\sum\limits_{\underset{l\neq i}{l=1}}^{p}}{\widehat{\sigma}}_{ll}+{\widehat{\sigma}}_{ii},{\widehat{\sigma}}_{ii}\bigg{)}\bigg{\}}$$

$$\displaystyle={\displaystyle\sum\limits_{i\neq j}^{p}}{\eta}{\mbox{Var}}({\widehat{\sigma}}_{ij})+{\displaystyle\sum\limits_{i=1}^{p}}\bigg{\{}(\eta+\gamma){\mbox{Var}}({\widehat{\sigma}}_{ii})+\gamma{\displaystyle\sum\limits_{\underset{l\neq i}{l=1}}^{p}}{\mbox{cov}}\bigg{(}{\widehat{\sigma}}_{ll},{\widehat{\sigma}}_{ii}\bigg{)}\bigg{\}}.$$

$$\widehat{\text{optimism}}$$

$$\displaystyle=2\bigg{\{}(a_{n}\eta+d_{n}\gamma)\bigg{(}\|{\widetilde{\Sigma}}\circ{\widetilde{\Sigma}}\|_{F}\;-\;\mbox{Tr}({\widetilde{\Sigma}}\circ{\widetilde{\Sigma}})\bigg{)}+(b_{n}\eta+c_{n}\gamma)$$

$$\displaystyle\bigg{(}{{\mbox{Tr}}({\widetilde{\Sigma}})}^{2}-{\mbox{Tr}}({\widetilde{\Sigma}}\circ{\widetilde{\Sigma}})\bigg{)}+c_{n}(\eta+\gamma){\mbox{Tr}}\bigg{(}{\widetilde{\Sigma}}\circ{\widetilde{\Sigma}}\bigg{)}\bigg{\}},$$

Full text

Keywords: Compressed sensing, Dimension Reduction, Low-rank, Factor model, Spiked covariance models, SURE

1. Introduction

Estimating a covariance matrix based on a sample of multivariate observations is a classical problem in statistics with applications across multitude of scientific disciplines. Time series analysis (Chen et al., 2013; Basu et al., 2015), portfolio optimization (Fan et al., 2008; Bai et al., 2011), gene networks (Butte et al., 2000; Schäfer et al., 2005), climate studies (Houtekamer and Mitchell, 2001; Hamill et al., 2001; Furrer et al., 2006), spatial data analysis (Kaufman et al., 2008), longitudinal data analysis (Smith and Kohn, 2002; Wu and Pourahmadi, 2003) among many other disciplines rely critically on the knowledge of the covariance structure. Recent efforts have focussed on high-dimensional data, where the dimension $p$ can be much larger than the sample size $n$ (Pourahmadi, 2011; Cai et al., 2016; Fan et al., 2016). An unstructured $p\times p$ covariance matrix has $O(p^{2})$ free entries. In moderate to high-dimensional situations, a general idea is to assume a restricted parameter space with much lower effective degrees of freedom. Examples include tapered covariance matrices (Furrer and Bengtsson, 2007), bandable covariance matrices (Wu and Pourahmadi, 2003; Bickel and Levina, 2008b), Toeplitz covariance matrices (Wu and Pourahmadi, 2009; McMurry and Politis, 2010; Xiao et al., 2012), sparse covariance matrices (Bickel and Levina, 2008a; Karoui, 2008a, b; Rothman et al., 2009; Cai and Liu, 2011), spiked sparse covariance matrices (Johnstone, 2001; Ma et al., 2013; Cai et al., 2015), covariances with a kronecker product structure (Werner et al., 2008), penalized likelihood estimation (Huang et al., 2006; d’Aspremont et al., 2008; Lam and Fan, 2009; Ravikumar et al., 2011) and regularization of principal components (Zou et al., 2006; Hoff, 2009; Johnstone and Lu, 2012; Cai et al., 2013). The challenge lies in the selection of appropriate structure or method for a specific data domain application.

An alternative approach to controlling the complexity in covariance matrix estimation is through low intrinsic dimensionality. Low intrinsic dimensionality posits that the dependencies between the variables are captured by a small number of latent components, also called factors, explicitly seperating the common variation from variable-specific noise in the observed variables. Factor models in the high-dimensional regime have been used in a myriad of applications in economics and finance (Engle and Watson, 1981; Goldfarb and Iyengar, 2003). There are two disparate literatures which express the covariance matrix as the sum of a low-rank and sparse matrix: (a) factor models (Bhattacharya and Dunson, 2011; Fan et al., 2013; Pati et al., 2014), and (b) spiked covariance model (Johnstone, 2001). It is worthwhile to point out that the literature on both factor models and spiked covariance models assume the low-rank matrix is also sparse.

While the sparsity assumptions on the covariance matrix or the low rank component are well motivated for specific applications (Zou et al., 2006; Shen and Huang, 2008; Witten et al., 2009; Bhattacharya and Dunson, 2011), recent studies reveal striking correlations between several genes/loci or gene networks and the features of the diseases they cause (Jimenez-Sanchez et al., 2001; Ideker and Sharan, 2008). The correlation between the attributes of complex disease genes is more extensive and stronger than previously thought. There is evidence that the etiology of many complex diseases involves, rather than a few genes/loci with large effects, many genes, each of which contributes a small risk, interacting with each other or with environmental risk factors to cause these complex diseases. These small effects are organized in networks/pathways that have distinct features.

In this article, we attempt to mitigate the aforementioned gaps in the literature by (i) providing an efficient way of estimating covariance matrices which admit a decomposition of the form that is written as a sum of low-rank plus sparse matrix albeit with dense to moderately sparse low-rank structures, (ii) providing a concrete principled way of choosing the dimension of the low-rank matrix using Stein’s Unbiased Risk Estimation (SURE) theory. We show via simulations that our approach is readily scalable to massive covariance matrices. The approach is based on a “compression-decompression” (C-D) mechanism. The C-D mechanism proceeds by projecting the high dimensional observations to a lower dimension to form compressed measurements and then decompress them back to the original dimension. The key idea is to use the sample covariance matrix of the compressed-decompressed data as an estimator of $\Sigma$ instead of the sample covariance matrix ${\widehat{\Sigma}}$. The covariance estimation problem considered in this paper is related to the covariance sketching problem considered in Dasarathy et al. (2015). In covariance sketching, the goal is to estimate the covariance matrix of high dimensional random vectors based on the low-dimensional projections where the same instance of projection matrix is used across all observations. The choice of the projection matrix, in our setting, varies across observations and is integrated out to obtain the final estimator.

Our simulation studies demonstrate the efficacy of our approach and show the C-D estimator to be highly competitive to the recent state-of-the-art covariance estimators such as POET (Fan et al., 2013) and adaptive thresholding (Cai and Liu, 2011).

2. Compression - Decompression Covariance Estimator

We now elucidate the compression-decompression (C-D) estimator. Let $X=[x_{1},\ldots,x_{n}]$ denote a $p\times n$ data matrix with the columns $x_{i}$ independent and identically distributed from a $p$-variate distribution whose covariance we wish to estimate. We assume the data to be column-centered, so that $\mathbb{E}(x_{i})=0$ for all $i$, and thereby define $\widehat{\Sigma}_{x}=n^{-1}\sum_{i=1}^{n}x_{i}x_{i}^{{\mathrm{\scriptscriptstyle T}}}=XX^{{\mathrm{\scriptscriptstyle T}}}/n$ as the sample covariance matrix. The mechanism proceeds by projecting the data to a lower-dimensional space to form compressed measurements, computing the sample covariance in the compressed space, and lifting back to the ambient space via a decompression operation.

Specifically, given $k<p$ and a $k\times p$ unitary matrix $\phi$ (with $\phi\phi^{*}=I_{k}$; $\phi^{*}$ denoting the complex conjugate of $\phi$), project the data $x_{i}$ from ${\mathbb{R}}^{p}\mapsto{\mathbb{R}}^{k}$ to create compressed data $w_{i}=\phi x_{i}$; let $W=[w_{1},\ldots,w_{n}]$ denote the corresponding $k\times n$ matrix. The $k\times k$ sample covariance matrix $\Sigma_{w}=n^{-1}\sum_{i=1}^{n}w_{i}w_{i}^{*}=WW^{*}/n$ is expected to be more stable compared to the $p\times p$ matrix $\widehat{\Sigma}_{x}$. To obtain a $p\times p$ covariance estimate for our original problem, we decompress $\Sigma_{w}$ using the transformation $\phi^{*}$ to define

[TABLE]

Observe that the regularization in this framework arises from the compression operation, which is entirely different from assuming $\ell_{q}$ type sparsity on the covariance matrix or its various decompositions. One of the motivations behind this approach comes from Theorem 8.1 of Pati et al. (2014), where it is shown that if the true $\Sigma_{0}=\Lambda_{0}{\Lambda_{0}}^{{\mathrm{\scriptscriptstyle T}}}+\sigma^{2}I_{p}$ with $\Lambda_{0}\in{\mathbb{R}}^{p\times k}$, and one chooses $\phi=\Lambda_{0}$, then ${\widehat{\Sigma}}$ concentrates around $\Sigma_{0}$ in operator norm with high probability when $k\ll p$ and $\sigma^{2}$ is bounded. It is important to mention here that a similar concentration for the sample covariance matrix in a spiked covariance model requires $\Sigma$ to be have a low effective rank (Bunea and Xiao, 2015), necessitating $\sigma^{2}=O(1/p)$, which is fairly restrictive.

Evidently, the estimator in (1) depends on two unknown parameters, the compression matrix $\phi$ and the projected dimension $k$. There has been recent work in the regression context (Guhaniyogi and Dunson, 2013) where very high-dimensional covariates are projected to a lower-dimensional space using one particular instance of a random sensing matrix. However, this approach of fixing $\phi$ has very poor performance in our setting. An alternative way may be to estimate $\phi$ from the data is computationally intensive as it requires estimating $p\times k$ many parameters. Here, instead of trying to estimate the high-dimensional parameter $\phi$, we average over the ensemble of unitary matrices which can be performed in closed-form using random matrix results (Marzetta et al., 2011). Specifically, let $h(\cdot)$ denote the Haar measure on the space of $k\times p$ unitary matrices satisfying $h(\phi\Psi)=h(\phi)$ for all non-stochastic $p\times p$ unitary matrices $\Psi$. Letting $\mathbb{E}_{\phi}$ denote the expectation with respect to $h$, define the C-D estimator

[TABLE]

Based on recent random matrix techniques as in Marzetta et al. (2011), the above expectation can be computed as

[TABLE]

where $\mbox{Tr}(A)$ denotes the trace of a matrix $A$. Clearly, averaging over the Haar distribution introduces an appropriate shrinkage on the sample covariance matrix resulting the estimator full rank even if $p\gg n$. As noted by Marzetta et al. (2011), (2) bears resemblance with shrinkage estimators of the type $a\widehat{\Sigma}_{x}+(1-a)I_{p}$, where $0<a<1$ and $I_{p}$ is a $p\times p$ identity matrix, originally proposed by Ledoit and Wolf (2004). However, a key difference is that the effect of shrinkage $a\widehat{\Sigma}_{x}$ is compensated by $\mbox{Tr}(\widehat{\Sigma}_{x})I_{p}$ in (2) instead of just $I_{p}$. This helps preserving the largest eigenvalues of the resulting estimator.

A fundamental principle of statistical decision theory is that there exists an interior optimum in the trade-off between bias and estimation error. One way of attaining this optimal trade-off is simply to take a properly weighted average of the biased and unbiased estimators. Refer, for example, to the seminal work by Stein et al. (1956), who showed that shrinking sample means towards a constant can, under certain circumstances, improve accuracy. The crux of our method is to shrink the unbiased but very variable sample covariance matrix towards the biased but less variable identity covariance matrix and to thereby obtain a more efficient estimator. In addition, the resulting estimator, in (2) is invertible and well-conditioned, which is of crucial importance in settings where one needs to estimate the inverse, for example, in portfolio selection (Ledoit et al., 2003), Gaussian graphical models (Meinshausen and Bühlmann, 2006; Yuan and Lin, 2007) among others.

3. Choice of the compressed dimension

In this section, we propose a data-driven method for choosing the size of the compressed dimension. The main idea is to consider the Frobenius risk of the estimator for each choice of compressed dimension, $k$, and then choose that value which minimizes the Frobenius risk. Since the risk depends on the truth, we derive an unbiased estimator of the Frobenius risk curve and find the minimizer of the estimated risk curve to choose a value of $k$. We use Stein’s Unbiased Risk Estimation (SURE) theory to find an unbiased estimator of the Frobenius risk associated with the C-D estimator.

SURE was originally proposed in Stein (1981) in deriving an unbiased estimator of the risk of James–Stein estimate. Efron (1986, 2004) applied SURE to prediction problems which was called covariance penalty method. Li and Zou (2016) proved the asymptotic properties of SURE information criterion for large bandable covariance matrices and proposed a family of generalized SURE ($\mbox{SURE}_{c}$) indexed by $c$ for covariance matrix estimation, where $c$ is some constant. For bandable covariance matrices, Li and Zou (2014) claimed that $\mbox{SURE}_{2}$ and $\mbox{SURE}_{\log{}{n}}$ can be regarded as AIC and BIC analogues, respectively, for covariance matrix estimation problems. Xiao and Bunea (2014) proposed an improved version of the banding estimator obtained in Bickel and Levina (2008b) and used SURE-type approach for selecting the bandwidth for the banding estimator.

Let $R(k)=\mathbb{E}\|\widehat{\Sigma}_{{\mathrm{\scriptscriptstyle CD}}}(k)-{\Sigma}_{0}\|_{F}^{2}$ be the Frobenius risk associated with the proposed estimator for a fixed $k$. We use the following risk identity proved in Yi and Zou (2013),

[TABLE]

where ${\widehat{\Sigma}}=[{\widehat{\sigma}}_{ij}]$ is the usual sample covariance matrix, $\widetilde{\Sigma}=[{\widetilde{\sigma}}_{ij}]$ is the maximum likelihood estimator of the covariance matrix, with ${\widehat{\Sigma}}=\frac{n}{n-1}{\widetilde{\Sigma}}$, and $\widehat{\Sigma}_{{\mathrm{\scriptscriptstyle CD}}}(k)=[{\widehat{\sigma}}_{ij}^{(k)}]$ is the proposed covariance estimator (2). The third term on the right hand side of (3) is referred to as the optimism (Efron, 2004). The second term on the right hand is the same for all estimators of $\Sigma_{0}$.

Standard results from multivariate statistics (Anderson, 1984) imply ${\text{var}}({\widehat{\sigma}}_{ij})=\frac{\sigma_{ij}^{2}+\sigma_{ii}\sigma_{jj}}{n-1},{\text{var}}({\widehat{\sigma}}_{ii})=\frac{2\sigma_{ii}^{2}}{n-1}$ and {\mbox{cov}}({\widehat{\sigma}}_{ll},{\widehat{\sigma}}_{ii})=\frac{2(n+2)}{n(n-1)}{\sigma_{il}}^{2}+\big{(}\frac{n+1}{n-1}-1\big{)}{\sigma_{ii}}{\sigma_{ll}}. Lemma 3.1 provides unbiased estimators to each of these quantities which enables us to provide an unbiased estimator of $R(k)$. Recall ${\widehat{\Sigma}}$ and ${\widetilde{\Sigma}}$ defined in the previous paragraph.

Lemma 3.1**.**

*Let ${\widehat{\cdot}}$ denote an unbiased estimator of the quantity of interest. Then *

[TABLE]

We now derive the SURE criterion for the Frobenius risk $R(k)$ defined in (3).

Theorem 3.2**.**

An unbiased estimator of the Frobenius risk associated with $\widehat{\Sigma}_{{\mathrm{\scriptscriptstyle CD}}}(k)$ is given by

[TABLE]

where $\eta=\frac{k(pk-1)}{p(p^{2}-1)}$, $\gamma=\frac{p-k}{p(p^{2}-1)}$, and $\Sigma_{1}\circ\Sigma_{2}$ denotes the Schur product between two matrices $\Sigma_{1}$ and $\Sigma_{2}$ of the same dimension.

Having obtained the SURE criterion, we can now minimize it with respect to $k$ to select the compressed dimension. Specifically, define

[TABLE]

Our simulation results show below that ${\widehat{k}^{\text{sure}}}$ provides an accurate estimate of the intrinsic dimensionality $k$ in the examples considered.

4. Experiments on Synthetic Data

In this section, we consider a number of simulation cases to compare our proposed approach in terms of (a) accuracy of the SURE method in estimating $\widehat{k}^{\text{sure}}$, (b) accuracy of covariance matrix estimation in operator norm, and (c) accuracy of covariance matrix estimation in frobenius norm. The Frobenius norm ($\|\cdot\|_{F}$) and the operator norm ($\|\cdot\|_{2}$) are defined in the usual way with $\|A\|_{F}=\sqrt{\mbox{trace}(A^{{\mathrm{\scriptscriptstyle T}}}A)}$ and $\|A\|_{2}=s_{{max}}(A)$ where $s_{{max}}(A)$ denotes the largest singular value of $A$.

We compare our method with Principal Orthognoal complement Thresholding (POET) of Fan et al. (2013) which is based on an additive decomposition of the covariance matrix in terms of a low rank matrix and a sparse residual covariance matrix. POET estimates the factors and the loadings by thresholding the principal components of the sample covariance matrix. We also compare with adaptive thresholding (AT) of Cai and Liu (2011) which thresholds the entries of the sample covariance matrix, with the resulting thresholded estimator ${\widehat{\Sigma}}$ being of the form ${\widehat{\Sigma}}_{jj^{\prime}}=S_{jj^{\prime}}1(|S_{jj^{\prime}}|>{\delta_{\kappa_{jj^{\prime}}}})$, where $\delta$ is the tuning parameter and $\kappa_{jj^{\prime}}$ is a threshold specific to the corresponding entry of $S$. We choose the tuning parameter $\delta$ by 5-fold cross-validation as suggested by Cai and Liu (2011).

The two simulation settings considered here are described below:

1. (1)

$y_{i}$, $i=1,\ldots,n$ are generated from $N_{p}(0,\Sigma_{0})$, where $\Sigma_{0}=\Lambda_{0}{\Lambda_{0}}^{{\mathrm{\scriptscriptstyle T}}}+\sigma_{0}^{2}{\mathrm{I}}_{p}$ and $\Lambda$ is a $p\times{\text{ktr}}$ matrix with $(1-s)\times 100\%$ non-zero entries. Here, $s\in(0,1)$, is the sparsity parameter. We choose different values of $s$ that lead to moderately sparse to dense covariance matrices. The nonzeros entries are independently drawn from a standard normal distribution. 2. (2)

This setting is designed to illustrate the performance of our approach under model misspecification. We let $\Sigma=\Lambda_{0}\Lambda_{0}^{{\mathrm{\scriptscriptstyle T}}}+\Omega_{0}$, where $\Lambda_{0}$ is as in simulation setting (1), but $\Omega_{0}$ is nondiagonal corresponding to the covariance matrix of an autoregressive sequence with pure error variance 0.4 and autoregressive coefficient 0.1.

For each simulation setting, we choose the sample size $n=100$, dimension $p=250,500,1000$ and the true number of factors ktr = 10,50. For each $(n,p,\text{ktr})$ triplet, we consider 100 simulation replicates. We consider two different values of the sparsity parameter $s$, $s\in\{0.1,0.5\}$. $s=0.1$ randomly sets $10\%$ entries in the factor loadings matrix to 0. This corresponds to the extreme situation in which there is not a single sparse entry in the true covariance matrix. Similarly, $s=0.5$ randomly sets $50\%$ entries in the factor loadings matrix to 0 and corresponds to the moderately sparse regime.

To evaluate the accuracy of $\widehat{k}^{\text{sure}}$ , we compare with $k^{\text{opt}}$, defined as the minimizer of the true predictive risk function which assumes knowledge of the truth. For different settings, results across simulation replicates are summarized for two two simulation settings to compare the matrix norm differences between the estimator resulting from different methods and the truth. In particular, normalized average operator norm error ($\|\cdot\|_{2}/p$), in the top panel, and normalized average frobenius norm error ($\|\cdot\|_{F}/p$), in the bottom panel, across 100 replicates is provided with standard error in paranthesis. The tables indicate that the SURE method is very accurate in estimating $k^{\text{opt}}$. From Tables 1, 2 it becomes evident as the number of model parameters increases and the sparsity reduces, the performance of AT and POET deteriorates, in terms of both operator and frobenius norms, due to the sparsity assumption on both estimators, while the C-D estimator has more robust performance. Even in Tables 3, 4, where the truth is misspecified for both C-D and AT, and in fact designed to favor POET, C-D performs better than its competitors. We also evaluate the performance of CD estimator across various levels of sparsity on the low-rank structure for fixed sample size $n$ and dimension $p$. Figure 1 displays the superior performance of CD estimator in dense ($s=0.1$) to moderately sparse ($s=0.7$) regimes.

5. Discussion

In this article, we developed a simple but useful method for estimating covariance matrices with dense low-rank structures under the assumption of low intrinsic dimensionality. We also provide a principled framework for choosing the size of the low- rank dimension using SURE theory. We observe excellent performances of the proposed method in terms of scalability to high dimensions and capability of dealing with model misspecification. From Figure 1, CD estimator outperforms AT and POET in dense $(s=0.1)$ to moderately sparse $(s=0.7$) regimes.