Sharper rates of convergence for the tensor graphical Lasso estimator

Matrix and tensor-valued data with complex dependencies are ubiquitous in modern statistics and machine learning, flowing from sources as diverse as medical and radar imaging modalities, spatial-temporal and meteorological data collected from sensor networks and weather stations, and biological, neuroscience and spatial gene expression data aggregated over trials and time points. Learning useful structures from these large scale, complex and high-dimensional data in the low sample regime is an important task. Undirected graphs are often used to describe high dimensional distributions. Under sparsity conditions, the graph can be estimated using $\ell_{1}$-penalization methods, such as the graphical Lasso (GLasso) [12] and multiple (nodewise) regressions [27]. Under suitable conditions, including the independence among samples, such approaches yield consistent (and sparse) estimation in terms of graphical structure and fast convergence rates with respect to the operator and Frobenius norm for the covariance matrix and its inverse. The independence assumptions as mentioned above substantially simplify mathematical derivations but they tend to be very restrictive.

To remedy this, recent work has demonstrated another regime where further improvements in the sample size lower bounds are possible under additional structural assumptions, which arise naturally in the above mentioned contexts for data with complex dependencies. For example, the matrix-normal model [7] as studied in [1, 24, 35, 43] restricts the topology of the graph to tensor product graphs where the precision matrix corresponds to a Kronecker product over two component graphs. In [43], the author showed that one can estimate the covariance and inverse covariance matrices well using only one instance from the matrix variate normal distribution. However, such a normality assumption is also not needed, as elaborated in a recent paper by the same author [44]. More specifically, while the precision matrix encodes conditional independence relations for Gaussian distributions, for the more general sub-gaussian matrix variate model, this no longer holds. However, the inverse covariance matrix still encodes certain zero correlation relations between the residual errors and the covariates in a regression model, analogous to the commonly known Gaussian regression model [23]. See [44], where such regression model is introduced for sub-gaussian matrix variate data. See also [17] and references therein for more recent applications of matrix variate models in genomics.

Along similar lines, the Bigraphical Lasso framework was proposed to parsimoniously model conditional dependence relationships of matrix variate data based on the Cartesian product of graphs. As pointed out in [19], the associativity of the Kronecker sum yields an approach to the modeling of datasets organized into 3 or higher-order tensors. In [16], we explored this possibility to a great extent, by (a) introducing the tensor graphical Lasso (TeraLasso) procedure for estimating sparse $L$-way decomposable inverse covariance matrices for all $L\geq 2$; and (b) showing the rates of convergence in the Frobenius and the operator norm for estimating this class of inverse covariance matrices for sub-gaussian tensor-valued data

Consider the mean zero $L$-order random tensor $\boldsymbol{\mathscr{X}}\in\mathbb{R}^{d_{1}\times\dots\times d_{L}}$, and assume that we are given $n$ independent samples $\boldsymbol{\mathscr{X}}_{1},\ldots,\boldsymbol{\mathscr{X}}_{n}\sim\boldsymbol{\mathscr{X}}$. Here $\sim$ represents that two vectors follow the same distribution. Denote by $\mathbf{p}=[d_{1},\dots,d_{L}]$ the vector of component dimensions and $p$ the product of $d_{j}$s. Hence

	$\displaystyle\mathrm{vec}(\boldsymbol{\mathscr{X}})\in\mathbb{R}^{p},\quad\text{where }\quad p=\prod\nolimits_{k}d_{k}\quad\text{and}$
	$\displaystyle m_{k}=\prod\nolimits_{i\neq k}d_{i}={p}/{d_{k}}$		(1)

is the effective sample size we have to estimate the relations among the $d_{k}$ features along the $k^{th}$ mode in the tensor model. It was shown that due to the element replication inherent in the Cartesian product structure, the precision matrix in the TeraLasso model can be accurately estimated from limited data samples of high dimensional variables with multiway coordinates such as space, time and replicates. In particular, single sample convergence was proved for $L>2$ and empirically observed for all $L$. In contrast, direct application of the models in [12, 45] and the analysis frameworks in [30, 45, 46] require the sample size $n$ to scale proportionally to $p$, when the total number of edges is $\propto p$, which is still often too large to be practical. As a result, it is common to assume certain axes of $\boldsymbol{\mathscr{X}}$ are i.i.d., often an overly simplistic model.

