Package “Tlasso”

Estimation: Tlasso Algorithm. Lyu et al. (2019)

A tensor \({\cal T} \in \mathbb R^{m_1 \times m_2 \times \cdots \times m_K}\) follows the tensor normal distribution with zero mean and covariance matrices \(\bf{\Sigma}_1, \ldots, \bf{\Sigma}_K\), denoted as \({\cal T} \sim \textrm{TN}({\bf0}; \bf{\Sigma}_1, \ldots, \bf{\Sigma}_K)\), if its probability density function is \[ p({\cal T}| \bf{\Sigma}_1,\ldots,\bf{\Sigma}_K) = (2\pi)^{-m/2} \biggl\{ \prod_{k=1}^K |\bf{\Sigma}_k|^{-m/(2m_k)} \biggr\} \exp \big(- \|{\cal T} \times \bf{\Sigma}^{-1/2}\|_F^2/2 \big), \] where \(m = \prod_{k=1}^K m_k\) and \(\bf{\Sigma}^{-1/2} := \{\bf{\Sigma}_1^{-1/2},\ldots,\bf{\Sigma}_K^{-1/2}\}\).

A standard approach to estimate precision matrix \(\bf{\Omega}_k^*\), \(k=1,\ldots,K\), is to use the maximum likelihood method. Up to a constant, the negative log-likelihood function of the tensor normal distribution is \[ \ell(\bf{\Omega}_1, \ldots, \bf{\Omega}_K) := \frac{1}{2}\textrm{tr}[\bf{S} (\bf{\Omega}_K \otimes \cdots \otimes \bf{\Omega}_1)] - \frac{1}{2}\sum_{k=1}^K \frac{m}{m_k} \log |\bf{\Omega}_k|, \] where \(\bf{S} := \frac{1}{n} \sum_{i=1}^n \textrm{vec}({\cal T}_i) \textrm{vec}({\cal T}_i)^{\top}\). To encourage the sparsity of each precision matrix in the high-dimensional scenario, a penalized log-likelihood estimator is proposed which minimizes \[ q_n(\bf{\Omega}_1, \ldots, \bf{\Omega}_K) := \frac{1}{m}\textrm{tr}[\bf{S} (\bf{\Omega}_K \otimes \cdots \otimes \bf{\Omega}_1)] - \sum_{k=1}^K \frac{1}{m_k} \log |\bf{\Omega}_k| + \sum_{k=1}^K P_{\lambda_k}(\bf{\Omega}_k), \] where \(P_{\lambda_k}(\bf{\Omega}_k) = \lambda_k \sum_{i\ne j} |[\bf{\Omega}_{k}]_{i,j}|\) is penalty function and \(\lambda_k\) is tuning parameter.

\(q_n(\bf{\Omega}_1, \ldots, \bf{\Omega}_K)\) is jointly non-convex with respect to \(\bf{\Omega}_1, \ldots, \bf{\Omega}_K\). Nevertheless, \(q_n(\bf{\Omega}_1, \ldots, \bf{\Omega}_K)\) is bi-convex problem since \(q_n(\bf{\Omega}_1, \ldots, \bf{\Omega}_K)\) is convex in \(\bf{\Omega}_k\) when the rest \(K-1\) precision matrices are fixed. According to its bi-convex property, this package solves this non-convex problem by alternatively update one precision matrix with other matrices fixed. Note that, for any \(k = 1,\ldots, K\), minimizing \(q_n(\bf{\Omega}_1, \ldots, \bf{\Omega}_K)\) with respect to \(\bf{\Omega}_k\) while fixing the rest \(K-1\) precision matrices is equivalent to minimizing \[ L(\bf{\Omega}_k) := \frac{1}{m_k}\textrm{tr}(\bf{S}_k \bf{\Omega}_k) - \frac{1}{m_k} \log |\bf{\Omega}_k| + \lambda_k \|\bf{\Omega}_k\|_{1,\textrm{off}}. \] Here \(\bf{S}_k := \frac{m_k}{n m}\sum_{i=1}^n \bf{V}_i^k \bf{V}_i^{k\top}\), where \(\bf{V}_i^k := \big[ {\cal T}_i \times \big\{\bf{\Omega}_1^{1/2},\ldots,\bf{\Omega}_{k-1}^{1/2}, 1_{m_k}, \bf{\Omega}_{k+1}^{1/2},\ldots,\bf{\Omega}_{K}^{1/2} \big\} \big]_{(k)}\) with \(\times\) the tensor product operation and \([\cdot]_{(k)}\) the mode-\(k\) matricization operation defined in . Note that minimizing \(L(\bf{\Omega}_k)\) corresponds to estimating vector-valued Gaussian graphical model and can be solved efficiently via the glasso algorithm Friedman et al. (2008). See Lyu et al. (2019) for detailed algorithm, named Tlasso.

Inference:

