Robust estimator of the correlation matrix with sparse Kronecker structure for a high-dimensional matrix-variate

Abstract

It is of interest in many applications to estimate the correlation matrix of a high dimensional matrix-variate $\mathbf{X}\in {\mathbb{R}}^{p\times q}$. Existing works usually impose strong assumptions on the distribution of $\mathbf{X}$ such as sub-Gaussian or strong moment conditions. These assumptions may be violated easily in practice and a robust estimator is desired. In this paper, we consider the case where the correlation matrix has a sparse Kronecker structure and propose a robust estimator based on Kendall’s $\tau$ correlation. The proposed estimator is extended further to tensor data. The theoretical properties of the estimator are established, showing that Kronecker structure actually increases the effective sample size and leads to a fast convergence rate. Finally, we apply the proposed estimator to bigraphical model, obtaining an estimator of better convergence rate than the existing results. Simulations and real data analysis confirm the competitiveness of the proposed method.

Introduction

With the rapid development of data-collection techniques, matrix-valued data are often encountered in biostatistics and medical research [29], such as Nuclear Magnetic Resonance (NMR) data [35], electroencephalograph (EEG) data [27] and functional Magnetic Resonance Imaging (fMRI) data [6], etc. For example, in imaging genetics, each subject has a weighted (or binary) adjacency matrix for characterizing structural (or functional) connectivity pattern, which can be jointly analyzed with genetic data for discovering genes of many neuropsychiatric and neurological disorders [14], [22], [25], [31]. Moreover, on study of an electroencephalography dataset of alcoholism, each subject was exposed to a stimulus. Voltage values were measured from 64 channels of electrodes placed on the subject’s scalp for 256 time points, which forms a 256 × 64 matrix. It is of scientific interest to study the association between alcoholism and the pattern of voltage over times and channels [13], [17], [39]. Suppose the matrix-variates ${\mathbf{X}}_k\in {\mathbb{R}}^{p\times q},k\in \{1,\dots ,n\}$ are independent copies of $\mathbf{X}\in {\mathbb{R}}^{p\times q}$. For each $k$, denote ${\mathbf{X}}_k={\left(X_{st,k}\right)}_{1\le s\le p,1\le t\le q}$ and $\mathrm{vec}\left({\mathbf{X}}_k\right)\in {\mathbb{R}}^{pq}$ to be the vector obtained by stacking the columns of ${\mathbf{X}}_k$.

Since the covariance matrix plays a critical role in statistic inference, estimating the covariance matrix of a matrix-variate has attracted much attention. Denote the covariance matrix of $\mathrm{vec}\left({\mathbf{X}}_k\right)$ as

$$\text{cov}\left(\mathrm{vec}\left({\mathbf{X}}_k\right)\right)=\mathbf{\Gamma }\in {\mathbb{R}}^{pq\times pq}$$

When both $p$ and $q$ are fixed, many works have been developed on estimating $\mathbf{\Gamma }$ [5], [7]. When $p$ and $q$ are diverging, estimating $\mathbf{\Gamma }$ can be challenging. Firstly, the number of the parameters in model (1) is $pq(pq+1)∕2$ by the symmetry of $\mathbf{\Gamma }$, which can be much larger than the sample size $n$. This will greatly reduce the computational efficiency and statistical accuracy. For example, when performing face recognition, Li and Yuan in [19] analyzed merely 400 observations with each of size 112 × 92, leading to more than 10,000 unknown parameters. Secondly, model (1) simply stacks a matrix into a vector, losing information brought by matrix structure. To overcome these difficulties, it is common to impose a Kronecker structure on $\mathbf{\Gamma }$ in the literature [11], [15], [16], [23], [32], [37], etc. For example, Teng and Huang in [30] proposed a Kronecker product matrix to model gene-experiment interactions. This gene expression matrix measured over multiple tissues is transposable, implying that potentially both the rows and/or columns are correlated. Similar idea has also been used in spatio-temporal data. Specifically,

$$\mathbf{\Gamma }=\mathbf{\Psi }\otimes \mathbf{\Sigma },$$

where $\mathbf{\Psi }\in {\mathbb{R}}^{q\times q}$ represents the dependence among the columns of ${\mathbf{X}}_k$, while $\mathbf{\Sigma }\in {\mathbb{R}}^{p\times p}$ denotes that of the rows. Both $\mathbf{\Psi }$ and $\mathbf{\Sigma }$ are symmetric matrices. Model (2) not only retains the dependence patterns among columns and rows, but also reduces the number of unknown parameters. Obviously, model (2) has $(p^2+q^2+p+q)∕2-1$ free parameters, much smaller than that of model (1). Many works have been developed for model (2) [15], [16], [23], [24], [32], [33], [37]. For example, under additional sparsity assumptions on $\mathbf{\Psi }$ and $\mathbf{\Sigma }$, Leng and Pan in [15] proposed their estimators, denoted as ${\hat{\mathbf{\Psi }}}_{LP}$ and ${\hat{\mathbf{\Sigma }}}_{LP}$, taking the advantages of both Kronecker structure and sparsity. It is shown that ${\hat{\mathbf{\Psi }}}_{LP}$ and ${\hat{\mathbf{\Sigma }}}_{LP}$ have the convergence rates of order $O_p\left(\sqrt{logq∕\left(np\right)}\right)$ and $O_p\left(\sqrt{logp∕\left(nq\right)}\right)$, respectively, in terms of the element-wise ${\ell }_{\infty }$ norm.