For a random variable $Z$, the subgaussian (or $\psi_{2}$) norm of $Z$ denoted by $\left\lVert Z\right\rVert_{\psi_{2}}$ is defined as $\left\lVert Z\right\rVert_{\psi_{2}}=\inf\{t>0\;:\;{\mathbb{E}}\exp(Z^{2}/t^{2})\leq 2\}.$ Consider the tensor-valued data $\boldsymbol{\mathscr{X}}$ generated from a subgaussian random vector $Z=(Z_{j})\in{\mathbb{R}}^{p}$ with independent mean-zero unit variance components whose $\psi_{2}$ norms are uniformly bounded:

	$\displaystyle\rm{vec}\{\,\boldsymbol{\mathscr{X}}\,\}$	$\displaystyle=$	$\displaystyle\Sigma_{0}^{1/2}Z,\;\;\text{ where }\;$		(4)
	$\displaystyle\mathbb{E}\left(Z_{j}\right)=0,$		$\displaystyle{\mathbb{E}}Z_{j}^{2}=1,\;\text{ and}\;\left\lVert Z_{j}\right\rVert_{\psi_{2}}\leq C_{0},\forall i,j.$

We refer to $\boldsymbol{\mathscr{X}}\in\mathbb{R}^{d_{1}\times\dots\times d_{L}}$ as an order-$L$ sub-gaussian random tensor with covariance $\Sigma_{0}\in\mathbb{R}^{p\times p}$ throughout this paper. Let $\Omega_{0}=\Sigma_{0}^{-1}$. We assume that the precision matrix $\Omega_{0}=\Psi_{1}\oplus\dots\oplus\Psi_{L}$ of $\boldsymbol{\mathscr{X}}$ is the $L$-way Kronecker sum of matrix components $\{\Psi_{k}\}_{k=1}^{L}$. As such, we have

	$\displaystyle\Omega_{0}$	$\displaystyle=$	$\displaystyle\sum_{k=1}^{L}I_{[d_{1:k-1}]}\otimes\Psi_{k}\otimes I_{[d_{k+1:L}]},$		(5)
	where		$\displaystyle{I}_{[d_{k:\ell}]}:=\underbrace{{I}_{d_{k}}\otimes\dots\otimes{I}_{d_{\ell}}}_{\ell-k+1\;\mathrm{factors}},$

where $\otimes$ denotes the Kronecker (direct) product and $\ell\geq k$. The precision matrix (5) has an immediate connection to the $L$ positive-semidefinite Gram matrices $S_{n}^{k}\succeq 0\in R^{d_{k}\times d_{k}}$ associated with each mode of the tensor $\boldsymbol{\mathscr{X}}$, through tensor unfolding. Let $\boldsymbol{\mathscr{X}}_{1},\ldots,\boldsymbol{\mathscr{X}}_{n}\in\mathbb{R}^{d_{1}\times\dots\times d_{L}}\sim\boldsymbol{\mathscr{X}}$ be $n$ i.i.d. random tensors following (4). We define the mode-$k$ Gram matrix $S_{n}^{k}$ and its corresponding factor-wise marginal covariance $\Sigma_{0}^{(k)}=\mathbb{E}[S^{k}]$ as follows:

	$\displaystyle S_{n}^{k}=\frac{1}{nm_{k}}\sum_{i=1}^{n}{\bf{X}}^{(k,i)}[{\bf{X}}^{(k,i)}]^{T}\quad\mathrm{and}$		(6)
	$\displaystyle\Sigma_{0}^{(k)}=\frac{1}{m_{k}}\mathbb{E}[{\bf{X}}^{(k)}{\bf{X}}^{(k)T}],\quad k=1,\dots,L;$		(7)

See [15]. Denote by $\mathcal{K}_{\mathbf{p}}^{\sharp}$ the set of positive definite matrices that are decomposable into a Kronecker sum of fixed factor dimensions $\mathbf{p}=[d_{1},\ldots,d_{L}]$:

	$\displaystyle\quad\quad\quad\mathcal{K}_{\mathbf{p}}^{\sharp}=\{A\succ 0|A\in\mathcal{K}_{\mathbf{p}}\subset\mathbb{R}^{p\times p}\},$		(8)
	$\displaystyle\mathcal{K}_{\mathbf{p}}=\{{A}:\exists\>{B}_{k}\in\mathbb{R}^{d_{k}\times d_{k}}\>\mathrm{s.t.}\>{A}={B}_{1}\oplus\dots\oplus{B}_{L}\}.$