Multiple testing method is established to testing the entries of precision matrix for each way of the tensor, Lyu et al. (2019). We focus on \(\bf{\Omega}_1^*\), and procedures on the rest \(K-1\) precision matrices are symmetric. Hypothesis for \(\bf{\Omega}_1^*\) is, \(\forall \,1 \le i < j \le m_1\), \[ H_{0 1 , i j } : \; [{\bf{\Omega}}_1^*]_{i,j} =0 \quad\quad \textrm{vs} \quad\quad H_{11 , ij} : \; [\bf{\Omega}_1^*]_{i,j} \ne 0 \]

Let index \(-i\) in \(k\)-th mode denote all elements of mode-\(k\) except the \(i\)-th one. From the distribution of \({\cal T }_{: , i_2 , \ldots , i_K}\), we have, for every \(i_k \in \{ 1 , \ldots , m_k\}\), \(k \in \{ 1 , \ldots , K\}\), \({\cal T}_{i_1 , i_2 , \ldots , i_K} | {\cal T}_{-i_1 , i_2 , \ldots , i_K} \sim \textrm{N} ( - [\bf{\Omega}_1^*]_{i_1,i_1}^{-1} [\bf{\Omega}_1^*]_{i_1, - i_1} {\cal T}_{-i_1 , i_2 , \ldots , i_K} ; [\bf{\Omega}_1^*]_{i_1, i_1 }^{-1} \prod_{k=2}^K[\bf{\Sigma}_k^*]_{i_k, i_k} )\). An equivalent linear model is \[ {\cal T }_{l ; i_1, i_2 , \ldots , i_K} = {\cal T }_{l ; - i_1, i_2 , \ldots , i_K}^\top \bf{\theta}_{i_1} + \xi_{l ; i_1, i_2 , \ldots , i_K}, \forall \, l \in \{ 1 , \ldots , n \} \] where \(\bf{\theta}_{i_1} = - [\bf{\Omega}_1^*]_{i_1,i_1}^{-1} [\bf{\Omega}_1^*]_{i_1, - i_1} \text{ , and }\, \xi_{l ;i_1, i_2 , \ldots , i_K} \sim \textrm{N} (0 \, ; [\bf{\Omega}_1^*]_{i_1, i_1 }^{-1} \prod_{k=2}^K [\bf{\Sigma}_k^*]_{i_k, i_k}).\) Note that intercept term is eliminated since zero mean of \({\cal T }_{: , i_2 , \ldots , i_K}\).

Test statistic is constructed by both bias and variance correction of sample covariance of residuals. Let \(\hat{\bf{\theta}}_{i_1} = (\hat{\theta}_{1, i_1} , \ldots , \hat{\theta}_{m_1-1, i_1})^\top = - [\hat{\bf{\Omega}}_1]_{i_1,i_1}^{-1} [\hat{\bf{\Omega}}_1]_{i_1, - i_1}\) be the estimator of \(\bf{\theta}_{i_1}\), where \(\hat{\bf{\Omega}}_1\) is the output of Tlasso Algorithm. Given \(\{ \hat{\bf{\theta}}_{i_1} \}_{i_1 = 1 }^{m_1}\), the residual of the linear model is defined as \[ \hat{\xi}_{l ; i_1, i_2 , \ldots , i_K} = {\cal T }_{l ; i_1, i_2 , \ldots , i_K} - \bar{{\cal T }}_{ i_1, i_2 , \ldots , i_K} - ( {\cal T }_{l ; - i_1, i_2 , \ldots , i_K} - \bar{{\cal T }}_{ - i_1, i_2 , \ldots , i_K} )^\top \hat{\bf{\theta}}_{i_1}, \] where \(\bar{{\cal T }} = \frac{1}{n}\sum \limits_{l=1}^{n} {\cal T }_{l }\). The sample covariance of residuals is
\[ \hat{\varrho}_{i,j}= \frac{m_1}{(n-1) m } \sum \limits_{l=1}^{n} \sum \limits_{i_2=1}^{m_2} \cdots \sum \limits_{i_K=1}^{m_K} \hat{\xi}_{l ; i, i_2 , \ldots , i_K} \hat{\xi}_{l ; j, i_2 , \ldots , i_K} . \] Test statistic is \[\tau_{i,j} = \frac{\hat{\varrho}_{i.j} + \mu_{i,j}}{\varpi}, \forall 1 \le i < j \le m_1.\] The bias correction term \(\mu_{i,j}=\hat{\varrho}_{i,i} \hat{\theta}_{i , j} + \hat{\varrho}_{j,j} \hat{\theta}_{j-1, i}\) translates the mean of \(\tau_{i,j}\) to zero under \(H_{01, ij}\). The variance correction term \[ \varpi^2 = \frac{m \cdot \|\widehat{\bf{S}}_2\|_F^2 \cdots \|\widehat{\bf{S}}_K\|_F^2}{m_1 \cdot (\textrm{tr}(\widehat{\bf{S}}_2))^2 \cdots (\textrm{tr}(\widehat{\bf{S}}_K))^2}, \] plays a significant role in rescaling \(\tau_{i,j}\) into an asymptotic standard normal distribution, where $ k := {i=1}^n _i _i^{} $ is the estimate of \(\bf{\Sigma}_k\) with \(\widehat{\bf{V}}_i := \big[ {\cal T}_i \times \bigl\{\widehat{\bf{\Omega}}_1^{1/2},\ldots,\widehat{\bf{\Omega}}_{k-1}^{1/2}, 1_{m_k}, \widehat{\bf{\Omega}}_{k+1}^{1/2},\ldots,\widehat{\bf{\Omega}}_{K}^{1/2} \bigr\} \big]_{(k)}\) and \(\widehat{\bf{\Omega}}_{k}\) from Tlasso Algorithm, \(k \in \{2 , \ldots , K \}\). Lyu et al. (2019) proves that, under certain conditions, \[ \sqrt{\frac{(n-1) m }{ m_1 \hat{\varrho}_{i,i} \hat{\varrho}_{j,j} }} \tau_{i,j} \rightarrow \textrm{N} ( 0 ;1 ) \] in distribution, as \(nm/m_1 \rightarrow \infty\).