However, there are some potential limitations on estimating the covariance matrix. As the covariance matrix is not scale invariant, researchers always standardize the variables to have zero mean and unit variance ahead. In this way, the estimation of covariance matrix is turned out to be that of correlation matrix. Without this prior information, more unknown parameters should be estimated and the convergence rate will be slowed down. Besides, correlation matrix is sufficient to reveal the dependence among variables without estimating the marginal variance. More discussion on the importance of high-dimensional correlation matrix estimation is referred to [4]. Thus, we focus on the estimation of the correlation matrix ${\mathbf{R}}^{\Gamma }\in {\mathbb{R}}^{pq\times pq}$ of $\mathrm{vec}\left({\mathbf{X}}_k\right)$ under the following Kronecker structure

$${\mathbf{R}}^{\Gamma }={\mathbf{R}}^{\Psi }\otimes {\mathbf{R}}^{\Sigma },$$

where ${\mathbf{R}}^{\Psi }\in {\mathbb{R}}^{q\times q}$ and ${\mathbf{R}}^{\Sigma }\in {\mathbb{R}}^{p\times p}$ are both symmetric matrices. Since the diagonal elements of ${\mathbf{R}}^{\Gamma }$ are all equal to one, it is easy to see that $\text{diag}\left({\mathbf{R}}^{\Psi }\right)=c{\mathbf{1}}_q$ and $\text{diag}\left({\mathbf{R}}^{\Sigma }\right)=c^{-1}{\mathbf{1}}_p$ for some $c>0$, where for any integer $m$, ${\mathbf{1}}_m$ denotes the vector of length $m$ with all elements being one. Note that multiplying ${\mathbf{R}}^{\Psi }$ by $\tilde{c}$ and ${\mathbf{R}}^{\Sigma }$ by ${\tilde{c}}^{-1}$ for any $\tilde{c}>0$, one gets the same ${\mathbf{R}}^{\Gamma }$. For identifiability, we set all the diagonal elements of both ${\mathbf{R}}^{\Psi }$ and ${\mathbf{R}}^{\Sigma }$ being one. That is, ${\mathbf{R}}^{\Psi }$ and ${\mathbf{R}}^{\Sigma }$ denote the correlations among the columns and rows of ${\mathbf{X}}_k$, respectively. In this paper, we focus on the case of $p$ and $q$ being much larger than $n$. Similar to [15], we assume that both ${\mathbf{R}}^{\Psi }$ and ${\mathbf{R}}^{\Sigma }$ are sparse matrices.

The estimator of covariance matrix can deliver an estimator of the correlation matrix in a way similar to [15]. However, most of the works mentioned above are developed under Gaussian, sub-Gaussian assumption or strong moment conditions on ${\mathbf{X}}_k$. As mentioned in [8], there are several disadvantages of these assumptions. Firstly, these estimators are not robust to outliers or heavy tailed distributions. Secondly, the Gaussian or sub-Gaussian assumption may not be realistic for many real-world applications, such as equities and market indices data [8], and gene expression data analysis. For example, when analyzing the atlas of gene expression in the mouse aging project dataset [38], Ning and Liu in [23] found out by the quantile–quantile plot that many gene expression levels may not be normally distributed. Therefore, it is desirable to develop a robust estimation under weaker assumptions on the distribution. Particularly, Han et al. in [20] and [8], [9] proposed a transelliptical family of distribution for random vector $\mathbf{y}={(y_1,\dots ,y_p)}^{\top }$. $\mathbf{y}$ follows a transelliptical distribution, denoted by $\mathbf{y}\sim T{E}_p({\Sigma }^0,\xi ;f_1,\dots ,f_p)$, if there exist some unknown univariate strictly increasing functions $(f_1,\dots ,f_p)$ such that ${(f_1\left(y_1\right),\dots ,f_p\left(y_p\right))}^{\top }$ follows an elliptical distribution $E{C}_p(0,{\Sigma }^0,\xi )$ with the location parameter zero and the scale parameter ${\Sigma }^0$, whose diagonal elements are one. ${\Sigma }^0$ is called latent generalized correlation matrix [8]. Particularly, if ${(f_1\left(y_1\right),\dots ,f_p\left(y_p\right))}^{\top }$ follows a normal distribution with correlation matrix ${\Sigma }^0$, then ${\Sigma }^0$ is called latent correlation matrix. For brevity, a unified term “latent correlation matrix” is used to denote “latent (generalized) correlation matrix” throughout this paper. Since each $f_j,j\in \{1,\dots ,p\}$ is a strictly increasing function, its inverse function exists. Consequently, viewing the nuisance parameters $(f_1,\dots ,f_p)$ as contaminations, $\mathbf{y}$ can be viewed as contaminated observation of some elliptical or normal variable with correlation matrix ${\Sigma }^0$, which is the parameter of interest.