The TeraLasso estimator [16] minimizes the negative $\ell_{1}$-penalized Gaussian log-likelihood function $Q(\Omega)$ over the domain $\mathcal{K}_{\mathbf{p}}^{\sharp}$ of precision matrices $\Omega$, where

	$\displaystyle\quad Q(\Omega):=-\log\left\lvert\Omega\right\rvert+\;\langle{\,\widehat{S},\Omega\,}\rangle\;+\sum_{k=1}^{L}m_{k}\rho_{n,k}\left\lvert{\Psi}_{k}\right\rvert_{1,{\rm off}},$		(9)
	$\displaystyle\quad\quad\quad\text{ where }\;\;\widehat{S}=\frac{1}{n}\sum_{i=1}^{n}\rm{vec}\{\,\boldsymbol{\mathscr{X}}_{i}^{T}\,\}(\rm{vec}\{\,\boldsymbol{\mathscr{X}}_{i}^{T}\,\})^{T},$		(10)
	$\displaystyle\text{and for each $k$},\quad\left\lvert{\Psi}_{k}\right\rvert_{1,{\rm off}}=\sum_{i\neq j}{\left\lvert\Psi\right\rvert_{k,ij}},\text{ and }\quad\rho_{n,k}>0$

is a penalty parameter to be specified. Here, the weight $m_{k}=p/d_{k}$ (1) for each $k$ is determined by the number of times for which a structure $\Psi_{k}$ is replicated in $\Omega_{0}$. Then for $\mathcal{K}_{\mathbf{p}}^{\sharp}$ as in (8), we have

	$\displaystyle\quad\quad\quad\text{(TeraLasso)}\quad\widehat{\Omega}:=\mathop{\text{arg\,min}\kern 0.86108pt}_{\Omega\in\mathcal{K}_{\mathbf{p}}^{\sharp}}Q(\Omega)=$
			$\displaystyle\mathop{\text{arg\,min}\kern 0.86108pt}_{\Omega\in\mathcal{K}_{\mathbf{p}}^{\sharp}}\big{(}-\log|\Omega|+\sum_{k=1}^{L}m_{k}\big{(}\langle S_{n}^{k},\Psi_{k}\rangle+\rho_{n,k}\left\lvert{\Psi}_{k}\right\rvert_{1,{\rm off}}\big{)}\big{)}.$

The objective function (2) depends on the training data via the coordinate-wise Gram matrices $S_{n}^{k}$, where we replace the trace term $\;\langle{\,\widehat{S},\Omega\,}\rangle\;$ in (9) with the weighted sum over component-wise trace terms in (2), in view of (5) and (6), cf. Lemma 2.1. Here and in [16], the set of penalty parameters $\{\rho_{n,k},k=1,\ldots,L\}$ are chosen to dominate the maximum of entrywise errors for estimating the population $\Sigma_{0}^{(k)}$ (7) with sample ${S}^{k}_{n}$ as in (6), for each $k\leq L$ on event $\mathcal{T}$; cf. (12). This choice works equally well for the subgaussian model (4). For $L=2$ and $\Omega_{0}=\Psi_{1}\oplus\Psi_{2}=\Psi_{1}\otimes I_{d_{1}}+I_{d_{2}}\otimes\Psi_{2}$, the objective function (2) is similar in spirit to the BiGLasso objective [19], where $S^{k}_{n},k=1,2$ correspond to the Gram matrices computed from row and column vectors of matrix variate samples $X_{1},\ldots,X_{n}\in{\mathbb{R}}^{d_{1}\times d_{2}}$ respectively. When $\Omega_{0}=\Psi_{1}\otimes\Psi_{2}$ is a Kronecker product rather than a Kronecker sum over the factors, the objective function (2) is also closely related to the Gemini estimators in [43], where $\log\left\lvert\Omega_{0}\right\rvert$ is a linear combination of $\log\left\lvert\Psi_{k}\right\rvert,k=1,2$. When $\boldsymbol{\mathscr{X}}$ follows a multivariate Gaussian distribution and the precision matrix $\Omega_{0}$ has a decomposition of the form (5), the sparsity pattern of $\Psi_{k}$ for each $k$ corresponds to the conditional independence graph across the $k^{\text{th}}$ dimension of the data. Similar to the graphical Lasso, incorporating an $\ell_{1}$-penalty promotes a sparse graphical structure in the $\Psi_{k}$ and by extension $\widehat{\Omega}$. See for example [6, 3, 42, 46, 29, 19, 43, 16] and references therein.