Let \(\tilde{\tau}_{i,j}= \sqrt{(n-1) m/ (m_1\hat{\varrho}_{i,i} \hat{\varrho}_{j,j}) } \tau_{i,j}\). Given the values of \(\tilde{\tau}_{ij}\) and thresholding level \(\varsigma\), we define a test procedure \(\varphi_{\varsigma} (\tilde{\tau}_{i,j})= 1 \{ |\tilde{\tau}_{i,j}| \ge \varsigma \}\) and reject \(H_{01, ij}\) if \(\varphi_{\varsigma} (\tilde{\tau}_{i,j})=1\). The definition of FDR/FDP in our notation is \[ \textrm{FDP} = \frac{| \{ (i,j) \in {\cal H}_0 : \varphi_{\varsigma} (\tilde{\tau}_{i,j})=1 \} | }{ | \{ (i,j) : 1 \le i < j \le m_1, \varphi_{\varsigma} (\tilde{\tau}_{i,j})=1 \} | \vee 1 } \, \text{ , and } \,\textrm{FDR} = \bf{E} (\textrm{FDP}). \] where \({\cal H}_0 = \{ (i,j) : [\bf{\Omega}_1^*]_{i,j} = 0 , 1 \le i < j \le m_1 \}\). To recover the support of \(\bf{\Omega}_1^*\), \((m_1-1)m_1/2\) tests are simultaneously conducted. Hence, a small enough \(\varsigma\) that enhances power while controls FDP under a pre-specific level \(\upsilon \in (0,1)\) is an ideal choice. In particular, \(\varsigma_{*} = \inf \{ \varsigma > 0: \text{FDP} \le \upsilon \}\). However, \(\varsigma_*\) is oracle since we have no access to the information of \({\cal H}_0\). Due to sparsity, \(w_0 =|{\cal H}_0|\) can be approximated by \(w= m_1(m_1 -1 )/2\). \(P(\varphi_{\varsigma} (\tilde{\tau}_{i,j})=1)\) is close to \(2(1 - \Phi( \varsigma))\) by test consistency. The approximation of \(\varsigma_*\) is \[ \hat{\varsigma}= \inf \bigg \{ \varsigma > 0: \frac{2(1-\Phi( \varsigma )) w}{ | \{ (i,j) : 1 \le i < j \le m_1, \varphi_{\varsigma} (\tilde{\tau}_{i,j})=1 \} | \vee 1 } \le \upsilon \bigg \}. \]

FDR/FDP control procedure is to infer the support of \(\bf{\Omega}_1^*\) at thresholding level \(\hat{\varsigma}\). In particular, FDR and FDP in the procedure are \[ \textrm{FDP}_1 = \frac{| \{ (i,j) \in {\cal H}_0 : \varphi_{\hat{\varsigma}} (\tilde{\tau}_{i,j})=1 \} | }{ | \{ (i,j) : 1 \le i < j \le m_1, \varphi_{\hat{\varsigma}} (\tilde{\tau}_{i,j})=1 \} | \vee 1 } \text{ , and } \,\textrm{FDR}_1 = \bf{E} (\textrm{FDP}_1). \] Note that we set \(\tilde{\tau}_{i,j} = \tilde{\tau}_{j,i}\) for \(1 \le i < j \le m_1\). The inference of \(\textrm{supp}(\bf{\Omega}_k^*) , \, k \in \{2 , \ldots , K \}\), is a symmetric procedure.

Lyu et al. (2019) gives the asymptotic control of \(\text{FDP}_1\) and \(\text{FDR}_1\) for the support of \(\bf{\Omega}_{1}^*\): under certain conditions, \[\textrm{FDP}_1 w / \upsilon w_0 \rightarrow 1 , \; \textrm{FDR}_1 w / \upsilon w_0 \rightarrow 1 \] in probability as \(nm/m_1 \rightarrow \infty\).

Back to Top