In this paper, we apply the frame of [8] to matrix-variates, assuming that $\mathrm{vec}\left({\mathbf{X}}_k\right),k\in \{1,\dots ,n\}$ are independently and identically distributed $(i.i.d.)$ variables from transelliptical family $T{E}_p({\mathbf{R}}^{\Gamma },\xi ;f_1,\dots ,f_p)$, where the latent correlation matrix ${\mathbf{R}}^{\Gamma }$ has the sparse Kronecker structure as in (3). A robust estimator of the latent correlation matrix is developed based on Kendall’s $\tau$ correlation. We extend our approach further to the tensor type data and apply the estimator to bigraphical model.

Our estimator is an extension of that of [15], relaxing the assumptions on ${\mathbf{X}}_k$. However, there are several major differences between two works. First, we focus on correlation matrix while the work of Leng and Pan [15] focuses on covariance matrix. The correlation matrix has a special property shown in (4) in Section 2, which is not owned by covariance matrix. As a result, the procedures of our estimators, utilizing this property, are significantly different from that of [15]. Second, our estimator is robust with weaker assumption on the distribution of the data. Third, because the Kendall’s $\tau$ estimator involves a high dimensional binary random vector, the techniques used for establishing the theoretical properties are different from that of [15]. Interestingly, the convergence rates obtained are similar to those of [15]. In addition, our work improves some existing works. For example, Niu and Zhao in [24] proposed a robust estimator for a general model, which included (3) as a special case. However, they did not assume the sparsity of the parameters, and the corresponding estimator has a slower convergence rate of the order $O_p\left(\sqrt{(p^2+q^2+logn)∕n+(pqlog\left(pq\right))∕n}\right)$. Our estimator improves those of [24] by utilizing the sparsity of the parameters. In addition, our estimator improves the results of [23] on the bigraphical model. Ning and Liu in [23] proposed an estimator for the nonparanormal bigraphical model based on Kendall’s $\tau$ correlation. However, they failed to utilize the Kronecker structure when estimating ${\mathbf{R}}^{\Gamma }$, obtaining only a slower convergence rate than ours. Details are referred to Section 3.

The main contents of the article are organized as follows. In Section 2, we propose the robust estimator for the latent correlation matrix with sparse Kronecker structure (3) based on Kendall’s $\tau$ correlation. The theoretical properties of the proposed estimators are established in Section 3.1. Then, our estimators are applied to bigraphical models in Section 3.2. The generalization to tensor data is presented in Section 4. Simulation results and real data analysis are presented in Section 5 and Section 6, respectively. All the proofs are presented in Appendix.

Finally, we introduce some notations. For any constant $a\in \mathbb{R}$, ${\left(a\right)}_{+}=max\{a,0\}$. For any vector $\mathbf{v}\in {\mathbb{R}}^m$, $\Vert \mathbf{v}\Vert$ denotes the Euclidean norm of $\mathbf{v}$. For any $m\times m$ matrix $\mathbf{M}=\left(M_{ij}\right)$, ${\Vert \mathbf{M}\Vert }_2$ denotes the operator norm and ${\Vert \mathbf{M}\Vert }_F$ denotes the Frobenius norm of $\mathbf{M}$. ${\Vert \mathbf{M}\Vert }_{max}$ stands for the element-wise ${\ell }_{\infty }$ norm, i.e. ${\Vert \mathbf{M}\Vert }_{max}={max}_{i,j}\left|M_{ij}\right|$. $\text{tr}\left(\mathbf{M}\right)$ denotes the trace of $\mathbf{M}$. $\mathbf{M}⪰0$ indicates that $\mathbf{M}$ is a positive semidefinite matrix, while $\mathbf{M}\succ 0$ denotes a positive definite matrix. For any symmetric matrix $\mathbf{M}$, ${\omega }_{min}\left(\mathbf{M}\right)$ and ${\omega }_{max}\left(\mathbf{M}\right)$ denote the minimum and maximum eigenvalues of $\mathbf{M}$, respectively. For clarity, for any random vector $\mathbf{y}={(y_1,\dots ,y_p)}^{\top }\in {\mathbb{R}}^p$, the Pearson correlation matrix and Kendall’s $\tau$ correlation matrix of $\mathbf{y}$ are denoted as $\mathrm{corr}\left(\mathbf{y}\right)\in {\mathbb{R}}^{p\times p}$ and ${\mathrm{corr}}^{\mathcal{K}}\left(\mathbf{y}\right)\in {\mathbb{R}}^{p\times p}$, respectively.