In the present work, due to the tighter error bound on the diagonal component of the loss function as stated in Lemma 2.2, we achieve the sharper rates of convergence in Theorem 3.1, which significantly improve upon earlier results in [16] as stated in Theorem 3.2. Specifically, we replace the $p\log p$ in the earlier factor with $p$ for the relative errors in the operator and Frobenius norm in Theorem 3.1 in the present work. Under suitable assumptions on the sparsity parameter $s$ and dimensions $d_{k},\forall k\in[L]$, consistency and the rate of convergence in the operator norm can be obtained for all $n$ and $L$. First, we state the set of assumptions, (A1), (A2) and (A3).

• (A1)

  Let $\min_{k}m_{k}\geq\log p$. Denote by $\delta_{n,k}\asymp\left\lVert\Sigma_{0}\right\rVert_{2}\sqrt{\frac{\log p}{nm_{k}}}$ for $k=1,\ldots,L$. Suppose $\rho_{n,k}=\delta_{n,k}/\varepsilon_{k}$, where $0<\varepsilon_{k}<1$ for all $k$.

• (A2)

  The smallest eigenvalue $\phi_{\min}({\Omega}_{0})=\sum_{k=1}^{L}\phi_{\min}({\Psi}_{k})\geq\underline{k}_{\Omega}>0$, and the largest eigenvalue $\phi_{\max}({\Omega}_{0})=\sum_{k=1}^{L}\phi_{\max}({\Psi}_{k})\leq\overline{k}_{\Omega}<\infty$.

• (A3)

  The sample size $n$ satisfies the following:

  $\displaystyle n(m_{\min})^{2}\geq C^{2}(L+1)\kappa(\Sigma_{0})^{4}(s\log p+Lp),$

  where $m_{\min}:=\min_{k}m_{k}$, $s=\sum_{k}m_{k}s_{k}$ is as in Definition 2.3, and $C$ is an absolute constant.

The condition number for $\Sigma_{0}=\Omega_{0}^{-1}$ is defined as

$\displaystyle\kappa(\Sigma_{0})=\kappa(\Omega_{0})=\left\lVert\Omega_{0}\right\rVert_{2}\left\lVert\Omega_{0}^{-1}\right\rVert_{2}=\frac{\sum_{k=1}^{L}\phi_{\max}({\Psi}_{k})}{\sum_{k=1}^{L}\phi_{\min}({\Psi}_{k})},$

where we have used the additivity of the eigenvalues of the Kronecker sum. Following [16], we focus on error bounds on the estimate of $\Omega_{0}$ itself, rather than the individually factors in the present work. We emphasize that we retain essentially the same error bound as that in [16] for the off-diagonal component of the trace terms in (2). Denote by $\mathcal{T}=\bigcap_{i=1}^{L}\mathcal{T}_{k}$, where

	$\displaystyle\text{ event }\;\mathcal{T}_{k}=\left\{\max_{i\not=j}\left\lvert{S}^{k}_{n,ij}-\Sigma_{0,ij}^{(k)}\right\rvert\leq\delta_{n,k}\right\},$		(12)
	$\displaystyle\text{ for }\delta_{n,k}\asymp\left\lVert\Sigma_{0}\right\rVert_{2}\sqrt{{\log p}/{(nm_{k})}}>0.$		(13)

Event $\mathcal{T}$ is needed to control the off-diagonal component of the loss function as in the supplementary Lemma A.3, cf. Proof of Lemma 12 [15].

Intuitively, we use $nm_{k}$ fibers to estimate relations between and among the $d_{k}$ features along the $k^{th}$ mode as encoded in $\Psi_{k}$ and this allows optimal statistical rates of convergence to be derived, in terms of entrywise errors for estimating $\Sigma_{0}^{(k)}$ (7) with ${S}^{k}_{n}$ (6). Correspondingly, events $\{\mathcal{T}_{k},k=1,\ldots,L\}$ in (12), where were originally defined in [16], are used in the present work that reflect this sample aggregation with $nm_{k}$ being the effective size for estimating $\Psi_{k}$. Indeed, as we will show in Theorem 3.2 [16], these entrywise error bounds already enabled a significant improvement in the sample size lower bound in order to estimate parameters and the associated conditional independence graphs along different coordinates such as space, time and experimental conditions. However, these entrywise error bounds are not sufficient to achieve the operator norm-type of bounds as in Theorem 3.1 for inverse covariance estimation. In fact, using the entrywise error bounds to control the diagonal components of the trace terms will result in an extra $\log p$ factor in the sample size lower bound and correspondingly a slower rate of convergence. This extraneous $\log p$ factor is undesirable since the diagonal component of the loss function dominates the overall rate of convergence in sparse settings for inverse covariance estimation.

This is consistent with the successful finite sample experiments in both [16] and the present work, where for $L=2$, bounded errors in the operator norm are observed as $p$ increases. As a result, our new bound supports the use of the TeraLasso estimator when $L=2$, so long as a small number of replicates are available, that is, when $n=o(\log p)$, in a way that the previous Theorem 3.2 cannot. More precisely, for finite sample settings, namely, when $n=O(1)$, the relative errors will still be bounded at $O_{p}(1)$ for $L=2$, for example, when the two dimensions are at the same order: $d_{1}\asymp d_{2}$, and rapidly converge to zero for $L>2$; cf. Theorem 3.4.

Under the bounded aspect ratio regime, the relative errors in the operator and Frobenius norm for estimating the precision matrix $\Omega_{0}$ depend on the decay of the worst aspect ratio ${d_{\max}}/{m_{\min}}$ and the average of ${s_{k}\log p}/{d_{k}}$ over all modes, which represents relative sparsity levels (sparsity / dimension) in an average sense. For $L>2$, typically the aspect ratio is much less than 1 and convergence happens rapidly. If the sparse support set is small relative to nominal dimension $d_{k}$ along each mode, for example, when $\frac{s_{k}\log p}{d_{k}}=O(1)$, this convergence is at the rate of decay of the worst aspect ratio. In this case, the diagonal component dominates the rate of convergence and this is essentially optimal, since in the largest component with dimension $d_{\max}$, it has $d_{\max}$ parameters to be estimated and $m_{\min}=p/d_{\max}$ effective samples for the task. Moreover, $L$ is needed in the bound since we estimate $L$ components all together using one sample in case $n=1$. As a first example, we consider the cubic setting where $d_{1}\asymp\ldots\asymp d_{L}\asymp p^{1/L}$.

The proof of Theorems 3.1, 3.4, and Corollary 3.3 appear in the supplementary Section A.

Throughout this proof, we assume $n=1$ for simplicity. We now provide outline for proving the upper bound on the diagonal component of the main result of the paper. Recall the true parameter $\Omega_{0}=\Psi_{1}\oplus\dots\oplus\Psi_{L}$, where $\Psi_{k}\in{\mathbb{R}}^{d_{k}\times d_{k}}$ (5). Let $\Omega\in\mathcal{K}_{\mathbf{p}}$. Since $\Omega_{0}\in\mathcal{K}_{\mathbf{p}}$, we have

$\displaystyle\Delta_{\Omega}$

$\displaystyle:=$

$\displaystyle\Omega-\Omega_{0}=\Delta_{\Psi_{1}}\oplus\Delta_{\Psi_{2}}\oplus\ldots\oplus\Delta_{\Psi_{L}},$

(16)

for some $\Delta_{\Psi_{k}}\in{\mathbb{R}}^{d_{k}\times d_{k}}$ whose off-diagonal (but not diagonal) elements are uniquely determined. For self-containment, we state Lemma 4.2, where we also state the notation we use throughout this section. Here we use the trace-zero convention which guarantees the uniqueness of the $\Delta^{\prime}_{\Psi_{k}}$ in (17). We will then restate Lemma 2.2 in Lemma 4.3. The off-diagonal component has been dealt with in [16], cf. the supplementary Lemma A.3. Hence in the rest of this section, we prove Lemma 4.3. First we show the following bounds on the $\varepsilon$-net of $\mathbb{S}^{d_{k}-1},\forall k$.