# Learning High-Dimensional Differential Graphs From Multi-Attribute Data

Jitendra K. Tugnait J.K. Tugnait is with the Department of Electrical &
Computer Engineering, 200 Broun Hall, Auburn University, Auburn, AL 36849,
USA. Email<EMAIL_ADDRESS>. This work was supported by the National
Science Foundation Grants ECCS-2040536 and CCF-2308473.

###### Abstract

We consider the problem of estimating differences in two Gaussian graphical
models (GGMs) which are known to have similar structure. The GGM structure is
encoded in its precision (inverse covariance) matrix. In many applications one
is interested in estimating the difference in two precision matrices to
characterize underlying changes in conditional dependencies of two sets of
data. Existing methods for differential graph estimation are based on single-
attribute (SA) models where one associates a scalar random variable with each
node. In multi-attribute (MA) graphical models, each node represents a random
vector. In this paper, we analyze a group lasso penalized D-trace loss
function approach for differential graph learning from multi-attribute data.
An alternating direction method of multipliers (ADMM) algorithm is presented
to optimize the objective function. Theoretical analysis establishing
consistency in support recovery and estimation in high-dimensional settings is
provided. Numerical results based on synthetic as well as real data are
presented.

###### Index Terms:

Sparse graph learning; differential graph estimation; undirected graph; multi-
attribute graphs.

##  I Introduction

Graphical models provide a powerful tool for analyzing multivariate data [1,
2]. In a statistical graphical model, the conditional statistical dependency
structure among $p$ random variables $x_{1},x_{1},\cdots,x_{p}$, is
represented using an undirected graph ${\cal G}=\left(V,{\cal E}\right)$,
where $V=\\{1,2,\cdots,p\\}=[p]$ is the set of $p$ nodes corresponding to the
$p$ random variables $x_{i}$s, and ${\cal E}\subseteq V\times V$ is the set of
undirected edges describing conditional dependencies among the components of
${\bm{x}}$. The graph ${\cal G}$ then is a conditional independence graph
(CIG) where there is no edge between nodes $i$ and $j$ (i.e.,
$\\{i,j\\}\not\in{\cal E}$) iff $x_{i}$ and $x_{j}$ are conditionally
independent given the remaining $p$-$2$ variables $x_{\ell}$, $\ell\in[p]$,
$\ell\neq i$, $\ell\neq j$. In particular, Gaussian graphical models (GGMs)
are CIGs where ${\bm{x}}$ is multivariate Gaussian. Suppose ${\bm{x}}$ has
positive-definite covariance matrix $\bm{\bm{\Sigma}}$ with inverse covariance
matrix $\bm{\Omega}=\bm{\bm{\Sigma}}^{-1}$. Then ${\Omega}_{ij}$, the
$(i,j)$-th element of $\bm{\Omega}$, is zero iff $x_{i}$ and $x_{j}$ are
conditionally independent. Such models for ${\bm{x}}$ have been extensively
studied. Given $n$ samples of ${\bm{x}}$, in high-dimensional settings where
$p\gg 1$ and/or $n$ is of the order of $p$, one estimates $\bm{\Omega}$ under
some sparsity constraints; see [3, 5, 6, 4].

More recently there has been increasing interest in differential network
analysis where one is interested in estimating the difference in two inverse
covariance matrices [8, 9, 7]. Given observations ${\bm{x}}$ and ${\bm{y}}$
from two groups of subjects, one is interested in the difference
${\bm{\Delta}}={\bm{\Omega}}_{y}-{\bm{\Omega}}_{x}$, where
${\bm{\Omega}}_{x}=(E\\{{\bm{x}}{\bm{x}}^{\top}\\})^{-1}$ and
${\bm{\Omega}}_{y}=(E\\{{\bm{y}}{\bm{y}}^{\top}\\})^{-1}$. The associated
differential graph is ${\cal G}_{\Delta}=\left(V,{\cal E}_{\Delta}\right)$
where $\\{i,j\\}\in{\cal E}_{\Delta}$ iff ${\bm{\Delta}}_{ij}\neq 0$. It
characterizes differences between the GGMs of the two sets of data. We use the
term differential graph as in [10, 11] ([8, 9, 7] use the term differential
network). As noted in [9], in biostatistics, the differential network/graph
describes the changes in conditional dependencies between components under
different environmental or genetic conditions. For instance, one may be
interested in the differences in the graphical models of healthy and impaired
subjects, or models under different disease states, given gene expression data
or functional MRI signals [3, 12, 13].

In the preceding graphs, each node represents a scalar random variable. In
many applications, there may be more than one random variable associated with
a node. This class of graphical models has been called multi-attribute (MA)
graphical models in [15, 16, 14, 17] and vector graphs or networks in [18, 19,
20, 21]. In a gene regulatory network, one may have different molecular
profiles available for a single gene, such as protein, DNA and RNA. Since
these molecular profiles are on the same set of biological samples, they
constitute multi-attribute data for gene regulatory graphical models in [16,
14]. Consider $p$ jointly Gaussian vectors ${\bm{z}}_{i}\in\mathbb{R}^{m}$,
$i\in[p]$. We associate ${\bm{z}}_{i}$ with the $i$th node of graph ${\cal
G}=\left(V,{\cal E}\right)$, $V=[p]$, ${\cal E}\subseteq V\times V$. We now
have $m$ attributes per node. Now $\\{i,j\\}\in{\cal E}$ iff vectors
${\bm{z}}_{i}$ and ${\bm{z}}_{j}$ are conditionally independent given the
remaining $p$-$2$ vectors $\\{{\bm{z}}_{\ell}\,,\ell\in
V\textbackslash\\{i_{,}j\\}\\}$. Let
${\bm{x}}=[{\bm{z}}_{1}^{\top}\;{\bm{z}}_{2}^{\top}\;\cdots\;{\bm{z}}_{p}^{\top}]^{\top}\in\mathbb{R}^{mp}$.
Let ${\bm{\Omega}}=(E\\{{\bm{x}}{\bm{x}}^{\top}\\})^{-1}$ assuming
$E\\{{\bm{x}}{\bm{x}}^{\top}\\}\succ{\bm{0}}$. Define the $m\times m$ subblock
${\bm{\Omega}}^{(ij)}$ of ${\bm{\Omega}}$ as
$[{\bm{\Omega}}^{(ij)}]_{rs}=[{\bm{\Omega}}]_{(i-1)m+r,(j-1)m+s}\,,\;r,s=1,2,\cdots,m$.
Then we have the following equivalence [16, Sec. 2.1]

$\\{i,j\\}\not\in{\cal E}\;\Leftrightarrow\;{\bm{\Omega}}^{(ij)}={\bm{0}}\,.$
(1)

This paper is concerned with estimation of differential graphs from multi-
attribute data. Given independent and identically distributed (i.i.d.) samples
${\bm{x}}(t)$, $t=1,2,\cdots,n_{x}$, of
${\bm{x}}=[{\bm{z}}_{1}^{\top}\;{\bm{z}}_{2}^{\top}\;\cdots\;{\bm{z}}_{p}^{\top}]^{\top}\in\mathbb{R}^{mp}$
where ${\bm{z}}_{i}\in\mathbb{R}^{m}$, $i\in[p]$, are jointly Gaussian, and
similarly given samples ${\bm{y}}(t)$, $t=1,2,\cdots,n_{y}$, of
${\bm{y}}\in\mathbb{R}^{mp}$, our objective is to estimate the difference
${\bm{\Delta}}={\bm{\Omega}}_{y}-{\bm{\Omega}}_{x}$, and determine the
differential graph ${\cal G}_{\Delta}=\left(V,{\cal E}_{\Delta}\right)$ with
edgeset ${\cal
E}_{\Delta}=\\{\\{k,\ell\\}\,:\,\|{\bm{\Delta}}^{(k\ell)}\|_{F}\neq 0\\}$.

###  I-A Related Work

All prior work on high-dimensional differential graph estimation from i.i.d.
samples addresses single-attribute (SA) models where each node represents a
scalar random variable. One naive approach would be to estimate the two
precision matrices separately by any existing estimator (see [5, 4] and
references therein) and then calculate their difference to estimate the
differential graph. (This approach is also applicable to MA graphs.) This
approach estimates twice the number of parameters, hence needs larger sample
sizes for same accuracy, and also imposes sparsity constraints on each
precision matrix for the methods to work. The same comment applies to methods
such as [3, 6], where the two precision matrices and their differences are
jointly estimated. A recent survey is in [22]. In these approaches, given
$K\geq 2$ related groups of data, each $p$-variate and sharing the same set of
nodes $V$, but possibly differing in connected edgesets, the objective is to
jointly estimate the $K$ precision matrices and their pairwise differences,
with sparsity constraints on each of the $K$ precision matrices and their
pairwise differences. These approaches require each of the $K$ precision
matrices to be sparse. If only the differences in the precision matrices is of
interest, alternative approaches exist where no sparsity constraints are
imposed on individual precision matrices. For instance, direct estimation of
the difference in the two precision matrices has been considered for SA graphs
in [23, 12, 25, 24, 8, 26, 27, 28, 9, 7], where only the difference is
required to be sparse, not the two individual precision matrices. In [23, 24,
8, 28, 9, 7] precision difference matrix estimators are based on a D-trace
loss [29], while [12] discusses a Dantzig selector type estimator. In [25, 26,
27] differential graph is estimated by directly modeling the ratio of the
probability densities of the random vectors under the two graphs.

Estimation of MA differential graphs has not been investigated before. The
work of [10, 11] is similar to an MA formulation except that in [10, 11],
${\bm{x}}(t)$ and ${\bm{y}}(t)$ are non-stationary (“functional” modeling),
and instead of a single record (sample) of ${\bm{x}}(t)$, $t=1,2,\cdots,n_{x}$
and ${\bm{y}}(t)$, $t=1,2,\cdots,n_{y}$, as in this paper, they assume
multiple independent observations of ${\bm{x}}(t)$, $t\in{\cal T}$ (a closed
subset of real line), and ${\bm{y}}(t)$, $t\in{\cal T}$. The objective
function in [11, Eqns. (10)-(11)] is the same as our objective function
(3)-(4), but consequent estimation of edges and theoretical analysis are
vastly different. We estimate edges as in (6), i.e., our threshold is set at
zero and this is the method analyzed in our Theorem 1(iv) for graph recovery
with high probability. In [10, 11], this threshold is set at a parameter
$\epsilon_{n}>0$ (see [11, Eqn. (13)]) which is a function of sample size $n$,
number of nonzero entries in true ${\bm{\Delta}}$, smallest eigenvalues of
true covariances ${\bm{\Omega}}_{y}^{-1}$ and ${\bm{\Omega}}_{x}^{-1}$ (in our
notation), and several other factors. That is, $\epsilon_{n}$ is unknowable
for practical implementation and it is used as a theoretical construct to
establish graph support recovery in [11, Theorem 10]. In simulations, [10, 11]
set $\epsilon_{n}=0$. That is, [10, 11] do not analyze what they implement
(the proof does not hold for $\epsilon_{n}=0$), and they do not implement what
they analyze ($\epsilon_{n}$ is unknowable). There is no counterpart to our
Theorem 1 in [10, 11], and the methodology of our Theorem 1 allows us to set
the edge detection threshold to zero. Our Theorem 2 follows the general
framework of [30] to bound the Frobenius norm of the error in estimating
${\bm{\Delta}}$, and [10, 11] also follow the general framework of [30] for
the same purpose. But their extension of this result to graph recovery does
not permit zero threshold for edge detection. We attempt no such extension.

###  I-B Our Contributions

In this paper, we analyze a group lasso penalized D-trace loss function
approach for differential graph learning from MA data, extending the SA
approach of [8, 28]. A two-block ADMM algorithm is presented to optimize the
objective function. The two-block ADMM is guaranteed to be convergent unlike
the three-block ADMM method used in [8]. Two different approaches to
theoretical analysis of the proposed approach in high-dimensional settings are
presented. Theorem 1 follows the approach(es) of [31, 16, 29, 8, 28] while
Theorem 2 follows the general framework of [30], not used in [8, 28]. The
general method of [31] requires an irrepresentability condition (see (20))
which is also required in [8, 28] for SA graphs, but is not needed by the
method of [30], hence in our Theorem 2. Numerical results based on synthetic
as well as real data are presented.

Preliminary version of parts of this paper appear in a conference paper [32].
Theorem 2, proof of Theorem 1 and real data example do not appear in [32].

###  I-C Outline and Notation

The rest of the paper is organized as follows. A group lasso penalized D-trace
loss function is presented in Sec. II for estimation of multi-attribute
differential graph. An ADMM algorithm is presented in Sec. III to optimize the
convex objective function. In Sec. IV we analyze the properties of the
estimator of the difference
${\bm{\Delta}}={\bm{\Omega}}_{y}-{\bm{\Omega}}_{x}$. Theorem 1 follows the
approach(es) of [31, 16, 29, 8, 28] while Theorem 2 follows the general
framework of [30]. The general method of [31] requires an irrepresentability
condition (see (20)) which is not needed by the method of [30]. On the other
hand, our Theorem 2 does not have a result like Theorem 1(ii), the oracle
property, nor does it have a result as in Theorem 1(iv), support recovery.
Numerical results based on synthetic as well as real data are presented in
Sec. V to illustrate the proposed approach. Proofs of Theorems 1 and 2 are
given in Appendices A and B, respectively.

For a set $V$, $|V|$ or $\mbox{card}(V)$ denotes its cardinality. Given
${\bm{A}}\in\mathbb{R}^{p\times p}$, we use $\phi_{\min}({\bm{A}})$,
$\phi_{\max}({\bm{A}})$, $|{\bm{A}}|$ and $\mbox{tr}({\bm{A}})$ to denote the
minimum eigenvalue, maximum eigenvalue, determinant and trace of ${\bm{A}}$,
respectively. For ${\bm{B}}\in\mathbb{R}^{p\times q}$, we define
$\|{\bm{B}}\|=\sqrt{\phi_{\max}({\bm{B}}^{\top}{\bm{B}})}$,
$\|{\bm{B}}\|_{F}=\sqrt{\mbox{tr}({\bm{B}}^{\top}{\bm{B}})}$,
$\|{\bm{B}}\|_{1}=\sum_{i,j}|B_{ij}|$, where $B_{ij}$ is the $(i,j)$-th
element of ${\bm{B}}$ (also denoted by $[{\bm{B}}]_{ij}$),
$\|{\bm{B}}\|_{\infty}=\max_{i,j}|B_{ij}|$ and
$\|{\bm{B}}\|_{1,\infty}=\max_{i}\sum_{j}|B_{ij}|$. The symbols $\otimes$ and
$\boxtimes$ denote Kronecker product and Tracy-Singh product [33],
respectively. In particular, given block partitioned matrices
${\bm{A}}=[{\bm{A}}_{ij}]$ and ${\bm{B}}=[{\bm{B}}_{k\ell}]$ with submatrices
${\bm{A}}_{ij}$ and ${\bm{B}}_{k\ell}$, Tracy-Singh product yields another
block partitioned matrix
${\bm{A}}\boxtimes{\bm{B}}=[{\bm{A}}_{ij}\boxtimes{\bm{B}}]_{ij}=[[{\bm{A}}_{ij}\otimes{\bm{B}}_{k\ell}]_{k\ell}]_{ij}$
[34]. Given ${\bm{A}}=[{\bm{A}}_{ij}]\in\mathbb{R}^{mp\times mp}$ with
${\bm{A}}_{ij}\in\mathbb{R}^{m\times m}$, ${\rm
vec}({\bm{A}})\in\mathbb{R}^{m^{2}p^{2}}$ denotes the vectorization of
${\bm{A}}$ which stacks the columns of the matrix ${\bm{A}}$, and

$\displaystyle{\rm bvec}({\bm{A}})=[({\rm vec}({\bm{A}}_{11}))^{\top}\;({\rm
vec}({\bm{A}}_{21}))^{\top}\;\cdots\;({\rm vec}({\bm{A}}_{p1}))^{\top}$
$\displaystyle\quad({\rm vec}({\bm{A}}_{12}))^{\top}\;\cdots\;({\rm
vec}({\bm{A}}_{p2}))^{\top}\;\cdots\;({\rm
vec}({\bm{A}}_{pp}))^{\top}]^{\top}.$

Let $S={\cal
E}_{\Delta}=\\{\\{k,\ell\\}\,:\,\|{\bm{\Delta}}^{(k\ell)}\|_{F}\neq 0\\}$
where ${\bm{\Delta}}=[{\bm{\Delta}}^{(k\ell)}]\in\mathbb{R}^{mp\times mp}$
with ${\bm{\Delta}}^{(k\ell)}\in\mathbb{R}^{m\times m}$ denoting the $(k,l)$th
$m\times m$ submatrix of ${\bm{\Delta}}$. Then ${\bm{\Delta}}_{S}$ denotes the
submatrix of ${\bm{\Delta}}$ with block rows and columns indexed by $S$, i.e.,
${\bm{\Delta}}_{S}=[{\bm{\Delta}}^{(k\ell)}]_{(k,\ell)\in S}$. Suppose
${\bm{\Gamma}}={\bm{A}}\boxtimes{\bm{B}}$ given block partitioned matrices
${\bm{A}}=[{\bm{A}}_{ij}]$ and ${\bm{B}}=[{\bm{B}}_{k\ell}]$. For any two
subsets $T_{1}$ and $T_{2}$ of $V\times V$, ${\bm{\Gamma}}_{T_{1},T_{2}}$
denotes the submatrix of ${\bm{\Gamma}}$ with block rows and columns indexed
by $T_{1}$ and $T_{2}$, i.e.,
${\bm{\Gamma}}_{T_{1},T_{2}}=[{\bm{A}}_{j\ell}\otimes{\bm{B}}_{kq}]_{(j,k)\in
T_{1},(\ell,q)\in T_{2}}$. Following [16], an operator
$\bm{\mathcal{C}}(\cdot)$ is used in Sec. IV. Consider
${\bm{A}}\in\mathbb{R}^{mp\times mp}$ with $(k,l)$th $m\times m$ submatrix
${\bm{A}}^{(k\ell)}$. Then $\bm{\mathcal{C}}(\cdot)$ operates on ${\bm{A}}$ as

$\displaystyle\begin{bmatrix}{\bm{A}}^{(11)}&\cdots&{\bm{A}}^{(1p)}\\\
\vdots&\ddots&\vdots\\\
{\bm{A}}^{(p1)}&\cdots&{\bm{A}}^{(pp)}\end{bmatrix}\overset{\bm{\mathcal{C}}(\cdot)}{\xrightarrow{\hskip
14.22636pt}}\begin{bmatrix}\|{\bm{A}}^{(11)}\|_{F}&\cdots&\|{\bm{A}}^{(1p)}\|_{F}\\\
\vdots&\ddots&\vdots\\\
\|{\bm{A}}^{(p1)}\|_{F}&\cdots&\|{\bm{A}}^{(pp)}\|_{F}\end{bmatrix}$

with $\bm{\mathcal{C}}({\bm{A}}^{(k\ell)})=\|{\bm{A}}^{(k\ell)}\|_{F}$ and
$\bm{\mathcal{C}}({\bm{A}})\in\mathbb{R}^{p\times p}$. Now consider
${\bm{A}}_{1},{\bm{A}}_{2}\in\mathbb{R}^{mp\times mp}$ with $(k,l)$th $m\times
m$ submatrices ${\bm{A}}_{1}^{(k\ell)}$ and ${\bm{A}}_{2}^{(k\ell)}$,
respectively, and Tracy-Singh product
${\bm{A}}_{1}\boxtimes{\bm{A}}_{2}\in\mathbb{R}^{(mp)^{2}\times(mp)^{2}}$.
Then $\bm{\mathcal{C}}(\cdot)$ operates on ${\bm{A}}_{1}\boxtimes{\bm{A}}_{2}$
as
$\bm{\mathcal{C}}({\bm{A}}_{1}\boxtimes{\bm{A}}_{2})\in\mathbb{R}^{p^{2}\times
p^{2}}$ with
$\bm{\mathcal{C}}({\bm{A}}_{1}^{(k_{1}\ell_{1})}\otimes{\bm{A}}_{2}^{(k_{2}\ell_{2})})=\|{\bm{A}}_{1}^{(k_{1}\ell_{1})}\otimes{\bm{A}}_{2}^{(k_{2}\ell_{2})}\|_{F}$
(=$\|{\bm{A}}_{1}^{(k_{1}\ell_{1})}\|_{F}\,\|{\bm{A}}_{2}^{(k_{2}\ell_{2})}\|_{F}$).
That is, each $m^{2}\times m^{2}$ submatrix
${\bm{A}}_{1}^{(k_{1}\ell_{1})}\otimes{\bm{A}}_{2}^{(k_{2}\ell_{2})}$ of
${\bm{A}}_{1}\boxtimes{\bm{A}}_{2}$ is mapped into its Frobenius norm.

##  II Group Lasso Penalized D-Trace Loss

Let
${\bm{x}}=[{\bm{z}}_{1x}^{\top}\;{\bm{z}}_{2x}^{\top}\;\cdots\;{\bm{z}}_{px}^{\top}]^{\top}\in\mathbb{R}^{mp}$
where ${\bm{z}}_{ix}\in\mathbb{R}^{m}$, $i\in[p]$, are zero-mean, jointly
Gaussian. Similarly, let
${\bm{y}}=[{\bm{z}}_{1y}^{\top}\;{\bm{z}}_{2y}^{\top}\;\cdots\;{\bm{z}}_{py}^{\top}]^{\top}\in\mathbb{R}^{mp}$
where ${\bm{z}}_{iy}\in\mathbb{R}^{m}$, $i\in[p]$, are zero-mean, jointly
Gaussian. Given i.i.d. samples ${\bm{x}}(t)$, $t=1,2,\cdots,n_{x}$, of
${\bm{x}}$, and similarly given i.i.d. samples ${\bm{y}}(t)$,
$t=1,2,\cdots,n_{y}$, of ${\bm{y}}\in\mathbb{R}^{mp}$, form the sample
covariance estimates

$\hat{\bm{\Sigma}}_{x}=\frac{1}{n_{x}}\sum_{t=1}^{n_{x}}{\bm{x}}(t){\bm{x}}^{\top}(t)\,,\;\;\hat{\bm{\Sigma}}_{y}=\frac{1}{n_{y}}\sum_{t=1}^{n_{y}}{\bm{y}}(t){\bm{y}}^{\top}(t)\,.$
(2)

and denote their true values as
${\bm{\Sigma}}_{x}^{\ast}={\bm{\Omega}}_{x}^{-\ast}(=({\bm{\Omega}}_{x}^{\ast})^{-1})$
and ${\bm{\Sigma}}_{y}^{\ast}={\bm{\Omega}}_{y}^{-\ast}$. Assume that
$\\{{\bm{x}}(t)\\}$ and $\\{{\bm{y}}(t)\\}$ are mutually independent
sequences. Assume ${\bm{\Sigma}}_{x}^{\ast}$ and ${\bm{\Sigma}}_{y}^{\ast}$
are positive definite. We wish to estimate
${\bm{\Delta}}={\bm{\Omega}}_{y}^{\ast}-{\bm{\Omega}}_{x}^{\ast}$ and graph
${\cal G}_{\Delta}=\left(V,{\cal E}_{\Delta}\right)$, based on
$\hat{\bm{\Sigma}}_{x}$ and $\hat{\bm{\Sigma}}_{y}$. Following the SA
formulation of [8] (see also [28, Sec. 2.1]), we will use a convex D-trace
loss function given by

$L({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})=\frac{1}{2}\mbox{tr}(\hat{\bm{\Sigma}}_{x}{\bm{\Delta}}\hat{\bm{\Sigma}}_{y}{\bm{\Delta}}^{\top})-\mbox{tr}({\bm{\Delta}}(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y}))$
(3)

where D-trace refers to difference-in-trace loss function, a term coined in
[29] in the context of graphical model estimation. The function
$L({\bm{\Delta}},{\bm{\Sigma}}_{x}^{\ast},{\bm{\Sigma}}_{y}^{\ast})$ is
strictly convex in ${\bm{\Delta}}$ (its Hessian w.r.t.
$\mbox{vec}({\bm{\Delta}})$ is
${\bm{\Sigma}}_{y}^{\ast}\otimes{\bm{\Sigma}}_{x}^{\ast}$), and has a unique
minimum at
${\bm{\Delta}}^{\ast}={\bm{\Omega}}_{y}^{\ast}-{\bm{\Omega}}_{x}^{\ast}$ [8,
28]. When we use sample covariances, we propose to estimate ${\bm{\Delta}}$ by
minimizing the group-lasso penalized loss function

$L_{\lambda}({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})=L({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})+\lambda\sum_{k,\ell=1}^{p}\|{\bm{\Delta}}^{(k\ell)}\|_{F}$
(4)

where $\lambda>0$ is a tuning parameter and $\|{\bm{\Delta}}^{(k\ell)}\|_{F}$
promotes blockwise sparsity in ${\bm{\Delta}}$ [35, 36, 37] where, if we
partition ${\bm{\Delta}}$ into $m\times m$ submatrices,
${\bm{\Delta}}^{(k\ell)}$ denotes its $(k,\ell)$th submatrix, associated with
edge $\\{k,\ell\\}$ of the differential graph ${\cal G}_{\Delta}=\left(V,{\cal
E}_{\Delta}\right)$.

For SA models ($m=1$), [28] has used the lasso-penalized loss function
$L_{J}({\bm{\Delta}})=L({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})+\lambda\sum_{k,\ell=1}^{p}\big{|}{\bm{\Delta}}_{k\ell}\big{|}$.
The cost $L_{J}({\bm{\Delta}})$ is optimized in [28] using a two-block ADMM
approach which is known to be convergent. The resulting estimator
$\hat{\bm{\Delta}}$ that minimizes the above cost is not necessarily
symmetric. To obtain a symmetric estimator for SA models, [8] proposes the
lasso-penalized loss function
$L_{Y}({\bm{\Delta}})=\frac{1}{4}\mbox{tr}(\hat{\bm{\Sigma}}_{x}{\bm{\Delta}}\hat{\bm{\Sigma}}_{y}{\bm{\Delta}}^{\top}+\hat{\bm{\Sigma}}_{y}{\bm{\Delta}}\hat{\bm{\Sigma}}_{x}{\bm{\Delta}}^{\top})-\mbox{tr}({\bm{\Delta}}(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y}))+\lambda\sum_{k,\ell=1}^{p}\big{|}{\bm{\Delta}}_{k\ell}\big{|}$.
In [8], cost $L_{Y}({\bm{\Delta}})$ is optimized using a three-block ADMM
method which is not necessarily convergent.

Suppose

$\hat{\bm{\Delta}}=\arg\min_{\bm{\Delta}}L_{\lambda}({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})\,.$
(5)

Even though ${\bm{\Delta}}$ is symmetric, $\hat{\bm{\Delta}}$ is not. We can
symmetrize it by setting
$\hat{\bm{\Delta}}_{sym}=\frac{1}{2}(\hat{\bm{\Delta}}+\hat{\bm{\Delta}}^{\top})$,
after obtaining $\hat{\bm{\Delta}}$. Then the differential graph edges are
estimated as

$\hat{\cal
E}_{\Delta}=\Big{\\{}\\{k,\ell\\}\,:\,\|\hat{\bm{\Delta}}_{sym}^{(k\ell)}\|_{F}>0\Big{\\}}\,.$
(6)

##  III Optimization

The objective function
$L_{\lambda}({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})$,
given by (4), is strictly convex. Several existing approaches such as an
alternating direction method of multipliers (ADMM) [38] or proximal gradient
descent (PGD) methods [39], can be followed to minimize (4). Note that [8, 28]
use ADMM while [9] uses a proximal gradient method, all for SA graphs. It is
stated in [8, Sec. 2.2] that in their simulation example, ADMM approach
yielded a slightly smaller value of the objective function compared to the PGD
approach. In [10, 11], similar to [9], a proximal gradient method is used for
an objective function similar to our (4). In this paper, motivated by [8], we
will develop an ADMM method. In a simulation example (Sec. V-A) we compare our
ADMM approach with ADMM and PGD approaches of [28] and [9, 11], respectively.

###  III-A ADMM Approach

Similar to [28] (also [8]), we use an ADMM approach [38] with variable
splitting. Using variable splitting, consider

$\displaystyle\min_{\bm{\Delta},{\bm{W}}}$
$\displaystyle\Big{\\{}L({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})+\lambda\sum_{k,\ell=1}^{p}\|{\bm{W}}^{(k\ell)}\|_{F}\Big{\\}}$
(7) $\displaystyle\quad\mbox{ subject to }\bm{\Delta}={\bm{W}}\,.$

The scaled augmented Lagrangian for this problem is [38]

$\displaystyle L_{\rho}=$ $\displaystyle
L({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})+\lambda\sum_{k,\ell=1}^{p}\|{\bm{W}}^{(k\ell)}\|_{F}$
$\displaystyle\quad+\frac{\rho}{2}\|{\bm{\Delta}}-{\bm{W}}+{\bm{U}}\|^{2}_{F}$
(8)

where ${\bm{U}}$ is the dual variable, and $\rho>0$ is the penalty parameter.
Given the results $\bm{\Delta}^{(i)},{\bm{W}}^{(i)},{\bm{U}}^{(i)}$ of the
$i$th iteration, in the $(i+1)$st iteration, an ADMM algorithm executes the
following three updates:

  * (a)

$\bm{\Delta}^{(i+1)}\leftarrow\arg\min_{\bm{\Delta}}\,L_{a}(\bm{\Delta}),\;\;L_{a}(\bm{\Delta}):=L({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})+\frac{\rho}{2}\|\bm{\Delta}-{\bm{W}}^{(i)}+{\bm{U}}^{(i)}\|^{2}_{F}$.

  * (b)

${\bm{W}}^{(i+1)}\leftarrow\arg\min_{{\bm{W}}}L_{b}({\bm{W}}),\;\;L_{b}({\bm{W}}):=\newline
\lambda\sum_{k,\ell=1}^{p}\|{\bm{W}}^{(k\ell)}\|_{F}+\frac{\rho}{2}\|\bm{\Delta}^{(i+1)}-{\bm{W}}+{\bm{U}}^{(i)}\|^{2}_{F}$.

  * (c)

${\bm{U}}^{(i+1)}\leftarrow{\bm{U}}^{(i)}+\left(\bm{\Delta}^{(i+1)}-{\bm{W}}^{(i+1)}\right)$.

Update (a): Differentiate $L_{a}(\bm{\Delta})$ w.r.t. ${\bm{\Delta}}$ to
obtain

$\displaystyle{\bm{0}}=$ $\displaystyle\frac{\partial
L_{a}(\bm{\Delta})}{\partial{\bm{\Delta}}}=\hat{\bm{\Sigma}}_{x}{\bm{\Delta}}\hat{\bm{\Sigma}}_{y}-(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})+\rho({\bm{\Delta}}-{\bm{W}}+{\bm{U}})$
(9) $\displaystyle\Rightarrow$
$\displaystyle(\hat{\bm{\Sigma}}_{y}\otimes\hat{\bm{\Sigma}}_{x}+\rho{\bm{I}}){\rm
vec}({\bm{\Delta}})={\rm
vec}(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y}+\rho({\bm{W}}-{\bm{U}}))$
(10)

Direct matrix inversion solution of (10) requires inversion of a
$(mp)^{2}\times(mp)^{2}$ matrix. A computationally cheaper solution is given
in [8, 28], as follows. Carry out eigendecomposition of
$\hat{\bm{\Sigma}}_{x}$ and $\hat{\bm{\Sigma}}_{y}$ as
$\hat{\bm{\Sigma}}_{x}={\bm{Q}}_{x}{\bm{D}}_{x}{\bm{Q}}_{x}^{\top}$,
${\bm{Q}}_{x}{\bm{Q}}_{x}^{\top}={\bm{I}}$ and
$\hat{\bm{\Sigma}}_{y}={\bm{Q}}_{y}{\bm{D}}_{y}{\bm{Q}}_{y}^{\top}$,
${\bm{Q}}_{y}{\bm{Q}}_{y}^{\top}={\bm{I}}$, where ${\bm{D}}_{x}$ and
${\bm{D}}_{y}$ are diagonal matrices of the respective eigenvalues. Then
$\hat{\bm{\Delta}}$ that minimizes $L_{a}(\bm{\Delta})$ is given by

$\displaystyle\hat{\bm{\Delta}}=$
$\displaystyle{\bm{Q}}_{x}\Big{[}{\bm{B}}\circ[{\bm{Q}}_{x}^{\top}\big{(}\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y}+\rho({\bm{W}}-{\bm{U}})\big{)}{\bm{Q}}_{y}]\Big{]}{\bm{Q}}_{y}^{\top}$
(11)

where the symbol $\circ$ denotes the Hadamard product and
${\bm{B}}\in\mathbb{R}^{mp\times mp}$ organizes the diagonal of
$({\bm{D}}_{y}\otimes{\bm{D}}_{x}+\rho{\bm{I}})^{-1}$ in a matrix with
${\bm{B}}_{jk}=1/([{\bm{D}}_{x}]_{jj}[{\bm{D}}_{y}]_{kk}+\rho)$. Note that the
eigendecomposition of $\hat{\bm{\Sigma}}_{x}$ and $\hat{\bm{\Sigma}}_{y}$ has
to be done only once. Thus

$\displaystyle{\bm{\Delta}}^{(i+1)}=$
$\displaystyle{\bm{Q}}_{x}\Big{[}{\bm{B}}\circ[{\bm{Q}}_{x}^{\top}\big{(}\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y}+\rho({\bm{W}}^{(i)}-{\bm{U}}^{(i)})\big{)}{\bm{Q}}_{y}]\big{]}{\bm{Q}}_{y}^{\top}$
(12)

Update (b): Here we have the group lasso solution [35, 36, 37]

$\displaystyle\quad\quad({\bm{W}}^{(k\ell)})^{(i+1)}$ $\displaystyle=$
$\displaystyle\Big{(}1-\frac{(\lambda/\rho)}{\|({\bm{\Delta}}^{(i+1)}+{\bm{U}}^{(i)})^{(k\ell)}\|_{F}}\Big{)}_{+}({\bm{\Delta}}^{(i+1)}+{\bm{U}}^{(i)})^{(k\ell)}$
(13)

where $(a)_{+}=\max(0,a)$.

Algorithm 1 ADMM Algorithm MA-ADMM

Input: Data $\\{{\bm{x}}(t)\\}_{t=1}^{n_{x}}$, ${\bm{x}}\in\mathbb{R}^{mp}$,
and $\\{{\bm{y}}(t)\\}_{t=1}^{n_{y}}$, ${\bm{y}}\in\mathbb{R}^{mp}$,
regularization and penalty parameters $\lambda$ and $\rho_{0}$, tolerances
$\tau_{abs}$ and $\tau_{rel}$, variable penalty factor $\mu$, maximum number
of iterations $i_{max}$.
Output: estimated $\hat{\bm{\Delta}}_{sym}$ and $\hat{\cal E}_{\Delta}$.

1:Calculate sample covariances
$\hat{\bm{\Sigma}}_{x}=\frac{1}{n_{x}}\sum_{t=1}^{n_{x}}{\bm{x}}(t){\bm{x}}^{\top}(t)$
and
$\hat{\bm{\Sigma}}_{y}=\frac{1}{n_{y}}\sum_{t=1}^{n_{y}}{\bm{y}}(t){\bm{y}}^{\top}(t)$.

2:Initialize: ${\bm{\Delta}}^{(0)}={\bm{U}}^{(0)}={\bm{W}}^{(0)}={\bm{0}}$,
where ${\bm{\Delta}},{\bm{U}},{\bm{W}}\in\mathbb{R}^{(mp)\times(mp)}$,
$\rho^{(0)}=\rho_{0}$.

3:Eigendecompose $\hat{\bm{\Sigma}}_{x}$ and $\hat{\bm{\Sigma}}_{y}$ as
$\hat{\bm{\Sigma}}_{x}={\bm{Q}}_{x}{\bm{D}}_{x}{\bm{Q}}_{x}^{\top}$ and
$\hat{\bm{\Sigma}}_{y}={\bm{Q}}_{y}{\bm{D}}_{y}{\bm{Q}}_{y}^{\top}$.

4:converged = FALSE, $i=0$

5:while converged = FALSE  AND  $i\leq i_{max}$, do

6: Construct ${\bm{B}}\in\mathbb{R}^{mp\times mp}$ with
${\bm{B}}_{jk}=1/([{\bm{D}}_{x}]_{jj}[{\bm{D}}_{y}]_{kk}+\rho^{(i)}$.

7: Set
${\bm{\Delta}}^{(i+1)}={\bm{Q}}_{x}\Big{[}{\bm{B}}\circ[{\bm{Q}}_{x}^{\top}\big{(}\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y}+\rho({\bm{W}}^{(i)}-{\bm{U}}^{(i)})\big{)}{\bm{Q}}_{y}]\big{]}{\bm{Q}}_{y}^{\top}\,$.

8: With $(a)_{+}:=\max(0,a)$,
${\bm{A}}=({\bm{\Delta}}^{(i+1)}+{\bm{U}}^{(i)})^{(k\ell)}$ and
$k,\ell\in[p]$, update $m\times m$ subblocks of ${\bm{W}}$ as
$\displaystyle({\bm{W}}^{(i+1)})^{(k\ell)}=$
$\displaystyle\Big{(}1-\frac{(\lambda/\rho)}{\|{\bm{A}}\|_{F}}\Big{)}_{+}{\bm{A}}^{(k\ell)}\,.$

9: Dual update
${\bm{U}}^{(i+1)}={\bm{U}}^{(i)}+\left(\bm{\Delta}^{(i+1)}-{\bm{W}}^{(i+1)}\right)$.

10: Check convergence. Set tolerances  $\displaystyle\tau_{pri}=$
$\displaystyle
mp\,\tau_{abs}+\tau_{rel}\,\max(\|{\bm{\Delta}}^{(i+1)}\|_{F},\|{\bm{W}}^{(i+1)}\|_{F})$
$\displaystyle\tau_{dual}=$ $\displaystyle
mp\,\tau_{abs}+\tau_{rel}\,\|{\bm{U}}^{(i+1)}\|_{F}/\rho^{(i)}\,.$ Define
$e_{p}=\|{\bm{\Delta}}^{(i+1)}-{\bm{W}}^{(i+1)}\|_{F}$ and
$e_{d}=\rho^{(i)}\|{\bm{W}}^{(i+1)}-{\bm{W}}^{(i)}\|_{F}$. If
$(e_{p}\leq\tau_{pri})\;AND\;(e_{d}\leq\tau_{dual})$, set converged = TRUE .

11: Update penalty parameter $\rho$ : If $e_{p}>\mu e_{d}$, set
$\rho^{(i+1)}=2\rho^{(i)}$, else if $e_{d}>\mu e_{p}$, set
$\rho^{(i+1)}=\rho^{(i)}/2$, otherwise $\rho^{(i+1)}=\rho^{(i)}$. We also need
to set ${\bm{U}}^{(i+1)}={\bm{U}}^{(i+1)}/2$ for $e_{p}>\mu e_{d}$ and
${\bm{U}}^{(i+1)}=2{\bm{U}}^{(i+1)}$ for $e_{d}>\mu e_{p}$.

12: $i\leftarrow i+1$

13:end while

14:Set $\hat{\bm{\Delta}}_{sym}=\frac{1}{2}({\bm{W}}+{\bm{W}}^{\top})$. If
$\|\hat{\bm{\Delta}}_{sym}^{(jk)}\|_{F}>0$, assign edge $\\{j,k\\}\in\hat{\cal
E}_{\Delta}$, else $\\{j,k\\}\not\in\hat{\cal E}_{\Delta}$.

Algorithm 2 ADMM Algorithm SA-ADMM

Input: As in Algorithm 1.
Output: Estimated $\hat{\bm{\Delta}}_{sym}$ and $\hat{\cal E}_{\Delta}$.

1:Follow lines 1-7 of Algorithm 1

2:With $(a)_{+}:=\max(0,a)$,
${\bm{A}}=({\bm{\Delta}}^{(i+1)}+{\bm{U}}^{(i)})^{(k\ell)}$ and
$k,\ell\in[mp]$, update $[{\bm{W}}]_{k\ell}\in\mathbb{R}$ as
$\displaystyle[{\bm{W}}^{(i+1)}]_{k\ell}=$
$\displaystyle\Big{(}1-\frac{(\lambda/\rho)}{\big{|}[{\bm{A}}]_{k\ell}\big{|}}\Big{)}_{+}[{\bm{A}}]_{k\ell}\,.$

3:Follow lines 9-14 of Algorithm 1

A pseudocode for the ADMM algorithm, MA-ADMM, used in this paper is given in
Algorithm 1 where we use the stopping (convergence) criterion following [38,
Sec. 3.3.1] and varying penalty parameter $\rho$ following [38, Sec. 3.4.1].
The stopping criterion is based on primal and dual residuals being small
where, in our case, at $(i+1)$st iteration, the primal residual is given by
$\bm{\Delta}^{(i+1)}-{\bm{W}}^{(i+1)}$ and the dual residual by
$\rho^{(i)}({\bm{W}}^{(i+1)}-{\bm{W}}^{(i)})$. Convergence criterion is met
when the norms of these residuals are below primary and dual tolerances
$\tau_{pri}$ and $\tau_{dual}$, respectively; see line 10 of Algorithm 1. In
turn, $\tau_{pri}$ and $\tau_{dual}$ are chosen using an absolute and relative
criterion as in line 10 of Algorithm 1 where $\tau_{abs}$ and $\tau_{rel}$ are
user chosen absolute and relative tolerances, respectively. Line 10 of
Algorithm 1 follows typical choices given in [38, Sec. 3.4.1]. For all
numerical results presented later, we used $\rho_{0}=2$, $\mu=10$, and
$\tau_{abs}=\tau_{rel}=10^{-4}$.

We will compare our approach with three other approaches in Sec. V-A. One of
them is the single attribute (SA) based ADMM approach (see [8, 28]. A
pseudocode of our implementation of this approach, SA-ADMM, is in Algorithm 2
which differs from in Algorithm 1 only in line 8 where we replace group lasso
with elementwise lasso.

###  III-B Proximal Gradient Descent Approach

It is a first-order method that is based on objective function values and
gradient evaluations. A pseudocode of the PGD method of [11], MA-proximal, is
in Algorithm 3, and that of [9], SA-proximal, is in Algorithm 4. Algorithm 4
differs from in Algorithm 3 only in line 8 where we replace group lasso with
element-wise lasso. For all numerical results presented later, we used
$\epsilon=10^{-3}$ in line 7.

Algorithm 3 PGD Algorithm MA-PGD

Input: Data $\\{{\bm{x}}(t)\\}_{t=1}^{n_{x}}$, ${\bm{x}}\in\mathbb{R}^{mp}$,
and $\\{{\bm{y}}(t)\\}_{t=1}^{n_{y}}$, ${\bm{y}}\in\mathbb{R}^{mp}$, tolerance
$\epsilon$, maximum number of iterations $i_{max}$
Output: Estimated $\hat{\bm{\Delta}}_{sym}$ and $\hat{\cal E}_{\Delta}$.

1:Calculate sample covariances
$\hat{\bm{\Sigma}}_{x}=\frac{1}{n_{x}}\sum_{t=1}^{n_{x}}{\bm{x}}(t){\bm{x}}^{\top}(t)$
and
$\hat{\bm{\Sigma}}_{y}=\frac{1}{n_{y}}\sum_{t=1}^{n_{y}}{\bm{y}}(t){\bm{y}}^{\top}(t)$.

2:Set
$\eta=1/\big{(}\phi_{max}(\hat{\bm{\Sigma}}_{x})\,\phi_{max}(\hat{\bm{\Sigma}}_{x})\big{)}$.
Initialize: ${\bm{\Delta}}^{(0)}={\bm{0}}$.

3:converged = FALSE, $i=0$

4:while converged = FALSE  AND  $i\leq i_{max}$, do

5: Set
${\bm{A}}={\bm{\Delta}}^{(i)}-\eta\Big{(}\hat{\bm{\Sigma}}_{x}{\bm{\Delta}}^{(i)}\hat{\bm{\Sigma}}_{y}-(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})\Big{)}$.

6: For $k,\ell\in[p]$, update $m\times m$ subblocks as
$\displaystyle({\bm{\Delta}}^{(i+1)})^{(k\ell)}=$
$\displaystyle\Big{(}1-\frac{\lambda\eta}{\|{\bm{A}}\|_{F}}\Big{)}_{+}{\bm{A}}^{(k\ell)}\,.$

7: If
$\frac{L_{\lambda}({\bm{\Delta}}^{(i+1)},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})-L_{\lambda}({\bm{\Delta}}^{(i)},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})}{L_{\lambda}({\bm{\Delta}}^{(i)},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})}\leq\epsilon$,
set converged = TRUE .

8: $i\leftarrow i+1$

9:end while

10:Set
$\hat{\bm{\Delta}}_{sym}=\frac{1}{2}({\bm{\Delta}}+{\bm{\Delta}}^{\top})$. If
$\|\hat{\bm{\Delta}}_{sym}^{(jk)}\|_{F}>0$, assign edge $\\{j,k\\}\in\hat{\cal
E}_{\Delta}$, else $\\{j,k\\}\not\in\hat{\cal E}_{\Delta}$.

Algorithm 4 PGD Algorithm SA-PGD

Input: As in Algorithm 3
Output: Estimated $\hat{\bm{\Delta}}$ and $\hat{\cal E}_{\Delta}$

1:Follow lines 1-5 of Algorithm 3

2:For $k,\ell\in[mp]$, update  $\displaystyle[{\bm{\Delta}}^{(i+1)}]_{k\ell}=$
$\displaystyle\Big{(}1-\frac{\lambda\eta}{\big{|}[{\bm{A}}]_{k\ell}\big{|}}\Big{)}_{+}[{\bm{A}}]_{k\ell}\,.$

3:Follow lines 7-10 of Algorithm 3

###  III-C Computational Complexity

The computational complexity of ADMM and PGD methods has been discussed in [9]
for SA differential graphs, and it is of the same order for MA graphs, because
the difference lies only in lasso versus group lasso, i.e., element-wise soft-
thresholding versus group-wise soft-thresholding. Noting that we have
$mp\times mp$ precision matrices, by [9], the computational complexity of the
ADMM approaches (our proposed and that of [8, 28]) is ${\cal O}((mp)^{3})$
while that of the PGD methods of [9, 11] is either ${\cal O}((mp)^{3})$ when
as implemented in Algorithms 3 and 4, or ${\cal O}((n_{x}+n_{y})(mp)^{2})$
when an alternative implementation of the cost gradient in line 5 of
Algorithms 3 and 4 is used (see [9, Sec. 2.2]). For $n_{x}+n_{y}\geq mp$,
there is no advantage to this alternative approach.

###  III-D Convergence of ADMM

The objective function (4), is strictly convex. It is also closed, proper and
lower semi-continuous. Hence, for any fixed $\rho>0$, the (two-block) ADMM
algorithm is guaranteed to converge [38, Sec. 3.2], in the sense that we have
primal residual convergence to 0, dual residual convergence to 0, and
objective function convergence to the optimal value. For varying $\rho$, the
convergence of ADMM has not been proven, but if we additionally impose
$\rho^{(i)}=\rho^{(i_{0})}>0$ for $i\geq i_{0}$ for some $i_{0}$, we have
convergence [38, Sec. 3.4.1].

###  III-E Model Selection

Following the lasso penalty work of [8] (who invokes [12]), we use the
following criterion for selection of group lasso penalty:

$\displaystyle BIC(\lambda)=$
$\displaystyle(n_{x}+n_{y})\,\|\hat{\bm{\Sigma}}_{x}\hat{\bm{\Delta}}\hat{\bm{\Sigma}}_{y}-(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})\|_{F}$
$\displaystyle\quad+\ln(n_{x}+n_{y})\,|\hat{\bm{\Delta}}|_{0}$ (14)

where $|{\bm{A}}|_{0}$ denotes number of nonzero elements in ${\bm{A}}$ and
$\hat{\bm{\Delta}}$ obeys (5). Choose $\lambda$ to minimize $BIC(\lambda)$.
Following [8] we term it BIC (Bayesian information criterion) even though the
cost function used is not negative log-likelihood although
$\ln(n_{x}+n_{y})\,|\hat{\bm{\Delta}}|_{0}$ penalizes over-parametrization as
in BIC. It is based on the fact that true ${\bm{\Delta}}^{\ast}$ satisfies
${\bm{\Sigma}}_{x}^{\ast}{\bm{\Delta}}^{\ast}{\bm{\Sigma}}_{y}^{\ast}-({\bm{\Sigma}}_{x}^{\ast}-{\bm{\Sigma}}_{y}^{\ast})={\bm{0}}$.
Since (14) is not scale invariant, we scale both $\hat{\bm{\Sigma}}_{x}$ and
$\hat{\bm{\Sigma}}_{y}$ (and $\hat{\bm{\Delta}}$ commensurately) by
$\bar{\bm{\Sigma}}^{-1}$ where
$\bar{\bm{\Sigma}}=\mbox{diag}\\{\hat{\bm{\Sigma}}_{x}\\}$ is a diagonal
matrix of diagonal elements of $\hat{\bm{\Sigma}}_{x}$.

In our simulations we search over $\lambda\in[\lambda_{\ell},\lambda_{u}]$,
where $\lambda_{\ell}$ and $\lambda_{u}$ are selected via a heuristic as in
[17]. Find the smallest $\lambda$, labeled $\lambda_{sm}$ for which we get a
no-edge model; then we set $\lambda_{u}=\lambda_{sm}/2$ and
$\lambda_{\ell}=\lambda_{u}/10$.

##  IV Theoretical Analysis

Here we analyze the properties of $\hat{\bm{\Delta}}$. Theorem 1 follows the
approach(es) of [31, 16, 29, 8, 28] while Theorem 2 follows the general
framework of [30]. The general method of [31] used in [16, 29, 8, 28] requires
an irrepresentability condition (see (20)) which is not needed by the method
of [30]. On the other hand, our Theorem 2 does not have a result like Theorem
1(ii), the oracle property, or support recovery Theorem 1(iv).

First some notation. Define the true differential edgeset

$\displaystyle S=$ $\displaystyle{\cal
E}_{\Delta^{\ast}}=\\{\\{k,\ell\\}\,:\,\|{\bm{\Delta}}^{\ast(k\ell)}\|_{F}\neq
0\\}\,,\quad s=|S|\,.$ (15)

Define

${\bm{\Gamma}}^{\ast}={\bm{\Sigma}}_{y}^{\ast}\boxtimes{\bm{\Sigma}}_{x}^{\ast}\,,\quad\hat{\bm{\Gamma}}=\hat{\bm{\Sigma}}_{y}\boxtimes\hat{\bm{\Sigma}}_{x}\,.$
(16)

Also, recall the operator $\bm{\mathcal{C}}(\cdot)$ defined in Sec. I-C. In
the rest of this section, we allow $p$, $s$ and $\lambda$ to be a functions of
sample size $n$, denoted as $p_{n}$, $s_{n}$ and $\lambda_{n}$, respectively.
Define

$\displaystyle M$
$\displaystyle=\max\\{\|\bm{\mathcal{C}}({\bm{\Sigma}}_{x}^{\ast})\|_{\infty}\,,\|\bm{\mathcal{C}}({\bm{\Sigma}}_{y}^{\ast})\|_{\infty}\\}\,,$
(17) $\displaystyle M_{\Sigma}$
$\displaystyle=\max\\{\|\bm{\mathcal{C}}({\bm{\Sigma}}_{x}^{\ast})\|_{1,\infty}\,,\|\bm{\mathcal{C}}({\bm{\Sigma}}_{y}^{\ast})\|_{1,\infty}\\}\,,$
(18) $\displaystyle\kappa_{\Gamma}$
$\displaystyle=\|\bm{\mathcal{C}}(\Gamma_{S,S}^{\ast})^{-1}\|_{1,\infty}\,,$
(19) $\displaystyle\alpha$ $\displaystyle=1-\max_{e\in
S^{c}}\|\bm{\mathcal{C}}({\bm{\Gamma}}_{e,S}^{\ast}({\bm{\Gamma}}_{S,S}^{\ast})^{-1})\|_{1}\,,$
(20) $\displaystyle\bar{\sigma}_{xy}$
$\displaystyle=\max\\{\max_{i}[{\bm{\Sigma}}_{x}^{\ast}]_{ii},\,\max_{i}[{\bm{\Sigma}}_{y}^{\ast}]_{ii}\\}$
(21) $\displaystyle{C}_{0}$
$\displaystyle=40\,m\,\bar{\sigma}_{xy}\sqrt{2\big{(}\tau+\ln(4m^{2})/\ln(p_{n})\big{)}}$
(22)

where $S$ and ${\bm{\Gamma}}^{\ast}$ have been defined in (15) and (16). In
(20), we require $0<\alpha<1$, and the expression

$\max_{e\in
S^{c}}\|\bm{\mathcal{C}}({\bm{\Gamma}}_{e,S}^{\ast}({\bm{\Gamma}}_{S,S}^{\ast})^{-1})\|_{1}\leq
1-\alpha$

for some $\alpha\in(0,1)$ is called the irrepresentability condition. Similar
conditions are also used in [31, 16, 29, 8, 28].

Let $\hat{\bm{\Delta}}$ be as in (5).
Theorem 1. For the system model of Sec. II, under (15) and the
irrepresentability condition (20) for some $\alpha\in(0,1)$, if

$\displaystyle\lambda_{n}=\max\Big{\\{}\frac{8}{\alpha}\,,\frac{3}{\alpha\bar{C}_{\alpha}}\,s_{n}\kappa_{\Gamma}MC_{M\kappa}\Big{\\}}{C}_{0}\sqrt{\frac{\ln(p_{n})}{n}}$
(23) $\displaystyle
n=\min(n_{x},n_{y})>\,\max\Big{\\{}\frac{1}{\min\\{M^{2},1\\}}\,,81M^{2}s_{n}^{2}\kappa_{\Gamma}^{2},$
$\displaystyle\quad\quad\quad\quad\frac{9s_{n}^{2}}{(\alpha\bar{C}_{\alpha})^{2}}(\kappa_{\Gamma}MC_{M\kappa})^{2}\Big{\\}}{C}_{0}^{2}\ln(p_{n})$
(24)

where $\bar{C}_{\alpha}=\frac{1-\alpha}{2(2M+1)-2\alpha M}$ and
$C_{M\kappa}=1.5\big{(}1+\kappa_{\Gamma}\min\\{s_{n}M^{2},M_{\Sigma}^{2}\\}\big{)}$,
then with probability $>1-2/p_{n}^{\tau-2}$, for any $\tau>2$, we have

  * (i)

$\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{\infty}\leq(C_{b1}+C_{b2}){C}_{0}\sqrt{\frac{\ln(p_{n})}{n}}$

$\displaystyle\mbox{where }\;C_{b1}=$ $\displaystyle
3\kappa_{\Gamma}\,\max\big{\\{}\frac{8}{\alpha}\,,\frac{3}{\alpha\bar{C}_{\alpha}}\,s_{n}\kappa_{\Gamma}MC_{M\kappa}\big{\\}}$
$\displaystyle C_{b2}=$ $\displaystyle 9s_{n}\kappa_{\Gamma}^{2}M^{2}\,.$

  * (ii)

$\hat{\bm{\Delta}}_{S^{c}}={\bm{0}}$.

  * (iii)

$\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{F}\leq\sqrt{s_{n}}\,\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{\infty}$
.

  * (iv)

Additionally, if $\min_{(k,\ell)\in
S}\|({\bm{\Delta}}^{\ast})^{(k\ell)}\|_{F}\geq\newline
2(C_{b1}+C_{b2}){C}_{0}\sqrt{\frac{\ln(p_{n})}{n}}$, then $P({\cal
G}_{\hat{\Delta}}={\cal G}_{{\Delta}^{\ast}})>1-2/p_{n}^{\tau-2}$ (support
recovery) $\quad\bullet$

The proof of Theorem 1 is given in Appendix A.

Now we present Theorem 2 that follows the general framework of [30]. Let
$\hat{\bm{\Delta}}$ be as in (5).
Theorem 2. For the system model of Sec. II, under (15), if

$\displaystyle\lambda_{n}\geq(4+6Ms_{n}C_{1}){C}_{0}\sqrt{\frac{\ln(p_{n})}{n}}$
(25) $\displaystyle
n=\min(n_{x},n_{y})>\,\max\Big{\\{}\frac{1}{M^{2}}\,,\big{(}\frac{96Ms_{n}}{\phi^{\ast}_{min}}\big{)}^{2}\Big{\\}}{C}_{0}^{2}\ln(p_{n})$
(26)

where $C_{1}=\max_{\\{k,\ell\\}\in V\times
V}\|({\bm{\Delta}}^{\ast})^{(k\ell)}\|_{F}$ and
$\phi^{\ast}_{min}=\phi_{min}({\bm{\Sigma}}_{x}^{\ast})\phi_{min}({\bm{\Sigma}}_{y}^{\ast})$,
then with probability $>1-2/p_{n}^{\tau-2}$, for any $\tau>2$, we have

$\displaystyle\|\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast}\|_{F}\leq\frac{12\sqrt{s_{n}}}{\phi^{\ast}_{min}}(4+6Ms_{n}C_{1}){C}_{0}\sqrt{\frac{\ln(p_{n})}{n}}\quad\bullet$
(27)

The proof of Theorem 2 is given in Appendix B. Note that
$\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{F}=\|\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast}\|_{F}$
when comparing Theorems 1 and 2.
Remark 1: Convergence Rate. If $M$, $M_{\Sigma}$ and $\kappa_{\Gamma}$ stay
bounded with increasing sample size $n$, we have
$\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{F}={\cal
O}_{P}(s_{n}^{1.5}\sqrt{\ln(p_{n})/n})$. Therefore, for
$\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{F}\rightarrow
0$ as $n\rightarrow\infty$, we must have
$s_{n}^{1.5}\sqrt{\ln(p_{n})/n}\rightarrow 0$. The SA results in [8] need
$s_{n}^{3.5}\sqrt{\ln(p_{n})/n}\rightarrow 0$ when we take into account the
dependence of various constants on $s_{n}$ in [8]. Notice that $M_{\Sigma}$
constraints covariances ${\bm{\Sigma}}_{x}^{\ast}$ and
${\bm{\Sigma}}_{y}^{\ast}$ which can be dense even if
${\bm{\Omega}}_{x}^{\ast}$ and ${\bm{\Omega}}_{y}^{\ast}$ are sparse (they
need not be sparse for differential estimation), making them possibly
unbounded with increasing sample size $n$. In this case we use
$\min\\{s_{n}M^{2},M_{\Sigma}^{2}\\}=s_{n}M^{2}$ in $C_{M\kappa}$ and
$C_{b1}$, with $M$ bounded, leading to
$\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{F}={\cal
O}_{P}(s_{n}^{2.5}\sqrt{\ln(p_{n})/n})$. On the other hand, in Theorem 2, we
always have
$\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{F}={\cal
O}_{P}(s_{n}^{1.5}\sqrt{\ln(p_{n})/n})$. $\;\;\Box$

Remark 2: Our results assume Gaussian data. Theorems 1 and 2 continue to hold
for more general sub-Gaussian distributions (Gaussian distribution is a sub-
Gaussian distribution) except that the zeros in the precision matrix (see
(1)), or in the difference of two precision matrices, no longer signify
conditional independence (or change in conditional independence) of the random
vectors associated with the respective nodes; they only imply zero partial
correlation. We state Lemma 1 in Appendix A for sub-Gaussian distributions,
following [31, Lemma 1]. Lemma 2 is then specialized to Gaussian distributions
by setting the sub-Gaussian parameter $\sigma_{sg}=1$. If $\sigma_{sg}\neq 1$,
then the only changes required in Theorems 1 and 2 is a scaling of $C_{0}$ in
(22) where one needs to replace the factor of 40 with $8(1+4\sigma_{sg}^{2})$.
$\;\;\Box$

Remark 3: Theorems 1 and 2 assume a constant number of attributes $m$, with
only $p$, $s$ and $\lambda$ allowed to be functions of sample size $n$, and
our Remark 1 reflects this fact. In terms of $m$, the bounds in Theorem 1(i)
and Theorem 2 are ${\cal O}(m^{3})$, which follows from $M={\cal O}(m)$ and
$C_{M\kappa}={\cal O}(m^{2})$ in Theorem 1, and $M={\cal O}(m)$, $C_{1}={\cal
O}(m)$ and $C_{0}={\cal O}(m)$ in Theorem 2. Therefore, for large $m$, one
would need much higher number of samples $n$, and it is not clear how to
circumvent this bottleneck. A different class of models based on matrix-valued
graphical modeling [40, 41, 42] is a potential solution. In a matrix graph
set-up, one has matrix-valued observations ${\bm{Z}}$, which in our context
would require attributes (components of ${\bm{z}}_{i}$) to be arranged along
rows and nodes $i$ along columns, with the covariance of
$\mbox{vec}({\bm{Z}})$ having a Kronecker-product structure. This structure
drastically reduces the number of unknowns from ${\cal O}((mp)^{2})$ to ${\cal
O}(m^{2}+p^{2})$. Prior reported work ([40, 41, 42]) is on matrix graph
estimation, with no reported work on differential matrix graphs. $\;\;\Box$

##  V Numerical Examples

We now present numerical results for both synthetic and real data to
illustrate the proposed approach. In synthetic data examples the ground truth
is known and this allows for assessment of the efficacy of various approaches.
In real data examples where the ground truth is unknown, our goal is
visualization and exploration of the differential conditional dependency
structures underlying the data.

###  V-A Synthetic Data: Erdös-Rènyi and Barabási-Albert Graphs

We consider two types of graphs: Erdös-Rènyi (ER) graph and Barabási-Albert
(BA) graph [43, 44]. The BA graphs are an example of scale-free graphs with
power law degree distribution [43]. In the ER graph, $p=100$ nodes are
connected to each other with probability $p_{er}=0.5$ and there are $m=3$
attributes per node whereas in the BA graph, we used $p=100$ and mean degree
of 2 to generate a BA graph using the procedure given in [44] (MATLAB function
BAmodel.m from https://github.com/ShanLu1984/Scale-Free-Network-Generation-
and-Comparison). In the upper triangular $\bm{\Omega}_{x}$, we set
$[\bm{\Omega}_{x}^{(jk)}]_{st}=0.5^{|s-t|}$ for $j=k=1,\cdots,p$,
$s,t=1,\cdots,m$. For $j\neq k$, if the two nodes are not connected in the
graph (ER or BA), we have $\bm{\Omega}^{(jk)}={\bm{0}}$, and if nodes $j$ and
$k$ are connected, then $[\bm{\Omega}^{(jk)}]_{st}$ is uniformly distributed
over $[-0.4,-0.1]\cup[0.1,0.4]$. Then add lower triangular elements to make
$\bm{\Omega}_{x}$ a symmetric matrix. To generate $\bm{\Omega}_{y}$, we follow
[8] and first generate a differential graph with
${\bm{\Delta}}\in\mathbb{R}^{(mp)\times(mp)}$ as an ER graph (regardless of
whether $\bm{\Omega}_{x}$ is based on ER or BA model), with connection
probability $p_{er}=0.05$ (sparse): if nodes $j$ and $k$ are connected in the
$\bm{\Omega}_{x}$ model, then each of $m^{2}$ elements of $\bm{\Delta}^{(jk)}$
is independently set to $\pm 0.9$ with equal probabilities. Then
$\bm{\Omega}_{y}=\bm{\Omega}_{x}+{\bm{\Delta}}$. Finally add $\gamma{\bm{I}}$
to $\bm{\Omega}_{y}$ and to $\bm{\Omega}_{x}$ and pick $\gamma$ so that
$\bm{\Omega}_{y}$ and $\bm{\Omega}_{x}$ are both positive definite. With
$\bm{\Phi}_{x}\bm{\Phi}_{x}^{\top}=\bm{\Omega}_{x}^{-1}$, we generate
${\bm{x}}=\bm{\Phi}{\bm{w}}$ with ${\bm{w}}\in\mathbb{R}^{mp}$ as zero-mean
Gaussian, with identity covariance, and similarly for ${\bm{y}}$. We generate
$n=n_{x}=n_{y}$ i.i.d. observations for ${\bm{x}}$ and ${\bm{y}}$, with $m=3$,
$p=100$, $n\in\\{100,200,300,400,800,1200,1600\\}$.

Figure 1: ROC curves for ER graph based on ADMM approaches. TPR=true positive rate, TNR=true negative rate Figure 2: ROC curves for BA graph based on ADMM approaches. TPR=true positive rate, TNR=true negative rate Approach |  $F_{1}$ score ($\pm\sigma$) | timing (s) ($\pm\sigma$) | TPR ($\pm\sigma$) | 1-TNR ($\pm\sigma$)
---|---|---|---|---
$n$ = 300
MA-ADMM | 0.6152 $\pm$0.0705 | 2.5044 $\pm$0.2939 | 0.6067 $\pm$0.1230 | 0.0184 $\pm$0.0080
MA-proximal | 0.6686 $\pm$0.0639 | 6.6931 $\pm$0.3743 | 0.6845 $\pm$0.1253 | 0.0184 $\pm$0.0085
SA-ADMM | 0.4549 $\pm$0.0332 | 0.1739 $\pm$0.0173 | 0.5795 $\pm$0.1132 | 0.0506 $\pm$0.0186
SA-proximal | 0.4772 $\pm$0.0328 | 0.3517 $\pm$0.0195 | 0.6263 $\pm$0.1157 | 0.0524 $\pm$0.0193
$n$ = 800
MA-ADMM | 0.8537 $\pm$0.0491 | 2.1911 $\pm$0.0639 | 0.9037 $\pm$0.0703 | 0.0111 $\pm$0.0041
MA-proximal | 0.8898 $\pm$0.0408 | 4.8277 $\pm$0.2647 | 0.9526 $\pm$0.0537 | 0.0009 $\pm$0.0041
SA-ADMM | 0.6336 $\pm$0.0055 | 0.1510 $\pm$0.0081 | 0.7468 $\pm$0.1017 | 0.0316 $\pm$0.0090
SA-proximal | 0.6612 $\pm$0.0401 | 0.2533 $\pm$0.0143 | 0.7917 $\pm$0.0931 | 0.0314 $\pm$0.0090
$n$ = 3000
MA-ADMM | 0.9795 $\pm$0.0152 | 2.1436 $\pm$0.0624 | 0.9928 $\pm$0.0018 | 0.0012 $\pm$0.0041
MA-proximal | 0.9914 $\pm$0.0106 | 3.8312 $\pm$0.3948 | 0.9964 $\pm$0.0121 | 0.0007 $\pm$0.0007
$n$ = 6000
MA-ADMM | 0.9914 $\pm$0.0087 | 1.9459 $\pm$0.0544 | 0.9997 $\pm$0.0013 | 0.0008 $\pm$0.0009
MA-proximal | 0.9979 $\pm$0.0037 | 3.517 $\pm$0.1829 | 0.9998 $\pm$0.0011 | 0.0002 $\pm$0.0004
TABLE I: Comparisons among various approaches: Erdös-Rènyi graph,
$n=300,\,800,\,3000,\,6000$, $p=100$, $m=3$. Tuning parameter $\lambda$ picked
to yield the highest $F_{1}$ score. Results based on 100 runs.

Simulation results based on 100 runs are shown in Figs. 1-4. By changing the
penalty parameter $\lambda$ and determining the resulting edges, we calculated
the true positive rate (TPR) and false positive rate 1-TNR (where TNR is the
true negative rate) over 100 runs. The receiver operating characteristic (ROC)
for ER graphs is shown in Fig. 1 for our MA-ADMM approach (labeled “MA”) as
well as for a SA-ADMM approach (labeled “SA”), based on [28], where we first
estimate an $mp$-node differential graph, and then use
$\|\hat{\bm{\Delta}}^{(k\ell)}\|_{F}\neq 0$ $\Leftrightarrow$
$\\{\\{k,\ell\\}\in{\cal E}_{\Delta}$. It is seen from Fig. 1 that our MA-ADMM
approach outperforms the SA-ADMM approach (that uses the same cost but
element-wise lasso penalty instead of group-lasso penalty). Fig. 2 is the
counterpart of Fig. 1 for BA graphs., and comments made regarding Fig. 1 apply
here too.

Figure 3: ROC curves for ER graph based on ADMM as well as proximal
approaches.

In Fig. 3 we compare ROC curves of our MA-ADMM approach with that for the MA-
proximal and SA-proximal approaches of [11] and [9], respectively, for $n=300$
and 800. It is seen that the MA-proximal approach outperforms our MA-ADMM
approach, while both significantly outperform the SA-proximal approach (and
the SA-ADMM approach whose ROC curves are in Fig. 1). In Table I, for the ER
graph ($n=300,800$, $p=100$, $m=3$) we compare the four approaches (MA-ADMM,
SA-ADMM, MA-proximal, SA-proximal) in terms of the $F_{1}$ score, execution
time (based on tic-toc functions in MATLAB), TPR and 1-TNR, for fixed penalty
parameter $\lambda$ selected from a grid of values (the same as for computing
the ROC curves) to maximize the $F_{1}$ score averaged over 100 runs. All
algorithms were run on a Window 10 Pro operating system with processor
Intel(R) Core(TM) i7-10700 CPU @2.90 GHz with 32 GB RAM, using MATLAB R2023a.
Notice that while the MA-proximal approach outperforms our MA-ADMM approach,
it also takes more than twice the computation time for the MA-ADMM approach.
Similarly, the SA-proximal approach takes more than twice the computation time
for the SA-ADMM approach. The latter observation is consistent with the
findings of [9].

In Table I, we also show results for $n=3000$ and $6000$ for MA-ADMM and MA-
proximal approaches in order to provide further empirical validation of the
theoretical results stated in Theorem 1. Theorem 1 states that for
sufficiently large sample size $n$, one can recover the differential graph
structure exactly w.h.p. It is seen that the $F_{1}$ score approaches one with
increasing $n$, implying graph support recovery, as claimed in Theorem 1(iv).
Note that Theorem 1(iv) holds w.h.p., not with probability one, implying a
nonzero probability of possibly inexact graph recovery, yielding an $F1$-score
less than 1 for $n=3000,6000$. The sample sizes of $n=300,800$ are not large
enough to yield an $F_{1}$-score close to 1. There is a lower bound (24) on
$n$ for Theorem 1(iv) to hold w.h.p. The sample sizes of $n=300,800$ are
apparently less than bound (which is not easily computable since it needs
$\alpha$ and $\kappa_{\Gamma}$).

In Fig. 4 we show the results based on 100 runs for our approach when BIC
parameter selection method (Sec. III-E) is applied in conjunction with the MA-
ADMM approach. Here we show the TPR, 1-TNR and $F_{1}$ score values along with
the $\pm\sigma$ error bars. The proposed approach works well both in terms of
$F_{1}$ score and TPR vs 1-TNR.

Figure 4: BIC based results for ER graph: $F_{1}$-scores, TPR and 1-TNR

(a) Proposed approach: $\|\hat{\bm{\Delta}}^{(k\ell)}\|_{F}$,
$k,\ell=1,2,\cdots,11$

(b) Proposed approach: edges
$\big{\\{}\\{k,\ell\\}\,:\,\|\hat{\bm{\Delta}}^{(k\ell)}\|_{F}>0\big{\\}}$

(c) $\|\hat{\bm{\Omega}}_{x}^{(k\ell)}\|_{F}$, $k,\ell=1,2,\cdots,11$ [17]

(d) $\|\hat{\bm{\Omega}}_{y}^{(k\ell)}\|_{F}$, $k,\ell=1,2,\cdots,11$ [17]

(e) $\|\hat{\bm{\Omega}}_{y}^{(k\ell)}-\hat{\bm{\Omega}}_{x}^{(k\ell)}\|_{F}$
([17])

(f) Estimated edges
$\big{\\{}\\{k,\ell\\}\,:\,\|\hat{\bm{\Omega}}_{y}^{(k\ell)}-\hat{\bm{\Omega}}_{x}^{(k\ell)}\|_{F}>0\big{\\}}$

Figure 5: Differential graphs comparing Beijing air-quality datasets [45] for
years 2013-14 and 2014-15: 8 monitoring stations and 11 features ($m=8$,
$p=11$, $n_{x}=n_{y}=365$). The features are numbered 1-11 beginning PM2.5
(PM2.5) and moving counter-clockwise in Fig. 5(b).

(a) Proposed approach: $\|\hat{\bm{\Delta}}^{(k\ell)}\|_{F}$,
$k,\ell=1,2,\cdots,11$

(b) Proposed approach: edges
$\big{\\{}\\{k,\ell\\}\,:\,\|\hat{\bm{\Delta}}^{(k\ell)}\|_{F}>0\big{\\}}$

(c) Estimated edges
$\big{\\{}\\{k,\ell\\}\,:\,\|\hat{\bm{\Omega}}_{x}^{(k\ell)}\|_{F}>0\big{\\}}$
[17]

(d) Estimated edges
$\big{\\{}\\{k,\ell\\}\,:\,\|\hat{\bm{\Omega}}_{y}^{(k\ell)}\|_{F}>0\big{\\}}$
[17]

(e) $\|\hat{\bm{\Omega}}_{y}^{(k\ell)}-\hat{\bm{\Omega}}_{x}^{(k\ell)}\|_{F}$
([17])

(f) Estimated edges
$\big{\\{}\\{k,\ell\\}\,:\,\|\hat{\bm{\Omega}}_{y}^{(k\ell)}-\hat{\bm{\Omega}}_{x}^{(k\ell)}\|_{F}>0\big{\\}}$

Figure 6: Differential graphs comparing Beijing air-quality datasets [45]
acquired from two sets of monitoring stations, 4 stations per set, year
2013-14: 4 monitoring stations and 11 features ($m=4$, $p=11$,
$n_{x}=n_{y}=365$).

###  V-B Real Data: Beijing air-quality dataset [45]

Here we consider Beijing air-quality dataset [45, 46], downloaded from
https://archive.ics.uci.edu/ml/datasets/Beijing+Multi-Site+Air-Quality+Data.
This data set includes hourly air pollutants data from 12 nationally-
controlled air-quality monitoring sites in the Beijing area from the Beijing
Municipal Environmental Monitoring Center, and meteorological data in each
air-quality site are matched with the nearest weather station from the China
Meteorological Administration. The time period is from March 1st, 2013 to
February 28th, 2017.The six air pollutants are PM2.5, PM10, SO2, NO2, CO, and
O3, and the meteorological data is comprised of five features: temperature,
atmospheric pressure, dew point, wind speed, and rain; we did not use wind
direction. Thus we have eleven features. We used data from 8 sites: 4
suburban/rural sites – Changping, Huairou, Shunyi, Dingling, and 4 urban area
stations – Aotizhongxin, Dongsi, Guanyuan, Gucheng [46, Fig. 1]. The data are
averaged over 24 hour period to yield daily averages. We used one year of
daily data resulting in $n_{x}=n_{y}=365$ days. The stations are used as
attributes, with $m=8$ for comparison between years 2013-14 and 2014-15, and
$m=4$ for comparison between suburban/rural sites and urban sites using
2013-14 year data.

We pre-process the data as follows. Given $i$th feature data
${\bm{z}}_{i}(t)\in\mathbb{R}^{m}$, we transform it to
$\bar{\bm{z}}_{i}(t)=\ln({\bm{z}}_{i}(t)/{\bm{z}}_{i}(t-1))$ and then detrend
it (i.e., remove the best straight-line fit using the MATLAB function
detrend). Finally, we scale the detrended scalar sequence to have a mean-
square value of one over $n_{x}$ or $n_{y}$ samples. The logarithmic
transformation and detrending of each feature sequence makes the sequence
closer to (univariate) stationary and Gaussian, while scaling “balances” the
possible wide variations in the scale of various feature measurements. All
temperatures were converted from Celsius to Kelvin to avoid negative numbers,
and if a value of a feature is zero (e.g., wind speed), we added a small
positive number to it, so that the logarithmic transformation is well-defined.

Fig. 5 shows the estimated differential graphs when comparing daily-averaged
data from 2013-14 ($x$-data) to that from 2014-15 ($y$-data), with air-quality
and meteorological variables as $p=11$ features measured at 8 monitoring sites
($m$=8). The objective is to visualize and explore differential conditional
dependency relationships among the 11 variables, comparing one year to
another, to investigate if pollution reduction measures have had any impact.
Our intuition is that one does not expect such rapid changes within a short
period of one year (see also [45]), therefore, our method should confirm our
intuition. Figs. 5(a)-(b) show estimated $\|\hat{\bm{\Delta}}^{(k\ell)}\|_{F}$
for various edges $\\{k,\ell\\}$, where it is unscaled in Fig. 5(a) but scaled
in Fig. 5(b) so that the largest $\|\hat{\bm{\Delta}}^{(k\ell)}\|_{F}$
(including $k=\ell$) is normalized to one. It is seen that differential graph
weights are essentially zero (very sparse), implying that there are no year-
to-year changes in the conditional dependency relationships among the 11
variables. This observation conforms to the findings of [45, 46]: no
significant year-to-year changes. We also estimated the MA graphs for each
year separately as $\|\hat{\bm{\Omega}}_{x}^{(k\ell)}\|_{F}$ and
$\|\hat{\bm{\Omega}}_{y}^{(k\ell)}\|_{F}$ (shown in Figs. 5(c)-(d)), using the
approach of [17], and based on the individual estimates, we computed the
differential graph
$\|\hat{\bm{\Omega}}_{y}^{(k\ell)}-\hat{\bm{\Omega}}_{x}^{(k\ell)}\|_{F}$
(shown in Figs. 5(e)-(f)). It is seen the separate-estimation approach does
not yield a sparse differential graph, even though the two individual graphs
in Figs. 5(c)-(d) are not all that different.

Fig. 6 shows the estimated differential graphs when comparing daily-averaged
data over the period 2013-14, from four suburban/rural sites ($x$-data) to
that from four urban sites ($y$-data), with air-quality and meteorological
variables as $p=11$ features measured at two sets of 4 monitoring sites
($m$=4). The objective again is to visualize and explore differential
conditional dependency relationships among the 11 variables, but in this case
comparing one subregion to another. There are significant differences in
meteorological conditions and pollutant sources, levels and mutual
interactions, among suburban and urban areas [45, 46]. The suburban areas
(located toward north) are less polluted than the urban areas (located toward
south) [45, 46]. Automobile exhaust is the main cause of NO2 which is likely
to undergo a chemical reaction with Ozone O3, thereby, lowering its
concentration [46]. Cold, dry air from the north reduces both dew point and
PM2.5 particle concentration in suburban areas while southerly wind brings
warmer and more humid air from the more polluted south that elevates the PM2.5
concentration [45]. The urban stations neighbor the south of Beijing which is
heavily installed with iron, steel and cement industries in Hebei province
[45]. Figs. 6(a)-(b) show estimated $\|\hat{\bm{\Delta}}^{(k\ell)}\|_{F}$ for
various edges $\\{k,\ell\\}$, where it is unscaled in Fig. 6(a) but scaled in
Fig. 6(b) so that the largest $\|\hat{\bm{\Delta}}^{(k\ell)}\|_{F}$ (including
$k=\ell$) is normalized to one. It is seen that quite a few of the
differential graph weights are significantly non-zero in Fig. 6(a), unlike
that in Fig. 5(a), implying significant differences in the conditional
dependency relationships among the 11 variables for suburban and urban areas.
This observation conforms to the findings of [45, 46]. The comments made
regarding Figs. 5(c)-(f) apply as well to Figs. 6(c)-(f).

##  VI Conclusions

A group lasso penalized D-trace loss function approach for differential graph
learning from multi-attribute data was presented. An ADMM algorithm was
presented to optimize the convex objective function. Theoretical analysis
establishing consistency of the estimator in high-dimensional settings was
performed. We tested the proposed approach on synthetic as well as real data.
In the synthetic data example, the multi-attribute approach is shown to
outperform a single-attribute approach in correctly detecting the differential
graph edges with ROC as the performance metric.

##  Appendix A Technical Lemmas and Proof of Theorem 1

In this Appendix, we provide a proof of Theorem 1. A necessary and sufficient
condition for minimization of convex
$L_{\lambda}({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})$ given
by (4) w.r.t. ${\bm{\Delta}}\in\mathbb{R}^{mp\times mp}$ is that
$\hat{\bm{\Delta}}$ minimizes (4) iff the zero matrix belongs to the sub-
differential of
$L_{\lambda}({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})$. That
is,

$\displaystyle{\bm{0}}=$ $\displaystyle\frac{\partial
L({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})}{\partial{\bm{\Delta}}}+\lambda{\bm{Z}}({\bm{\Delta}})\,\Big{|}_{{\bm{\Delta}}=\hat{\bm{\Delta}}}$
$\displaystyle=$
$\displaystyle\hat{\bm{\Sigma}}_{x}\hat{\bm{\Delta}}\hat{\bm{\Sigma}}_{y}-(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})+\lambda{\bm{Z}}(\hat{\bm{\Delta}})$
(28)

where
${\bm{Z}}({\bm{\Delta}})\in\partial\sum_{k,\ell=1}^{p}\|{\bm{\Delta}}^{(k\ell)}\|_{F}\in\mathbb{R}^{mp\times
mp}$, the sub-differential of group lasso penalty term, is given by [36, 37]

$\displaystyle({\bm{Z}}({\bm{\Delta}}))^{(k\ell)}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\frac{{\bm{\Delta}}^{(k\ell)}}{\|{\bm{\Delta}}^{(k\ell)}\|_{F}}&\mbox{if
}\,\|{\bm{\Delta}}^{(k\ell)}\|_{F}\neq 0\\\ {\bm{V}}\in\mathbb{R}^{m\times
m},\;\|{\bm{V}}\|_{F}\leq 1,&\mbox{if
}\,\|{\bm{\Delta}}^{(k\ell)}\|_{F}=0\end{array}\right.\,.$ (31)

In terms of $m\times m$ submatrices of ${\bm{\Delta}}$,
$\hat{\bm{\Sigma}}_{x}$, $\hat{\bm{\Sigma}}_{y}$ and ${\bm{Z}}({\bm{\Delta}})$
corresponding to various graph edges, using
$\mbox{bvec}({\bm{A}}{\bm{D}}{\bm{B}})=({\bm{B}}^{\top}\boxtimes{\bm{A}})\mbox{bvec}({\bm{D}})$
[33, Lemma 1], we may rewrite (28) as

$\displaystyle(\hat{\bm{\Sigma}}_{y}\boxtimes\hat{\bm{\Sigma}}_{x})\mbox{bvec}(\hat{\bm{\Delta}})-\mbox{bvec}(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})+\lambda\,\mbox{bvec}({\bm{Z}}(\hat{\bm{\Delta}}))={\bm{0}}$
(32)

Then (32) can be rewritten as

$\displaystyle\begin{bmatrix}\hat{\bm{\Gamma}}_{S,S}&\hat{\bm{\Gamma}}_{S,S^{c}}\\\
\hat{\bm{\Gamma}}_{S^{c},S}&\hat{\bm{\Gamma}}_{S^{c},S^{c}}\end{bmatrix}\begin{bmatrix}\mbox{bvec}(\hat{\bm{\Delta}}_{S})\\\
\mbox{bvec}(\hat{\bm{\Delta}}_{S^{c}})\end{bmatrix}-\begin{bmatrix}\mbox{bvec}((\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})_{S})\\\
\mbox{bvec}((\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})_{S^{c}})\end{bmatrix}$
$\displaystyle\quad+\lambda\,\begin{bmatrix}\mbox{bvec}({\bm{Z}}(\hat{\bm{\Delta}}_{S}))\\\
\mbox{bvec}({\bm{Z}}(\hat{\bm{\Delta}}_{S^{c}}))\end{bmatrix}=\begin{bmatrix}{\bm{0}}\\\
{\bm{0}}\end{bmatrix}\,.$ (33)

The general approach of [31] (followed in [16, 29, 8, 28]) is to first solve
the hypothetical constrained optimization problem with known edgeset $S$

$\tilde{\bm{\Delta}}=\arg\min_{{\bm{\Delta}}:{\bm{\Delta}}_{S^{c}}={\bm{0}}}L_{\lambda}({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})$
(34)

where $S^{c}$ is the complement of $S$. Since, by construction,
$\tilde{\bm{\Delta}}_{S^{c}}={\bm{0}}$, in this case (33) reduces to

$\displaystyle\hat{\bm{\Gamma}}_{S,S}\mbox{bvec}(\tilde{\bm{\Delta}}_{S})-\mbox{bvec}((\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})_{S})+\lambda\,\mbox{bvec}({\bm{Z}}(\tilde{\bm{\Delta}}_{S}))={\bm{0}}\,.$
(35)

In the approach of [31], one investigates conditions under which the solution
$\hat{\bm{\Delta}}$ to (4) is the same as the solution $\tilde{\bm{\Delta}}$
to (34). This is done by showing that $\hat{\bm{\Delta}}$ satisfies (33). The
choice $\hat{\bm{\Delta}}=\tilde{\bm{\Delta}}$ implies that
$\hat{\bm{\Delta}}_{S^{c}}={\bm{0}}$ and (35) is true with
$\tilde{\bm{\Delta}}$ replaced with $\hat{\bm{\Delta}}$. In order to satisfy
(33), it remains to show that for any edge $e\in S^{c}$, we have strict
feasibility

$\displaystyle\|\hat{\bm{\Gamma}}_{e,S}\mbox{bvec}(\tilde{\bm{\Delta}}_{S})-\mbox{bvec}((\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})_{e})\|_{2}<\lambda\,,$
(36)

where for ${\bm{a}}\in\mathbb{R}^{q}$,
$\|{\bm{a}}\|_{2}=\sqrt{{\bm{a}}^{\top}{\bm{a}}}$. This requires a set of
sufficient conditions stated in Theorem 1.

Lemma 1 follows from [31, Lemma 1]. It is stated for more general sub-Gaussian
distributions as in [31, Lemma 1], but will be used later for Gaussian
distributions, a subset of sub-Gaussian distributions.
Lemma 1. Suppose
$\hat{\bm{\Sigma}}=(1/n)\sum_{t=1}^{n}{\bm{x}}(t){\bm{x}}^{\top}(t)$, given
$n$ i.i.d. samples $\\{{\bm{x}}(t)\\}_{t=1}^{n}$ of zero-mean sub-Gaussian
${\bm{x}}\in\mathbb{R}^{mp}$ with covariance ${\bm{\Sigma}}^{\ast}$ such that
each component $x_{i}/\sqrt{\Sigma_{ii}^{\ast}}$ is sub-Gaussian with
parameter $\sigma_{sg}$. Define $\sigma_{max}=\max_{1\leq i\leq
mp_{n}}\Sigma_{ii}^{\ast}$ and

$\tilde{C}_{0}=8(1+4\sigma_{sg}^{2})m\sigma_{max}\sqrt{2\big{(}\tau+\ln(4m^{2})/\ln(p_{n})\big{)}}.$
(37)

Then

$P\Big{(}\|\bm{\mathcal{C}}(\hat{\bm{\Sigma}}-{\bm{\Sigma}}^{\ast})\|_{\infty}>\tilde{C}_{0}\sqrt{\ln(p_{n})/n}\Big{)}\leq
1/p_{n}^{\tau-2}$ (38)

for any $\tau>2$ and $n>2m^{2}(\ln(4m^{2})+\tau\ln(p_{n}))$. $\quad\bullet$
Proof. By [31, Lemma 1], with $b=8(1+4\sigma_{sg}^{2})$, we have

$P\Big{(}|[\hat{\bm{\Sigma}}-{\bm{\Sigma}}^{\ast}]_{ij}|>\delta\Big{)}\leq
4\,\exp(-c_{\ast}n\delta^{2})$ (39)

for any $\delta\in(0,b\,\sigma_{max})$ where
$c_{\ast}^{-1}=2b^{2}\sigma_{max}^{2}$. For any edge $\\{k,\ell\\}$ of the MA
graph, with $m^{2}$ edges $\\{i,j\\}$ of the corresponding SA graph associated
with $\\{k,\ell\\}$, using the union bound, we have

$\displaystyle P$
$\displaystyle\Big{(}|[\bm{\mathcal{C}}(\hat{\bm{\Sigma}}-{\bm{\Sigma}}^{\ast})]_{kl}|>\delta\Big{)}$
$\displaystyle\leq
P\Big{(}\max_{\\{i,j\\}\in\\{k,\ell\\}}([\hat{\bm{\Sigma}}-{\bm{\Sigma}}^{\ast}]_{ij})^{2}>\frac{\delta^{2}}{m^{2}}\Big{)}$
$\displaystyle\leq
m^{2}P\Big{(}|[\hat{\bm{\Sigma}}-{\bm{\Sigma}}^{\ast}]_{ij}|>\frac{\delta}{m}\Big{)}=4m^{2}\exp\big{(}-c_{\ast}n\frac{\delta^{2}}{m^{2}}\big{)}\,.$
(40)

Applying the union bound once more over all $p_{n}^{2}$ entries

$\displaystyle P$
$\displaystyle\Big{(}\|\bm{\mathcal{C}}(\hat{\bm{\Sigma}}-{\bm{\Sigma}}^{\ast})\|_{\infty}>\delta\Big{)}\leq
4(mp_{n})^{2}\exp\big{(}-c_{\ast}n\frac{\delta^{2}}{m^{2}}\big{)}=:P_{tb}\,.$
(41)

Choose
$\delta=\tilde{C}_{0}\sqrt{\ln(p_{n})/n}=bm\sigma_{max}\sqrt{2\ln(4p_{n}^{\tau}m^{2})/n}\,$.
Then we have

$\displaystyle P_{tb}=$ $\displaystyle
4(mp_{n})^{2}\exp\big{(}\ln(4p_{n}^{\tau}m^{2})^{-1}\big{)}=1/p_{n}^{\tau-2}$
(42)

provided $\delta\in(0,b\sigma_{max})$. Therefore, we need to have
$\tilde{C}_{0}\sqrt{\ln(p_{n})/n}<b\sigma_{max}$ requiring
$n>2m^{2}(\ln(4m^{2})+\tau\ln(p_{n}))$. This completes the proof. $\quad\Box$

Using the union bound, Lemma 1 and Gaussian assumption, we have Lemma 2.
Lemma 2. Let $\hat{\bm{\Sigma}}_{x}$ and $\hat{\bm{\Sigma}}_{y}$ be as in (2),
$\bar{\sigma}_{xy}$ as in (21), ${C}_{0}$ as in (22) and assume data are
Gaussian. Define $n=\min(n_{x},n_{y})$ and

$\displaystyle{\mathcal{A}}=$
$\displaystyle\max\big{\\{}\|\bm{\mathcal{C}}(\hat{{\bm{\Sigma}}_{x}}-{\bm{\Sigma}}_{x}^{\ast})\|_{\infty}\,,\|\bm{\mathcal{C}}(\hat{{\bm{\Sigma}}_{y}}-{\bm{\Sigma}}_{y}^{\ast})\|_{\infty}\big{\\}}\,.$
(43)

Then for any $\tau>2$ and $n>2m^{2}\ln(4m^{2}p_{n}^{\tau})$,

$\displaystyle P$
$\displaystyle\Big{(}{\mathcal{A}}>{C}_{0}\sqrt{\ln(p_{n})/n}\Big{)}\leq
2/p_{n}^{\tau-2}\quad\bullet$ (44)

Proof. For Gaussian distribution, the sub-Gaussian parameter $\sigma_{sg}$ of
Lemma 1 equals 1. Then $8(1+4\sigma_{sg}^{2})=40$. Let
$C_{0x}=40m(\max_{i}\Sigma_{x,ii}^{\ast})\sqrt{2\big{(}\tau+\ln(4m^{2})/\ln(p_{n})\big{)}}$
and
$C_{0y}=40m(\max_{i}\Sigma_{y,ii}^{\ast})\sqrt{2\big{(}\tau+\ln(4m^{2})/\ln(p_{n})\big{)}}$
(where $\Sigma_{x,ii}^{\ast}=[{\bm{\Sigma}}_{x}^{\ast})]_{ii}$, etc.). Then
using Lemma 1 and union bound,

$\displaystyle P$
$\displaystyle\Big{(}{\mathcal{A}}>{C}_{0}\sqrt{\ln(p_{n})/n}\Big{)}$
$\displaystyle\leq$ $\displaystyle
P\Big{(}\|\bm{\mathcal{C}}(\hat{\bm{\Sigma}}_{x}-{\bm{\Sigma}}_{x}^{\ast})\|_{\infty}>C_{0}\sqrt{\ln(p_{n})/n}\Big{)}$
$\displaystyle\;+P\Big{(}\|\bm{\mathcal{C}}(\hat{\bm{\Sigma}}_{y}-{\bm{\Sigma}}_{y}^{\ast})\|_{\infty}>C_{0}\sqrt{\ln(p_{n})/n}\Big{)}$
$\displaystyle\leq$ $\displaystyle 2/p_{n}^{\tau-2}$ (45)

since $C_{0}\geq C_{0x}$ and $C_{0}\geq C_{0y}$. $\quad\Box$

Recall (16)-(20) and define

$\displaystyle{\bm{\Delta}}_{x}=$
$\displaystyle\hat{\bm{\Sigma}}_{x}-{\bm{\Sigma}}_{x}^{\ast}\,,\;{\bm{\Delta}}_{y}=\hat{\bm{\Sigma}}_{y}-{\bm{\Sigma}}_{y}^{\ast}\,,\;{\bm{\Delta}}_{\Gamma}=\hat{\bm{\Gamma}}-{\bm{\Gamma}}^{\ast}\,,$
(46) $\displaystyle{\bm{\Delta}}_{\Sigma}=$
$\displaystyle{\bm{\Delta}}_{x}-{\bm{\Delta}}_{y}\,,\;{\epsilon}_{x}=\|\bm{\mathcal{C}}({\bm{\Delta}}_{x})\|_{\infty}\,,$
(47) $\displaystyle{\epsilon}_{y}=$
$\displaystyle\|\bm{\mathcal{C}}({\bm{\Delta}}_{y})\|_{\infty}\,,\;{\epsilon}>\max\\{{\epsilon}_{x},{\epsilon}_{y}\\}\,.$
(48)

Lemma 3. Assume that

$\displaystyle\kappa_{\Gamma}$
$\displaystyle<\frac{1}{3s_{n}(\epsilon^{2}+2M\epsilon)}\,.$ (49)

Let (${\bm{\Gamma}}_{S,S}^{-\ast}$ denotes
$({\bm{\Gamma}}_{S,S}^{\ast})^{-1}$)

$\displaystyle R({\bm{\Delta}}_{\Gamma})=$
$\displaystyle\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast}+{\bm{\Gamma}}_{S,S}^{-\ast}({\bm{\Delta}}_{\Gamma})_{S,S}{\bm{\Gamma}}_{S,S}^{-\ast}\,.$
(50)

Then we have

$\displaystyle\|\bm{\mathcal{C}}(R({\bm{\Delta}}_{\Gamma}))\|_{\infty}\leq\frac{3}{2}\kappa_{\Gamma}^{3}s_{n}(\epsilon^{2}+2M\epsilon)^{2}\,,$
(51)
$\displaystyle\|\bm{\mathcal{C}}(R({\bm{\Delta}}_{\Gamma}))\|_{1,\infty}\leq\frac{3}{2}\kappa_{\Gamma}^{3}s_{n}^{2}(\epsilon^{2}+2M\epsilon)^{2}\,,$
(52)
$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{\infty}$
$\displaystyle\;\;\leq\kappa_{\Gamma}^{2}(\epsilon^{2}+2M\epsilon)\big{(}1+1.5s_{n}\kappa_{\Gamma}(\epsilon^{2}+2M\epsilon)\big{)}\,,$
(53)
$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}\leq
s_{n}\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{\infty}\quad\bullet$
(54)

Proof. We have

$\displaystyle{\bm{\Delta}}_{\Gamma}=$
$\displaystyle\hat{\bm{\Sigma}}_{y}\boxtimes\hat{\bm{\Sigma}}_{x}-{\bm{\Sigma}}_{y}^{\ast}\boxtimes{\bm{\Sigma}}_{x}^{\ast}$
$\displaystyle=$
$\displaystyle{\bm{\Delta}}_{y}\boxtimes{\bm{\Delta}}_{x}+{\bm{\Sigma}}_{y}^{\ast}\boxtimes{\bm{\Delta}}_{x}+{\bm{\Delta}}_{y}\boxtimes{\bm{\Sigma}}_{x}^{\ast}.$
(55)

By [16, Lemma 14],

$\displaystyle\|\bm{\mathcal{C}}({\bm{A}}{\bm{B}})\|_{1,\infty}$
$\displaystyle\leq\|\bm{\mathcal{C}}({\bm{A}})\|_{1,\infty}\|\bm{\mathcal{C}}({\bm{B}})\|_{1,\infty}$
(56)

and by [16, Lemma 15],

$\displaystyle\|\bm{\mathcal{C}}({\bm{A}}{\bm{B}})\|_{\infty}$
$\displaystyle\leq\|\bm{\mathcal{C}}({\bm{A}})\|_{\infty}\|\bm{\mathcal{C}}({\bm{B}}^{\top})\|_{1,\infty}\,.$
(57)

Since $\|{\bm{A}}\otimes{\bm{B}}\|_{F}=\|{\bm{A}}\|_{F}\,\|{\bm{B}}\|_{F}$ and
${\bm{A}}\boxtimes{\bm{B}}=[{\bm{A}}_{ij}\boxtimes{\bm{B}}]_{ij}=[[{\bm{A}}_{ij}\otimes{\bm{B}}_{k\ell}]_{k\ell}]_{ij}$,
we have

$\displaystyle\|\bm{\mathcal{C}}({\bm{A}}\boxtimes{\bm{B}})\|_{\infty}$
$\displaystyle\leq\|\bm{\mathcal{C}}({\bm{A}})\|_{\infty}\|\bm{\mathcal{C}}({\bm{B}})\|_{\infty}\,.$
(58)

From (17), (48), (55) and (58),

$\displaystyle\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Gamma})\|_{\infty}$
$\displaystyle\leq{\epsilon}_{x}{\epsilon}_{y}+M{\epsilon}_{x}+M{\epsilon}_{y}<\epsilon^{2}+2M\epsilon$
(59)

and since $|S|=s_{n}$,

$\displaystyle\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{1,\infty}$
$\displaystyle\leq
s_{n}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{\infty}\leq
s_{n}\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Gamma})\|_{\infty}$
$\displaystyle<s_{n}(\epsilon^{2}+2M\epsilon)\,.$ (60)

By assumption (49),

$\displaystyle\kappa_{\Gamma}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{1,\infty}$
$\displaystyle=\|\bm{\mathcal{C}}((\Gamma_{S,S}^{\ast})^{-1})\|_{1,\infty}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{1,\infty}$
$\displaystyle<\frac{1}{3}\,.$ (61)

By (61) we can invoke [31, Lemma 5] to have

$\displaystyle R({\bm{\Delta}}_{\Gamma})=$
$\displaystyle{\bm{\Gamma}}_{S,S}^{-\ast}({\bm{\Delta}}_{\Gamma})_{S,S}{\bm{\Gamma}}_{S,S}^{-\ast}({\bm{\Delta}}_{\Gamma})_{S,S}{\bm{J}}{\bm{\Gamma}}_{S,S}^{-\ast}$
(62)

where
${\bm{J}}=\sum_{k=0}^{\infty}(-1)^{k}\big{(}{\bm{\Gamma}}_{S,S}^{-\ast}({\bm{\Delta}}_{\Gamma})_{S,S}\big{)}^{k}$.
Using (56), (57) and (62), we have

$\displaystyle\|\bm{\mathcal{C}}(R({\bm{\Delta}}_{\Gamma}))\|_{\infty}\;\leq\;\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast}({\bm{\Delta}}_{\Gamma})_{S,S})\|_{\infty}$
$\displaystyle\quad\quad\quad\quad\times\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast}({\bm{\Delta}}_{\Gamma})_{S,S}J{\bm{\Gamma}}_{S,S}^{-\ast})^{\top}\|_{1,\infty}$
$\displaystyle\leq\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}^{3}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{\infty}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{1,\infty}$
$\displaystyle\quad\quad\quad\quad\times\|\bm{\mathcal{C}}({\bm{J}}^{\top})\|_{1,\infty}\,.$
(63)

Now using (61),

$\displaystyle\|\bm{\mathcal{C}}({\bm{J}}^{\top})\|_{1,\infty}\;\leq\;\sum_{k=0}^{\infty}\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}^{k}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{1,\infty}^{k}$
$\displaystyle=\frac{1}{1-\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{1,\infty}}$
$\displaystyle=\frac{1}{1-\kappa_{\Gamma}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{1,\infty}}\overset{(\ref{aeqn429})}{<}\frac{1}{1-(1/3)}=\frac{3}{2}\,.$
(64)

Using (59), (61), (63) and (64), we have

$\displaystyle\|\bm{\mathcal{C}}(R({\bm{\Delta}}_{\Gamma}))\|_{\infty}\;\leq\;\frac{3}{2}\kappa_{\Gamma}^{3}s_{n}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{\infty}^{2}$
$\displaystyle\quad\quad<\frac{3}{2}\kappa_{\Gamma}^{3}s_{n}(\epsilon^{2}+2M\epsilon)^{2}\,.$
(65)

This proves (51), from which (52) immediately follows.

Using (19), (50), (56), (57) and (59) we have

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{\infty}\leq\|\bm{\mathcal{C}}(R({\bm{\Delta}}_{\Gamma}))\|_{\infty}$
$\displaystyle\quad\quad+\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast}({\bm{\Delta}}_{\Gamma})_{S,S}{\bm{\Gamma}}_{S,S}^{-\ast})\|_{\infty}$
$\displaystyle\leq\|\bm{\mathcal{C}}(R({\bm{\Delta}}_{\Gamma}))\|_{\infty}+\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}^{2}\|\bm{\mathcal{C}}(({\bm{\Delta}}_{\Gamma})_{S,S})\|_{\infty}$
$\displaystyle<\kappa_{\Gamma}^{2}(\epsilon^{2}+2M\epsilon)\big{(}1+1.5s_{n}\kappa_{\Gamma}(\epsilon^{2}+2M\epsilon)\big{)}\,.$
(66)

This proves (53). The claim (54) follows by noting that $|S|=s_{n}$. This
completes the proof. $\quad\Box$

Lemma 4. Assume (49) and the following conditions:

$\displaystyle 0<\alpha<1\mbox{ where }\alpha\mbox{ is as in
}(\ref{eqn335})\,,$ (67)
$\displaystyle\epsilon<\min\left\\{M,\frac{\alpha\lambda_{n}}{2(2-\alpha)}\right\\}\,,$
(68) $\displaystyle\alpha C_{\alpha}\min\\{\lambda_{n},1\\}\geq 3s_{n}\epsilon
M\kappa_{\Gamma}B_{s}$ (69)

where

$\displaystyle C_{\alpha}$
$\displaystyle=\frac{\alpha\lambda_{n}+2\epsilon\alpha-4\epsilon}{2M\alpha\lambda_{n}+\alpha\lambda_{n}+2\epsilon\alpha}\,,$
(70) $\displaystyle B_{s}$
$\displaystyle=\Big{[}1+\kappa_{\Gamma}\Big{(}3s_{n}\epsilon
M+\min\\{s_{n}M^{2},M_{\Sigma}^{2}\\}\Big{)}$
$\displaystyle\quad\quad\times\big{(}4.5s_{n}\epsilon
M\kappa_{\Gamma}+1\big{)}\Big{]}\,.$ (71)

Then we have

  * (i)

$\mbox{bvec}(\hat{\bm{\Delta}}_{S^{c}})={\bm{0}}$.

  * (ii)

$\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{\infty}\leq
2\lambda_{n}\kappa_{\Gamma}+3s_{n}\epsilon
M\kappa_{\Gamma}^{2}\big{(}4.5s_{n}\epsilon
M\kappa_{\Gamma}+1\big{)}\big{(}2M+2\lambda_{n}\big{)}$ $\quad\bullet$

Proof. To establish part (i), we need to show that (36) is true. Let $d$
denote the left-side of (36). It follow from (35) that

$\displaystyle\mbox{bvec}(\tilde{\bm{\Delta}}_{S})=$
$\displaystyle\hat{\bm{\Gamma}}_{S,S}^{-1}\Big{(}\mbox{bvec}((\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})_{S})-\lambda\,\mbox{bvec}({\bm{Z}}(\tilde{\bm{\Delta}}_{S}))\Big{)}\,.$
(72)

Substitute (72) in the left-side of (36) to yield

$\displaystyle d=$
$\displaystyle\|\hat{\bm{\Gamma}}_{e,S}\big{[}\hat{\bm{\Gamma}}_{S,S}^{-1}\big{(}\mbox{bvec}((\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})_{S})-\lambda\,\mbox{bvec}({\bm{Z}}(\tilde{\bm{\Delta}}_{S}))\big{)}\big{]}$
$\displaystyle\quad-\mbox{bvec}((\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})_{e})\|_{2}\,.$
(73)

At the true values we have

$\displaystyle{\bm{0}}=$ $\displaystyle\frac{\partial
L_{\lambda}({\bm{\Delta}},{\bm{\Sigma}}_{x}^{\ast},{\bm{\Sigma}}_{y}^{\ast})}{\partial{\bm{\Delta}}}\Big{|}_{{\bm{\Delta}}={\bm{\Delta}}^{\ast}}={\bm{\Sigma}}_{x}^{\ast}{\bm{\Delta}}^{\ast}{\bm{\Sigma}}_{y}^{\ast}-({\bm{\Sigma}}_{x}^{\ast}-{\bm{\Sigma}}_{y}^{\ast})$

implying

$\displaystyle{\bm{\Gamma}}^{\ast}\mbox{bvec}({\bm{\Delta}}^{\ast})-\mbox{bvec}({\bm{\Sigma}}_{x}^{\ast}-{\bm{\Sigma}}_{y}^{\ast})={\bm{0}}\,,$
(74)

which, noting that $({\bm{\Delta}}^{\ast})_{S^{c}}={\bm{0}}$, can be rewritten
as (cf. (33))

$\displaystyle{\bm{\Gamma}}_{S,S}^{\ast}\mbox{bvec}({\bm{\Delta}}_{S}^{\ast})=$
$\displaystyle\mbox{bvec}({\bm{\Sigma}}_{x}^{\ast})_{S}-\mbox{bvec}({\bm{\Sigma}}_{y}^{\ast})_{S}\,,$
(75)
$\displaystyle{\bm{\Gamma}}_{e,S}^{\ast}\mbox{bvec}({\bm{\Delta}}_{S}^{\ast})=$
$\displaystyle\mbox{bvec}({\bm{\Sigma}}_{x}^{\ast})_{e}-\mbox{bvec}({\bm{\Sigma}}_{y}^{\ast})_{e}\,.$
(76)

Therefore, (${\bm{A}}^{-\ast}=({\bm{A}}^{\ast})^{-1}$),

$\displaystyle{\bm{\Gamma}}_{e,S}^{\ast}$
$\displaystyle{\bm{\Gamma}}_{S,S}^{-\ast}\big{(}\mbox{bvec}({\bm{\Sigma}}_{x}^{\ast})_{S}-\mbox{bvec}({\bm{\Sigma}}_{y}^{\ast})_{S}\big{)}$
$\displaystyle\quad\quad-\mbox{bvec}({\bm{\Sigma}}_{x}^{\ast})_{e}+\mbox{bvec}({\bm{\Sigma}}_{y}^{\ast})_{e}={\bm{0}}\,.$
(77)

Recalling (46) and using (77) in (73),

$\displaystyle d=$
$\displaystyle\|\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}\mbox{bvec}(({\bm{\Delta}}_{\Sigma})_{S})$
$\displaystyle+\big{(}\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast}\big{)}\big{(}\mbox{bvec}({\bm{\Sigma}}_{x}^{\ast})_{S}-\mbox{bvec}({\bm{\Sigma}}_{y}^{\ast})_{S}\big{)}$
$\displaystyle-\lambda\,\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}\mbox{bvec}({\bm{Z}}(\tilde{\bm{\Delta}}_{S}))-\mbox{bvec}(({\bm{\Delta}}_{\Sigma})_{e})\|_{2}\,.$
(78)

We now bound various terms in (78). Note that
$\hat{\bm{\Gamma}}_{e,S}\in\mathbb{R}^{m^{2}\times(m^{2}s_{n})}$,
$\hat{\bm{\Gamma}}_{S,S}^{-1}\in\mathbb{R}^{(m^{2}s_{n})\times(m^{2}s_{n})}$,
and $\mbox{bvec}(({\bm{\Delta}}_{\Sigma})_{S})\in\mathbb{R}^{m^{2}s_{n}}$
where ${\bm{\Delta}}_{\Sigma}\in\mathbb{R}^{(mp_{n})\times(mp_{n})}$. Consider
${\bm{A}}_{e,S}\in\mathbb{R}^{m^{2}\times(m^{2}s_{n})}$. Then

$\displaystyle\|{\bm{A}}_{e,S}\mbox{bvec}(({\bm{\Delta}}_{\Sigma})_{S})\|_{2}$
$\displaystyle=\|\sum_{f\in
S}{\bm{A}}_{e,f}\mbox{vec}(({\bm{\Delta}}_{\Sigma})_{f})\|_{2}$ (79)

where edge $f\in S$, ${\bm{A}}_{e,f}\in\mathbb{R}^{m^{2}\times m^{2}}$ and
$({\bm{\Delta}}_{\Sigma})_{f}\in\mathbb{R}^{m\times m}$. By the triangle
inequality

$\displaystyle\|{\bm{A}}_{e,S}\mbox{bvec}(({\bm{\Delta}}_{\Sigma})_{S})\|_{2}$
$\displaystyle\leq\sum_{f\in
S}\|{\bm{A}}_{e,f}\mbox{vec}(({\bm{\Delta}}_{\Sigma})_{f})\|_{2}\,.$ (80)

With ${\bm{B}}_{i.}$ denoting the $i$th row of matrix ${\bm{B}}$ and using
Cauchy-Schwartz inequality, we have

$\displaystyle\|{\bm{A}}_{e,f}\mbox{vec}(({\bm{\Delta}}_{\Sigma})_{f})\|_{2}=\Big{(}\sum_{i=1}^{m^{2}}\big{(}[{\bm{A}}_{e,f}]_{i.}\mbox{vec}(({\bm{\Delta}}_{\Sigma})_{f})\big{)}^{2}\Big{)}^{1/2}$
$\displaystyle\quad\leq\Big{(}\sum_{i=1}^{m^{2}}\|[{\bm{A}}_{e,f}]_{i.}\|_{2}^{2}\,\|\mbox{vec}(({\bm{\Delta}}_{\Sigma})_{f})\|_{2}^{2}\Big{)}^{1/2}$
$\displaystyle\quad=\|\mbox{vec}(({\bm{\Delta}}_{\Sigma})_{f})\|_{2}\Big{(}\sum_{i=1}^{m^{2}}\|[{\bm{A}}_{e,f}]_{i.}\|_{2}^{2}\Big{)}^{1/2}$
$\displaystyle\quad=\|\mbox{vec}(({\bm{\Delta}}_{\Sigma})_{f})\|_{2}\,\|{\bm{A}}_{e,f}\|_{F}\,.$
(81)

Therefore, using (80) and (81),

$\displaystyle\|{\bm{A}}_{e,S}\mbox{bvec}(({\bm{\Delta}}_{\Sigma})_{S})\|_{2}\leq\Big{(}\sum_{f\in
S}\|{\bm{A}}_{e,f}\|_{F}\Big{)}\max_{g\in
S}\|\mbox{vec}(({\bm{\Delta}}_{\Sigma})_{g})\|_{2}$
$\displaystyle\quad\quad\quad=\|\bm{\mathcal{C}}({\bm{A}}_{e,S})\|_{1}\,\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Sigma})\|_{\infty}\,.$
(82)

Using (82), we have

$\displaystyle\|\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}\mbox{bvec}(({\bm{\Delta}}_{\Sigma})_{S})\|_{2}$
$\displaystyle\quad\leq\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1})\|_{1}\,\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Sigma})\|_{\infty}\,,$
(83)
$\displaystyle\|\big{(}\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast}\big{)}\big{(}\mbox{bvec}({\bm{\Sigma}}_{x}^{\ast})_{S}-\mbox{bvec}({\bm{\Sigma}}_{y}^{\ast})_{S}\big{)}\|_{2}$
$\displaystyle\quad\leq\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1}\,\|\bm{\mathcal{C}}({\bm{\Sigma}}_{x}^{\ast}-{\bm{\Sigma}}_{y}^{\ast})\|_{\infty}\,,$
(84)
$\displaystyle\|\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}\mbox{bvec}({\bm{Z}}(\tilde{\bm{\Delta}}_{S}))\|_{2}$
$\displaystyle\quad\leq\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1})\|_{1}\,\|\bm{\mathcal{C}}({\bm{Z}}(\tilde{\bm{\Delta}}_{S}))\|_{\infty}\,,$
(85)
$\displaystyle\|\mbox{bvec}(({\bm{\Delta}}_{\Sigma})_{e})\|_{2}\leq\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Sigma})\|_{\infty}\,.$
(86)

By (31), (46) and (48)

$\displaystyle\|\bm{\mathcal{C}}({\bm{\Sigma}}_{x}^{\ast}-{\bm{\Sigma}}_{y}^{\ast})\|_{\infty}$
$\displaystyle\leq\|\bm{\mathcal{C}}({\bm{\Sigma}}_{x}^{\ast})\|_{\infty}+\|\bm{\mathcal{C}}({\bm{\Sigma}}_{y}^{\ast})\|_{\infty}\leq
2M\,,$ (87)
$\displaystyle\|\bm{\mathcal{C}}({\bm{Z}}(\tilde{\bm{\Delta}}_{S}))\|_{\infty}$
$\displaystyle\leq 1\,,$ (88)
$\displaystyle\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Sigma})\|_{\infty}$
$\displaystyle\leq\|\bm{\mathcal{C}}({\bm{\Delta}}_{x})\|_{\infty}+\|\bm{\mathcal{C}}({\bm{\Delta}}_{y})\|_{\infty}<2\epsilon\,.$
(89)

Using (78) and (83)-89),

$\displaystyle d\;<\;$ $\displaystyle
2\epsilon\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1})\|_{1}+2M\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1}$
$\displaystyle+\lambda\,\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1})\|_{1}+2\epsilon\,.$
(90)

Therefore, $d<\lambda$ for any edge $e\in S^{c}$ if

$\displaystyle U_{b1}:=\max_{e\in
S^{c}}2M\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1}$
$\displaystyle\quad\quad\quad+2\epsilon(1+\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1})\|_{1})\leq\alpha\lambda_{n}(1-C_{\alpha})\,,$
(91) $\displaystyle U_{b2}:=\max_{e\in
S^{c}}\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1})\|_{1}\leq
1-(1-C_{\alpha})\alpha\,.$ (92)

It remains to show that (91) and (92) are true under the assumptions of Lemma
4. Since

$\displaystyle\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast}=(\hat{\bm{\Gamma}}_{e,S}-{\bm{\Gamma}}_{e,S}^{\ast}){\bm{\Gamma}}_{S,S}^{-\ast}$
$\displaystyle+{\bm{\Gamma}}_{e,S}^{\ast}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})+(\hat{\bm{\Gamma}}_{e,S}-{\bm{\Gamma}}_{e,S}^{\ast})(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\,,$
(93)

we have

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast})\|_{\infty}$
$\displaystyle\quad\leq\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}-{\bm{\Gamma}}_{e,S}^{\ast})\|_{\infty}\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}$
$\displaystyle\quad\quad+\|\bm{\mathcal{C}}({\bm{\Gamma}}_{e,S}^{\ast})\|_{\infty}\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}$
$\displaystyle\quad\quad+\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}-{\bm{\Gamma}}_{e,S}^{\ast})\|_{\infty}\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}\,.$
(94)

With edge $e=\\{i,k\\}\in S^{c}$ and edge $f=\\{j,\ell\\}\in S$, consider

$\displaystyle\hat{\bm{\Gamma}}_{e,f}-{\bm{\Gamma}}_{e,f}^{\ast}=\hat{\bm{\Gamma}}_{ik,j\ell}-{\bm{\Gamma}}_{ik,j\ell}^{\ast}$
$\displaystyle=\hat{\bm{\Sigma}}_{y}^{(ij)}\otimes\hat{\bm{\Sigma}}_{x}^{(k\ell)}-{\bm{\Sigma}}_{y}^{\ast(ij)}\otimes{\bm{\Sigma}}_{x}^{\ast(k\ell)}$
$\displaystyle={\bm{\Delta}}_{y}^{(ij)}\otimes{\bm{\Delta}}_{x}^{(k\ell)}+{\bm{\Sigma}}_{y}^{\ast(ij)}\otimes{\bm{\Delta}}_{x}^{(k\ell)}+{\bm{\Delta}}_{y}^{(ij)}\otimes{\bm{\Sigma}}_{x}^{\ast(k\ell)}.$
(95)

It then follows that

$\displaystyle|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,f}-{\bm{\Gamma}}_{e,f}^{\ast})|\leq\|{\bm{\Delta}}_{y}^{(ij)}\|_{F}\,\|{\bm{\Delta}}_{x}^{(k\ell)}\|_{F}$
$\displaystyle\quad+\|{\bm{\Sigma}}_{y}^{\ast(ij)}\|_{F}\,\|{\bm{\Delta}}_{x}^{(k\ell)}\|_{F}+\|{\bm{\Delta}}_{y}^{(ij)}\|_{F}\,\|{\bm{\Sigma}}_{x}^{\ast(k\ell)}\|_{F}$
$\displaystyle\leq\epsilon_{y}\epsilon_{x}+M\epsilon_{x}+M\epsilon_{y}<\epsilon^{2}+2M\epsilon.$
(96)

Hence

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}-{\bm{\Gamma}}_{e,S}^{\ast})\|_{\infty}<\;$
$\displaystyle\epsilon^{2}+2M\epsilon\,,$ (97)
$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}-{\bm{\Gamma}}_{e,S}^{\ast})\|_{1}<\;$
$\displaystyle s_{n}(\epsilon^{2}+2M\epsilon).$ (98)

Since
${\bm{\Gamma}}_{e,f}^{\ast}={\bm{\Sigma}}_{y}^{\ast(ij)}\otimes{\bm{\Sigma}}_{x}^{\ast(k\ell)}$,
we have $|{\bm{\Gamma}}_{e,f}^{\ast}|\leq M^{2}$ and

$\displaystyle\|\bm{\mathcal{C}}({\bm{\Gamma}}_{e,S}^{\ast})\|_{\infty}$
$\displaystyle\leq
M^{2},\quad\|\bm{\mathcal{C}}({\bm{\Gamma}}_{e,S}^{\ast})\|_{1}\leq
s_{n}M^{2}\,.$ (99)

Alternatively, with $e=\\{i,k\\}\in S^{c}$ and $f=\\{j,\ell\\}\in S$,

$\displaystyle\|\bm{\mathcal{C}}({\bm{\Gamma}}_{e,S}^{\ast})\|_{1}=\sum_{f\in
S}|\bm{\mathcal{C}}({\bm{\Sigma}}_{y}^{\ast(ij)}\otimes{\bm{\Sigma}}_{x}^{\ast(k\ell)})|$
$\displaystyle\quad\leq\sum_{\\{j,\ell\\}\in
S}\|{\bm{\Sigma}}_{y}^{\ast(ij)}\|_{F}\,\|{\bm{\Sigma}}_{x}^{\ast(k\ell)}\|_{F}$
$\displaystyle\quad\leq(\sum_{j=1}^{p}\|{\bm{\Sigma}}_{y}^{\ast(ij)}\|_{F})(\sum_{\ell=1}^{p}\|{\bm{\Sigma}}_{x}^{\ast(k\ell)}\|_{F})$
$\displaystyle\quad\leq\|\bm{\mathcal{C}}({\bm{\Sigma}}_{y}^{\ast})\|_{1,\infty}\,\|\bm{\mathcal{C}}({\bm{\Sigma}}_{x}^{\ast})\|_{1,\infty}\leq
M_{\Sigma}^{2}\,.$ (100)

From (93)-(100) and Lemma 3, we have

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1}$
$\displaystyle\quad\leq\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}-{\bm{\Gamma}}_{e,S}^{\ast})\|_{1}\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}$
$\displaystyle\quad\quad+\|\bm{\mathcal{C}}({\bm{\Gamma}}_{e,S}^{\ast})\|_{1}\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}$
$\displaystyle\quad\quad+\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}-{\bm{\Gamma}}_{e,S}^{\ast})\|_{1}\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}$
$\displaystyle\quad\leq
s_{n}(\epsilon^{2}+2M\epsilon)\kappa_{\Gamma}+\big{[}\min\\{s_{n}M^{2},M_{\Sigma}^{2}\\}+s_{n}(\epsilon^{2}+2M\epsilon)\big{]}$
$\displaystyle\quad\quad\quad\times\big{[}s_{n}(\epsilon^{2}+2M\epsilon)\kappa_{\Gamma}^{2}\big{]}\big{(}1+1.5s_{n}(\epsilon^{2}+2M\epsilon)\kappa_{\Gamma}\big{)}$
$\displaystyle\quad\overset{\epsilon<M}{\leq}3s_{n}\epsilon
M\kappa_{\Gamma}B_{s}\,\leq\,\alpha C_{\alpha}\min\\{\lambda_{n},1\\}\,,$
(101)

where $B_{s}$ is as in (71) and we used $\epsilon<M$ to infer
$\epsilon^{2}+2M\epsilon<3M\epsilon$. Using the triangle inequality
$|a|-|b|\leq|a-b|\leq|a|+|b|$, we have
$\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1}$
$\geq$
$\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1})\|_{1}$
$-\|\bm{\mathcal{C}}({\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1}$
, which, using (67) and (101), leads to

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{e,S}\hat{\bm{\Gamma}}_{S,S}^{-1})\|_{1}\leq\|\bm{\mathcal{C}}({\bm{\Gamma}}_{e,S}^{\ast}{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1}+\alpha
C_{\alpha}\min\\{\lambda_{n},1\\}$ $\displaystyle\quad\leq 1-\alpha+\alpha
C_{\alpha}\min\\{\lambda_{n},1\\}\leq 1-(1-C_{\alpha})\alpha\,.$ (102)

This establishes (92). To show (91), using (101)-(102),

$\displaystyle U_{b1}\;\leq\;$ $\displaystyle 2M\alpha
C_{\alpha}\min\\{\lambda_{n},1\\}+2\epsilon(1+1-(1-C_{\alpha})\alpha)$
$\displaystyle\leq\;$ $\displaystyle 2M\alpha
C_{\alpha}\lambda_{n}+2\epsilon(2-(1-C_{\alpha})\alpha)\overset{(\ref{aeqn504})}{=}\alpha\lambda_{n}(1-C_{\alpha})\,.$
(103)

This proves (91), and thus, part (i) of Lemma 4.

We now turn to the proof of Lemma 4(ii). Since
$\hat{\bm{\Delta}}=\tilde{\bm{\Delta}}$, for any edge $\\{k,\ell\\}\in S$, we
have

$\displaystyle\|(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})^{(k\ell)}\|_{F}=\|(\tilde{\bm{\Delta}}-{\bm{\Delta}}^{\ast})^{(k\ell)}\|_{F}$
$\displaystyle\quad=\|\mbox{vec}(\tilde{\bm{\Delta}}^{(k\ell)})-\mbox{vec}(({\bm{\Delta}}^{\ast})^{(k\ell)})\|_{2}\,.$
(104)

Using (35) and (75)

$\displaystyle\mbox{bvec}((\tilde{\bm{\Delta}}-{\bm{\Delta}}^{\ast})_{S})=\hat{\bm{\Gamma}}_{S,S}^{-1}\mbox{bvec}(({\bm{\Delta}}_{\Sigma})_{S})+(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})$
$\displaystyle\quad\times\mbox{bvec}(({\bm{\Sigma}}_{x}^{\ast}-{\bm{\Sigma}}_{y}^{\ast})_{S})-\lambda_{n}\hat{\bm{\Gamma}}_{S,S}^{-1}\,\mbox{bvec}({\bm{Z}}(\tilde{\bm{\Delta}}_{S}))\,.$
(105)

Since
$\hat{\bm{\Gamma}}_{S,S}^{-1}=\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast}+{\bm{\Gamma}}_{S,S}^{-\ast}\,,$

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1})\|_{1,\infty}\leq\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}+\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}\,.$
(106)

By (105), for any edge $f=\\{k,\ell\\}\in S$, we have

$\displaystyle\|\mbox{vec}((\tilde{\bm{\Delta}}-{\bm{\Delta}}^{\ast})^{(k\ell)})\|_{2}\leq\|(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})_{f,S}$
$\displaystyle\quad\quad\quad\times\mbox{bvec}\big{(}({\bm{\Delta}}_{\Sigma})_{S}+({\bm{\Sigma}}_{x}^{\ast}-{\bm{\Sigma}}_{y}^{\ast})_{S}-\lambda_{n}{\bm{Z}}(\tilde{\bm{\Delta}}_{S})\big{)}\|_{2}$
$\displaystyle\quad+\|({\bm{\Gamma}}_{S,S}^{-\ast})_{f,S}\,\mbox{bvec}\Big{(}({\bm{\Delta}}_{\Sigma})_{S}-\lambda_{n}{\bm{Z}}(\tilde{\bm{\Delta}}_{S})\Big{)}\|_{2}$
$\displaystyle\leq\|\bm{\mathcal{C}}(\hat{\bm{\Gamma}}_{S,S}^{-1}-{\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}\,\Big{(}\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Sigma})\|_{\infty}+\|\bm{\mathcal{C}}({\bm{\Sigma}}_{x}^{\ast}-{\bm{\Sigma}}_{y}^{\ast})\|_{\infty}$
$\displaystyle\quad\quad+\lambda_{n}\Big{)}+\|\bm{\mathcal{C}}({\bm{\Gamma}}_{S,S}^{-\ast})\|_{1,\infty}\,\big{(}\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Sigma})\|_{\infty}+\lambda_{n}\big{)}$
$\displaystyle\leq
s_{n}\kappa_{\Gamma}^{2}(\epsilon^{2}+2M\epsilon)\big{(}1+1.5s_{n}(\epsilon^{2}+2M\epsilon)\kappa_{\Gamma}\big{)}$
$\displaystyle\quad\quad\quad\times(2\epsilon+2M+\lambda_{n})+\kappa_{\Gamma}(2\epsilon+\lambda_{n})=:U_{b3}\,.$
(107)

By (68), for $0<\alpha<1$, we have
$2\epsilon<\alpha\lambda_{n}/(2-\alpha)<\alpha\lambda_{n}<\lambda_{n}$.
Therefore,
$\kappa_{\Gamma}(2\epsilon+\lambda_{n})<2\kappa_{\Gamma}\lambda_{n}$ and
$2\epsilon+2M+\lambda_{n}<2M+2\lambda_{n}$. Since $\epsilon<M$ by (68), we
also have $\epsilon^{2}+2M\epsilon<3M\epsilon$. Using these relations and
(107), it follows that

$\displaystyle U_{b3}\;\leq\;$ $\displaystyle 3s_{n}\epsilon
M\kappa_{\Gamma}^{2}\big{(}1+4.5s_{n}\epsilon
M\kappa_{\Gamma}\big{)}\big{(}2M+2\lambda_{n}\big{)}+2\lambda_{n}\kappa_{\Gamma}\,.$

Finally,

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{\infty}=$
$\displaystyle\max_{f=\\{k,\ell\\}\in
S}\|\mbox{vec}((\tilde{\bm{\Delta}}-{\bm{\Delta}}^{\ast})^{(k\ell)})\|_{2}\,,$

proving the desired result. $\quad\Box$

Proof of Theorem 1. Here we first show that under the sufficient conditions of
Theorem 1, the assumptions of Lemmas 2-4 holds true. We pick
$\epsilon=C_{0}\sqrt{ln(p_{n})/n}$, implying, by Lemma 2, that
$\|\bm{\mathcal{C}}(\hat{\bm{\Sigma}}_{x}-{\bm{\Sigma}}_{x}^{\ast})\|_{\infty}\leq\epsilon$
and
$\|\bm{\mathcal{C}}(\hat{\bm{\Sigma}}_{y}-{\bm{\Sigma}}_{y}^{\ast})\|_{\infty}\leq\epsilon$.
with probability $\geq 1-2/p_{n}^{\tau-2}$, $\tau>2$. We first show that with
this choice of $\epsilon$, condition (68) of Lemma 4 holds true. By the choice
of $\lambda_{n}$, we have $\lambda_{n}\geq 8\epsilon/\alpha$. Clearly, for
$0<\alpha<1$, $(\alpha/8)<\alpha/(4(2-\alpha))$. Therefore

$\displaystyle\epsilon$
$\displaystyle\leq\frac{\alpha\lambda_{n}}{8}<\frac{\alpha\lambda_{n}}{4(2-\alpha)}<\frac{\alpha\lambda_{n}}{2(2-\alpha)}\,.$
(108)

By a choice of $n$ in (24), we have $n>C_{0}^{2}\ln(p_{n})/M^{2}$. Hence,
$\epsilon=C_{0}\sqrt{ln(p_{n})/n}<M$. Thus, (68) of Lemma 4 holds true. Next
we show that condition (49) of Lemma 3 holds. By a choice of $n$ in (24), we
have $n>81C_{0}^{2}\ln(p_{n})M^{2}s_{n}^{2}\kappa_{\Gamma}^{2}$. Therefore,

$\displaystyle\kappa_{\Gamma}$
$\displaystyle<\frac{1}{C_{0}}\sqrt{\frac{n}{\ln(p_{n})}}\times\frac{1}{9s_{n}M}=\frac{1}{9s_{n}M\epsilon}<\frac{1}{3s_{n}(\epsilon^{2}+2M\epsilon)}$
(109)

since $\epsilon<M$. This proves (49) of Lemma 3 holds.

Now we show that (69) of Lemma 4 holds. Since
$\epsilon<\alpha\lambda_{n}/(4(2-\alpha))$, by (70), we have

$\displaystyle C_{\alpha}$
$\displaystyle=\frac{\alpha\lambda_{n}+2\epsilon\alpha-4\epsilon}{2M\alpha\lambda_{n}+\alpha\lambda_{n}+2\epsilon\alpha}>\frac{\alpha\lambda_{n}-4\epsilon}{2M\alpha\lambda_{n}+\alpha\lambda_{n}+2\epsilon\alpha}$
$\displaystyle>\frac{\alpha\lambda_{n}-\alpha\lambda_{n}/(2-\alpha)}{\alpha\lambda_{n}(2M+1)+2\epsilon\alpha}=\frac{1-\alpha}{(2-\alpha)\big{[}2M+1+\frac{2\epsilon\alpha}{\lambda_{n}\alpha}\big{]}}$
$\displaystyle>\frac{1-\alpha}{(2-\alpha)(2M+1)+\alpha}=\bar{C}_{\alpha}$
(110)

where in the last inequality above, we used
$\epsilon<\alpha\lambda_{n}/(2(2-\alpha))$ from (108). Consider the right-side
$3s_{n}\epsilon M\kappa_{\Gamma}B_{s}$ of (69). From (109),
$\kappa_{\Gamma}<1/(9s_{n}M\epsilon)$. Therefore,

$\displaystyle B_{s}<$ $\displaystyle
1+\kappa_{\Gamma}\big{(}\min\\{s_{n}M^{2},M_{\Sigma}^{2}\\}+\frac{1}{3\kappa_{\Gamma}}\big{)}\big{(}\frac{1}{2}+1\big{)}$
$\displaystyle=$ $\displaystyle
1.5+1.5\kappa_{\Gamma}\min\\{s_{n}M^{2},M_{\Sigma}^{2}\\}=C_{M\kappa}\,.$
(111)

By (23), $\lambda_{n}\geq
4.5\epsilon(\alpha\bar{C}_{\alpha})^{-1}s_{n}M\kappa_{\Gamma}(1+\kappa_{\Gamma}\min\\{s_{n}M^{2},M_{\Sigma}^{2}\\})$.
Hence, using (110), we have

$\displaystyle 3s_{n}\epsilon M\kappa_{\Gamma}B_{s}<$
$\displaystyle\alpha\bar{C}_{\alpha}\lambda_{n}<\alpha{C}_{\alpha}\lambda_{n}\,,$
(112)

proving part of (69) for our choice of $\lambda_{n}$. To show that we also
have $3s_{n}\epsilon M\kappa_{\Gamma}B_{s}<\alpha\bar{C}_{\alpha}$, consider
the choice of $n$ in (24) given by

$\displaystyle n>\;$
$\displaystyle{C}_{0}^{2}\ln(p_{n})\,\frac{9s_{n}^{2}}{(\alpha\bar{C}_{\alpha})^{2}}(\kappa_{\Gamma}MC_{M\kappa})^{2}\,.$
(113)

Then

$\displaystyle\epsilon=\;$
$\displaystyle{C}_{0}\sqrt{\frac{\ln(p_{n})}{n}}\,<\;\frac{\alpha\bar{C}_{\alpha}}{3s_{n}M\kappa_{\Gamma}C_{M\kappa}}\,,$
(114)

and from (112),

$\displaystyle 3s_{n}\epsilon M\kappa_{\Gamma}B_{s}<\;$
$\displaystyle\alpha\bar{C}_{\alpha}<\alpha{C}_{\alpha}\,.$ (115)

Thus, all assumptions of Lemma 4 hold true.

Therefore, Lemma 4(i) applies, proving Theorem 1(ii). By (109),
$9s_{n}\epsilon M\kappa_{\Gamma}<1$. Using this fact in Lemma 4(ii),

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{\infty}\leq
2\lambda_{n}\kappa_{\Gamma}+9s_{n}\epsilon
M\kappa_{\Gamma}^{2}\big{(}M+\lambda_{n}\big{)}$ $\displaystyle\quad\quad\leq
3\lambda_{n}\kappa_{\Gamma}+9s_{n}\epsilon M^{2}\kappa_{\Gamma}^{2}\,.$ (116)

Since, by (23),
$3\lambda_{n}\kappa_{\Gamma}={C}_{0}\sqrt{\ln(p_{n})/n}\,C_{b1}$ and we picked
$\epsilon={C}_{0}\sqrt{\ln(p_{n})/n}$, we have

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{\infty}\leq(C_{b1}+C_{b2}){C}_{0}\sqrt{\ln(p_{n})/n}\,,$

proving Theorem 1(i). To prove part (iii), since
$\hat{\bm{\Delta}}_{S^{c}}=\tilde{\bm{\Delta}}_{S^{c}}={\bm{\Delta}}_{S^{c}}^{\ast}={\bm{0}}$,
we have

$\displaystyle\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{F}=\Big{(}\sum_{\\{k,\ell\\}\in
S}\|\hat{\bm{\Delta}}^{(k\ell)}-({\bm{\Delta}}^{\ast})^{(k\ell)}\|_{F}^{2}\Big{)}^{1/2}$
$\displaystyle\quad\quad\quad\leq\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{\infty}\,\sqrt{s_{n}}\,.$
(117)

Finally, to establish part (iv), note that parts (i)-(iii) hold with
probability $>1-2/p_{n}^{\tau-2}$ (with high probability (w.h.p.)). Recall
that ${\cal G}_{\Delta}=\left(V,{\cal E}_{\Delta}\right)$ denotes the MA
differential graph with edgeset ${\cal
E}_{\Delta}=\\{\\{k,\ell\\}\,:\,\|{\bm{\Delta}}^{(k\ell)}\|_{F}>0\\}$. Let
${\cal G}_{\Delta^{\ast}}$ and ${\cal G}_{\hat{\Delta}}$ denoted true and
estimated graphs based on ${\bm{\Delta}}^{\ast}$ and $\hat{\bm{\Delta}}$,
respectively. If $\min_{(k,\ell)\in
S}\|({\bm{\Delta}}^{\ast})^{(k\ell)}\|_{F}\geq
2\|\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{\infty}$, then
$\bm{\mathcal{C}}(\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})\|_{\infty}=\bm{\mathcal{C}}((\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast})_{S})\|_{\infty}\leq(1/2)\min_{(k,\ell)\in
S}\|({\bm{\Delta}}^{\ast})^{(k\ell)}\|_{F}$, therefore, $\min_{(k,\ell)\in
S}\|(\hat{\bm{\Delta}}_{S})^{(k\ell)}\|_{F}\geq(1/2)\min_{(k,\ell)\in
S}\|({\bm{\Delta}}^{\ast})^{(k\ell)}\|_{F}>0$, while
$\hat{\bm{\Delta}}_{S^{c}}={\bm{0}}$ w.h.p. $\quad\Box$

##  Appendix B Technical Lemmas and Proof of Theorem 2

In order to invoke [30], we first vectorize (3), using
${\bm{\theta}}=\mbox{bvec}({\bm{\Delta}})\in\mathbb{R}^{m^{2}p^{2}}$, as (cf.
(32))

${\cal
L}({\bm{\theta}})=\frac{1}{2}{\bm{\theta}}^{\top}(\hat{\bm{\Sigma}}_{y}\boxtimes\hat{\bm{\Sigma}}_{x}){\bm{\theta}}-{\bm{\theta}}^{\top}\mbox{bvec}(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})$
(118)

where previous $L({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})$
is now ${\cal L}({\bm{\theta}})$. To include sparse-group penalty, recall that
the submatrix ${\bm{\Delta}}^{(k\ell)}$ of ${\bm{\Delta}}$ corresponds to the
edge $\\{k,\ell\\}$ of the MA graph. We denote its vectorized version as
${\bm{\theta}}_{Gt}\in\mathbb{R}^{m^{2}}$ (subscript $G$ for grouped variables
[30]) with index $t=1,2,\cdots,p^{2}$. Then
${\bm{\theta}}_{Gt}=\mbox{vec}({\bm{\Delta}}^{(k\ell)})$ where
$t=(k-1)p+\ell$, $\ell=t\mod{p}$, and $k=\lfloor t/p\rfloor+1$. Using this
notation, the penalty
$\lambda\sum_{k,\ell=1}^{p}\|{\bm{\Delta}}^{(k\ell)}\|_{F}=\lambda\sum_{t=1}^{p^{2}}\|{\bm{\theta}}_{Gt}\|_{2}$.
In the notation of [30], the regularization penalty

${\cal R}({\bm{\theta}})=\|{\bm{\theta}}\|_{\bar{\cal
G},2}:=\sum_{t=1}^{p^{2}}\|{\bm{\theta}}_{Gt}\|_{2}$ (119)

where the index set $\\{1,2,\cdots,(mp)^{2}\\}$ is partitioned into a set of
$N_{G}=p^{2}$ disjoint groups $\bar{\cal
G}=\\{G_{1},G_{2},\cdots,G_{p^{2}}\\}$. As shown in [30, Sec. 2.2], ${\cal
R}({\bm{\theta}})$ is a norm. The counterpart to penalized
$L_{\lambda}({\bm{\Delta}},\hat{\bm{\Sigma}}_{x},\hat{\bm{\Sigma}}_{y})$ of
(4) is (we denote $\lambda$ by $\lambda_{n}$, as in Appendix A)

$\displaystyle{\cal L}_{\lambda}({\bm{\theta}})=$ $\displaystyle{\cal
L}({\bm{\theta}})+\lambda_{n}{\cal R}({\bm{\theta}})\,.$ (120)

As discussed in [30, Sec. 2.2], w.r.t. the usual Euclidean inner product
$\langle{\bm{u}},{\bm{v}}\rangle={\bm{u}}^{\top}{\bm{v}}$ for
${\bm{u}},{\bm{v}}\in\mathbb{R}^{m^{2}p^{2}}$ and given any subset
$S_{\bar{\cal G}}\subseteq\\{1,2,\cdots,N_{G}\\}$ of group indices, define the
subspace

${\mathcal{M}}=\\{{\bm{\theta}}\in\mathbb{R}^{m^{2}p^{2}}\,|\,{\bm{\theta}}_{Gt}={\bm{0}}\mbox{
for all }t\not\in S_{\bar{\cal G}}\\}$ (121)

and its orthogonal complement

${\mathcal{M}}^{\perp}=\\{{\bm{\theta}}\in\mathbb{R}^{m^{2}p^{2}}\,|\,{\bm{\theta}}_{Gt}={\bm{0}}\mbox{
for all }t\in S_{\bar{\cal G}}\\}\,.$ (122)

The chosen ${\cal R}({\bm{\theta}})$ is decomposable w.r.t.
$({\mathcal{M}},{\mathcal{M}}^{\perp})$ since ${\cal
R}({\bm{\theta}}^{(1)}+{\bm{\theta}}^{(2)})={\cal
R}({\bm{\theta}}^{(1)})+{\cal R}({\bm{\theta}}^{(2)})$ for any
${\bm{\theta}}^{(1)}\in{\mathcal{M}}$ and
${\bm{\theta}}^{(2)}\in{\mathcal{M}}^{\perp}$ [30, Sec. 2.2, Example 2].

In order to invoke [30], we need the dual norm ${\cal R}^{\ast}$ of
regularizer ${\cal R}$ w.r.t. the inner product
$\langle{\bm{u}},{\bm{v}}\rangle={\bm{u}}^{\top}{\bm{v}}$. It is given by [30,
Eqns. (14)-(15)]

$\displaystyle{\cal R}^{\ast}({\bm{v}})=$ $\displaystyle\sup_{{\cal
R}({\bm{u}})\leq 1}\langle{\bm{u}},{\bm{v}}\rangle\;=\;\max_{t=1,2,\cdots
p^{2}}\|{\bm{u}}_{Gt}\|_{2}\,.$ (123)

We also need the subspace compatibility index [30], defined as

$\displaystyle\Psi({\mathcal{M}})=$
$\displaystyle\sup_{{\bm{u}}\in{\mathcal{M}}\backslash\\{0\\}}\frac{{\cal
R}({\bm{u}})}{\|{\bm{u}}\|_{2}}\,.$ (124)

For group lasso penalty, $\Psi({\mathcal{M}})=\sqrt{s_{n}}\,$ [30, Sec. 9.2
(Supplementary)], where $s_{n}=|S_{\bar{\cal G}}|=$ number of edges in the
true MA differential graph. We need to establish a restricted strong convexity
condition [30] on ${\cal L}({\bm{\theta}})$. With
${\bm{\theta}}^{\ast}=\mbox{bvec}({\bm{\Delta}}^{\ast})$ denoting the true
value, and ${\bm{\theta}}={\bm{\theta}}^{\ast}+\tilde{\bm{\theta}}$, define

$\displaystyle\delta{\cal L}(\tilde{\bm{\theta}},{\bm{\theta}}^{\ast})=$
$\displaystyle{\cal L}({\bm{\theta}}^{\ast}+\tilde{\bm{\theta}})-{\cal
L}({\bm{\theta}}^{\ast})-\langle\nabla{\cal
L}({\bm{\theta}}^{\ast}),\tilde{\bm{\theta}}\rangle$ (125)

where the gradient $\nabla{\cal L}({\bm{\theta}}^{\ast})$ at
${\bm{\theta}}={\bm{\theta}}^{\ast}$ is

$\displaystyle\nabla{\cal L}({\bm{\theta}}^{\ast})=$
$\displaystyle(\hat{\bm{\Sigma}}_{y}\boxtimes\hat{\bm{\Sigma}}_{x}){\bm{\theta}}^{\ast}-\mbox{bvec}(\hat{\bm{\Sigma}}_{x}-\hat{\bm{\Sigma}}_{y})\,.$
(126)

Hence (125) simplifies to

$\displaystyle\delta{\cal L}(\tilde{\bm{\theta}},{\bm{\theta}}^{\ast})=$
$\displaystyle\tilde{\bm{\theta}}^{\top}(\hat{\bm{\Sigma}}_{y}\boxtimes\hat{\bm{\Sigma}}_{x})\tilde{\bm{\theta}}=\tilde{\bm{\theta}}^{\top}\hat{\bm{\Gamma}}\tilde{\bm{\theta}}\,,$
(127)

which may be rewritten as

$\displaystyle\delta{\cal L}(\tilde{\bm{\theta}},{\bm{\theta}}^{\ast})=$
$\displaystyle\tilde{\bm{\theta}}^{\top}{\bm{\Gamma}}^{\ast}\tilde{\bm{\theta}}+\tilde{\bm{\theta}}^{\top}(\hat{\bm{\Gamma}}-{\bm{\Gamma}}^{\ast})\tilde{\bm{\theta}}\,.$
(128)

By the sparsity assumption,
${\bm{\theta}}^{\ast}={\bm{\theta}}^{\ast}_{\mathcal{M}}$, hence,
${\bm{\theta}}^{\ast}_{{\mathcal{M}}^{\perp}}={\bm{0}}$, where
${\bm{\theta}}_{\mathcal{M}}$ and ${\bm{\theta}}_{{\mathcal{M}}^{\perp}}$
denote projection of ${\bm{\theta}}$ on subspaces ${\mathcal{M}}$ and
${\mathcal{M}}^{\perp}$, respectively.

Similar to (5), suppose

$\hat{\bm{\theta}}=\arg\min_{\bm{\theta}}\big{\\{}{\cal
L}({\bm{\theta}})+\lambda_{n}{\cal R}({\bm{\theta}})\big{\\}}\,,$ (129)

and we consider (125) and (127) with
$\hat{\bm{\theta}}={\bm{\theta}}^{\ast}+\tilde{\bm{\theta}}$. Then

$\displaystyle\hat{\bm{\theta}}$
$\displaystyle-{\bm{\theta}}^{\ast}=\hat{\bm{\theta}}_{\mathcal{M}}-{\bm{\theta}}^{\ast}+\hat{\bm{\theta}}_{{\mathcal{M}}^{\perp}}=\tilde{\bm{\theta}}_{\mathcal{M}}+\tilde{\bm{\theta}}_{{\mathcal{M}}^{\perp}}\,.$
(130)

By [30, Lemma 1],

$\displaystyle{\cal R}(\tilde{\bm{\theta}}_{{\mathcal{M}}^{\perp}})\;\leq\;$
$\displaystyle 3{\cal R}(\tilde{\bm{\theta}}_{{\mathcal{M}}})+4{\cal
R}({\bm{\theta}}^{\ast}_{{\mathcal{M}}^{\perp}})\,,$ (131)

if

$\displaystyle\lambda_{n}\geq$ $\displaystyle 2{\cal R}^{\ast}(\nabla{\cal
L}({\bm{\theta}}^{\ast}))\,.$ (132)

Since in our case ${\bm{\theta}}^{\ast}_{{\mathcal{M}}^{\perp}}={\bm{0}}$, we
have ${\cal R}({\bm{\theta}}^{\ast}_{{\mathcal{M}}^{\perp}})=0$.
Lemma 5. Under (15) and using the notation of Appendix A,

$\displaystyle{\cal R}^{\ast}(\nabla{\cal
L}({\bm{\theta}}^{\ast}))\;\leq\;(\epsilon^{2}+2M\epsilon)s_{n}\,\max_{t=1,\cdots,p^{2}}\|{\bm{\theta}}_{Gt}^{\ast}\|_{2}+2\epsilon\,\quad\bullet$

Proof. Using (46), (74) and (126), we have

$\displaystyle\nabla{\cal L}({\bm{\theta}}^{\ast})=$
$\displaystyle{\bm{\Delta}}_{\Gamma}{\bm{\theta}}^{\ast}+\mbox{bvec}({\bm{\Delta}}_{y})-\mbox{bvec}({\bm{\Delta}}_{x})\,.$
(133)

Expressing it group-wise, with groups $t$ and $t_{1}$ corresponding to edges
$\\{j,k\\}$ and $\\{\ell,q\\}$, respectively,

$\displaystyle(\nabla{\cal L}({\bm{\theta}}^{\ast}))_{Gt_{1}}=$
$\displaystyle\sum_{t=1}^{p^{2}}({\bm{\Delta}}_{\Gamma})_{Gt_{1},Gt}{\bm{\theta}}^{\ast}_{Gt}+\mbox{bvec}({\bm{\Delta}}_{y})_{Gt_{1}}$
$\displaystyle\quad-\mbox{bvec}({\bm{\Delta}}_{x})_{Gt_{1}}\,.$ (134)

Therefore, by the Cauchy-Schwartz inequality, and using (47), (48) and (59),
we have

$\displaystyle\|(\nabla{\cal
L}({\bm{\theta}}^{\ast}))_{Gt_{1}}\|_{2}\leq\sum_{t=1}^{p^{2}}\|({\bm{\Delta}}_{\Gamma})_{Gt_{1},Gt}\|_{F}\,\|{\bm{\theta}}^{\ast}_{Gt}\|_{2}$
$\displaystyle\quad\quad\quad+\|\mbox{bvec}({\bm{\Delta}}_{y})_{Gt_{1}}\|_{2}+\|\mbox{bvec}({\bm{\Delta}}_{x})_{Gt_{1}}\|_{2}$
$\displaystyle\quad\leq\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Gamma})\|_{\infty}\,\sum_{t=1}^{p^{2}}\|{\bm{\theta}}^{\ast}_{Gt}\|_{2}+\|{\bm{\Delta}}_{y}^{(\ell
q)}\|_{\infty}+\|{\bm{\Delta}}_{y}^{(\ell q)}\|_{\infty}$
$\displaystyle\quad\leq(\epsilon^{2}+2M\epsilon)s_{n}\,\max_{t=1,\cdots,p^{2}}\|{\bm{\theta}}_{Gt}^{\ast}\|_{2}+\epsilon+\epsilon\,.$
(135)

By (123) and (135) we have the desired result. $\quad\Box$

Lemma 6. Under (15) and the notation of Appendix A,

$\displaystyle\delta{\cal
L}(\tilde{\bm{\theta}},{\bm{\theta}}^{\ast})\geq\kappa_{\cal
L}\,\|\tilde{\bm{\theta}}\|_{2}^{2}$ (136)

where $\kappa_{\cal
L}=\frac{1}{2}\phi_{min}^{\ast}-8s_{n}(\epsilon^{2}+2M\epsilon)$.
$\quad\bullet$
Proof. We have

$\displaystyle\tilde{\bm{\theta}}^{\top}(\hat{\bm{\Gamma}}-{\bm{\Gamma}}^{\ast})\tilde{\bm{\theta}}=\sum_{t_{1}=1}^{p^{2}}\sum_{t_{2}=1}^{p^{2}}\tilde{\bm{\theta}}^{\top}_{Gt_{1}}({\bm{\Delta}}_{\Gamma})_{Gt_{1},Gt_{2}}\tilde{\bm{\theta}}_{Gt_{2}}\,.$
(137)

Therefore,

$\displaystyle|\tilde{\bm{\theta}}^{\top}(\hat{\bm{\Gamma}}-{\bm{\Gamma}}^{\ast})\tilde{\bm{\theta}}|\leq\sum_{t_{1}=1}^{p^{2}}\sum_{t_{2}=1}^{p^{2}}|\tilde{\bm{\theta}}^{\top}_{Gt_{1}}({\bm{\Delta}}_{\Gamma})_{Gt_{1},Gt_{2}}\tilde{\bm{\theta}}_{Gt_{2}}|$
$\displaystyle\quad\leq\sum_{t_{1}=1}^{p^{2}}\sum_{t_{2}=1}^{p^{2}}\|\tilde{\bm{\theta}}_{Gt_{1}}\|_{2}\,\|({\bm{\Delta}}_{\Gamma})_{Gt_{1},Gt_{2}}\|_{F}\,\|\tilde{\bm{\theta}}_{Gt_{2}}\|_{2}$
$\displaystyle\quad\leq\|\bm{\mathcal{C}}({\bm{\Delta}}_{\Gamma})\|_{\infty}\sum_{t_{1}=1}^{p^{2}}\sum_{t_{2}=1}^{p^{2}}\|\tilde{\bm{\theta}}_{Gt_{1}}\|_{2}\,\|\tilde{\bm{\theta}}_{Gt_{2}}\|_{2}$
$\displaystyle\quad\leq(\epsilon^{2}+2M\epsilon)\,\|\tilde{\bm{\theta}}\|_{\bar{\cal
G},2}^{2}\;,$ (138)

where we used (119). We have

$\displaystyle\|\tilde{\bm{\theta}}\|_{\bar{\cal
G},2}^{2}=\|\tilde{\bm{\theta}}_{\mathcal{M}}+\tilde{\bm{\theta}}_{{\mathcal{M}}^{\perp}}\|_{\bar{\cal
G},2}^{2}=(\|\tilde{\bm{\theta}}_{\mathcal{M}}\|_{\bar{\cal
G},2}+\|\tilde{\bm{\theta}}_{{\mathcal{M}}^{\perp}}\|_{\bar{\cal G},2})^{2}$
$\displaystyle\quad\overset{(\ref{aeqn734})}{\leq}16\,\|\tilde{\bm{\theta}}_{\mathcal{M}}\|_{\bar{\cal
G},2}^{2}\overset{(\ref{aeqn716})}{\leq}16s_{n}\|\tilde{\bm{\theta}}_{\mathcal{M}}\|_{2}^{2}\leq
16s_{n}\|\tilde{\bm{\theta}}\|_{2}^{2}\,.$ (139)

Noting that
$\tilde{\bm{\theta}}^{\top}{\bm{\Gamma}}^{\ast})\tilde{\bm{\theta}}\geq\phi_{min}^{\ast}\,\|\tilde{\bm{\theta}}\|_{2}^{2}$
and using (128), (138) and (139), we have

$\displaystyle\delta{\cal
L}(\tilde{\bm{\theta}},{\bm{\theta}}^{\ast})\geq\big{(}\frac{1}{2}\phi_{min}^{\ast}-8s_{n}(\epsilon^{2}+2M\epsilon)\big{)}\|\tilde{\bm{\theta}}\|_{2}^{2}=\kappa_{\cal
L}\,\|\tilde{\bm{\theta}}\|_{2}^{2}\,,$

proving the desired result. $\quad\Box$

Proof of Theorem 2. First choose $\epsilon$ to make $\kappa_{\cal L}>0$ in
Lemma 6. For instance, suppose we take
$8s_{n}(\epsilon^{2}+2M\epsilon)\leq\phi_{min}^{\ast}/4$. Then $\kappa_{\cal
L}>\phi_{min}^{\ast}/4$. Now pick

$\displaystyle\epsilon=$ $\displaystyle
C_{0}\sqrt{ln(p_{n})/n}\;\leq\;\min\big{\\{}M,\frac{\phi_{min}^{\ast}}{96s_{n}M}\big{\\}}\,,$
(140)

leading to $8s_{n}(\epsilon^{2}+2M\epsilon)\leq
24s_{n}M\epsilon\leq\phi_{min}^{\ast}/4$. These upper bounds can be ensured by
picking appropriate lower bounds to sample size $n$ and invoking Lemma 2. The
choice of $n$ specified in (26) satisfies (140) with probability
$>1-2/p_{n}^{\tau-2}$. Using $\epsilon=C_{0}\sqrt{ln(p_{n})/n}\leq M$, the
lower bound on $\lambda_{n}$ given in (25) satisfies (132) with ${\cal
R}^{\ast}(\nabla{\cal L}({\bm{\theta}}^{\ast}))$ as in Lemma 5. By [30,
Theorem 1], $\hat{\bm{\theta}}$ given by (129) satisfies

$\displaystyle\|\hat{\bm{\theta}}-{\bm{\theta}}^{\ast}\|_{2}\,\leq\,$
$\displaystyle\frac{3\lambda_{n}}{\kappa_{\cal L}}\Psi({\mathcal{M}})\,.$
(141)

The left-side of (141) equals $\|\hat{\bm{\Delta}}-{\bm{\Delta}}^{\ast}\|_{F}$
while the right-side of (141) equals right-side of (27) using
$\Psi({\mathcal{M}})=\sqrt{s_{n}}$, $\kappa_{\cal L}>\phi_{min}^{\ast}/4$, and
noting that
$\max_{t=1,\cdots,p^{2}}\|{\bm{\theta}}_{Gt}^{\ast}\|_{2}=\max_{\\{k,\ell\\}\in
V\times V}\|({\bm{\Delta}}^{\ast})^{(k\ell)}\|_{F}$. This proves Theorem 2.
$\quad\Box$

## References

  * [1] S.L. Lauritzen, Graphical models. Oxford, UK: Oxford Univ. Press, 1996.
  * [2] J. Whittaker, Graphical Models in Applied Multivariate Statistics. New York: Wiley, 1990.
  * [3] P. Danaher, P. Wang and D.M. Witten, “The joint graphical lasso for inverse covariance estimation across multiple classes,” J. Royal Statistical Society, Series B (Methodological), vol. 76, pp. 373-397, 2014.
  * [4] J. Friedman, T. Hastie and R. Tibshirani, “Sparse inverse covariance estimation with the graphical lasso,” Biostatistics, vol. 9, no. 3, pp. 432-441, July 2008.
  * [5] N. Meinshausen and P. Bühlmann, “High-dimensional graphs and variable selection with the Lasso,” Ann. Statist., vol. 34, no. 3, pp. 1436-1462, 2006.
  * [6] K. Mohan, P. London, M. Fazel, D. Witten and S.I. Lee, “Node-based learning of multiple Gaussian graphical models,” J. Machine Learning Research, vol. 15, pp. 445-488, 2014.
  * [7] Y. Wu, T. Li, X. Liu and L.| Chen, “Differential network inference via the fused D-trace loss with cross variables,” Electronic J. Statistics, vol. 14, pp. 1269-1301, 2020.
  * [8] H. Yuan, R. Xi, C. Chen and M. Deng, “Differential network analysis via lasso penalized D-trace loss,” Biometrika, vol. 104, pp. 755-770, 2017.
  * [9] Z. Tang, Z. Yu and C. Wang, “A fast iteraive algorithm for high-dimensional differential network,” Computational Statistics, vol. 35, pp. 95-109, 2020.
  * [10] B. Zhao, Y.S. Wang and M. Kolar, “Direct estimation of differential functional graphical models,” in Proc. 33rd Conf. Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada, 2019.
  * [11] B. Zhao, Y.S. Wang and M. Kolar, “FuDGE: A method to estimate a functional differential graph in a high-dimensional setting,” J. Machine Learning Research, vol. 23, pp. 1-82, 2022.
  * [12] S.D. Zhao, T.T. Cai and H. Li, “Direct estimation of differential networks,” Biometrika, vol. 101, pp. 253-268, June 2014.
  * [13] E. Belilovsky, G. Varoquaux and M.B. Blaschko, “Hypothesis testing for differences in Gaussian graphical models: Applications to brain connectivity,” Advances in Neural Information Processing Systems (NIPS 2016), vol. 29, Dec. 2016.
  * [14] J. Chiquet, G. Rigaill and M. Sundquist, “A multiattribute Gaussian graphical model for inferring multiscale regulatory networks: an application in breast cancer.” In: Sanguinetti G., Huynh-Thu V. (eds), Gene Regulatory Networks. Methods in Molecular Biology, vol 1883. Humana Press, New York, NY, 2019.
  * [15] M. Kolar, H. Liu and E.P. Xing, “Markov network estimation from multi-attribute data,” in Proc. 30th Intern. Conf. Machine Learning (ICML), Atlanta, GA, 2013.
  * [16] M. Kolar, H. Liu and E.P. Xing, “Graph estimation from multi-attribute data,” J. Machine Learning Research, vol. 15, pp. 1713-1750, 2014.
  * [17] J.K. Tugnait, “Sparse-group lasso for graph learning from multi-attribute data,” IEEE Trans. Signal Process., vol. 69, pp. 1771-1786, 2021. (Corrections, vol. 69, p. 4758, 2021.)
  * [18] G. Marjanovic and V. Solo, “Vector $l_{0}$ sparse conditional independence graphs,” in Proc. IEEE ICASSP 2018, pp. 2731-2735, 2018.
  * [19] Z. Yue, P. Sundaram and V. Solo, “Fast block-sparse estimation for vector networks,” in Proc. IEEE ICASSP 2020, pp. 5505-5509, 2020.
  * [20] P. Sundaram, M. Luessi, M. Bianciardi, S. Stufflebeam, M. Hämäläinen and V. Solo, “Individual resting-state brain networks enabled by massive multivariate conditional mutual information,” IEEE Trans. Med. Imaging, vol. 39, pp. 1957-1966, 2020.
  * [21] Z. Yue and V. Solo, “Comparing vector networks via frequency domain persistent homology,” in Proc. IEEE CDC, pp. 126-131, Dec. 2021.
  * [22] K. Tsai, O. Koyejo, and M. Kolar, “Joint Gaussian graphical model estimation: A survey,” WIREs Computational Statistics, vol. 14, no. 6, pp. e1582, Nov/Dec. 2022.
  * [23] B. Zhang, H. Li, R.B. Riggins, M. Zhan, J. Xuan, Z. Zhang, E.P. Hoffman, R. Clarke and Y. Wang, “Differential dependency network analysis to identify condition-specific topological changes in biological networks,” Bioinformatics, vol. 25, no. 4, pp. 526-532, 2009.
  * [24] P. Xu and Q. Gu, “Semiparametric differential graph models,” in Proc. 30th Conference on Neural Information Processing Systems (NIPS 2016),Barcelona, Spain, 2016.
  * [25] S. Liu, J.A. Quinn, M.U. Gutmann, T. Suzuki and M. Sugiyama, “Direct learning of sparse changes in Markov networks by density ratio estimation,” Neural Computation, 26(6):1169-1197, 2014.
  * [26] S. Liu, T. Suzuki, R. Relator, J. Sese, M. Sugiyama, and K. Fukumizu, “Support consistency of direct sparse-change learning in Markov networks” Annals Statistics, vol. 45, no. 3, pp. 959-990, 2017.
  * [27] S. Liu, K. Fukumizu and T. Suzuki, “Learning sparse structural changes in high-dimensional Markov network: A review on methodologies and theories,” Behaviormetrika, vol. 44, pp. 265-286, 2017.
  * [28] B. Jiang, X. Wang and C. Leng, “A direct approach for sparse quadratic discriminant analysis,” J. Machine Learning Research, vol. 19, pp. 1-37, 2018.
  * [29] T. Zhang and H. Zou, “Sparse precision matrix estimation via lasso penalized D-trace loss,” Biometrika, vol. 101, pp. 103-120, 2014.
  * [30] S.N. Negahban, P. Ravikumar, M.J. Wainwright and B. Yu, “A unified framework for high-dimensional analysis of M-estimators with decomposable regularizers,” Statistical Science, vol. 27, No. 4, pp. 538-557, 2012.
  * [31] P. Ravikumar, M.J. Wainwright, G. Raskutti and B. Yu, “High-dimensional covariance estimation by minimizing $\ell_{1}$-penalized log-determinant divergence,” Electronic J. Statistics, vol. 5, pp. 935-980, 2011.
  * [32] J.K. Tugnait, “Estimation of high-dimensional differential graphs from multi-attribute data,” in Proc. 2023 IEEE Intern. Conf. Acoustics, Speech & Signal Processing (ICASSP 2023), pp. 1-5, Rhodes Island, Greece, June 4-9, 2023.
  * [33] D.S. Tracy and K.G. Jinadasa, “Partitioned Kronecker products of matrices and applications,” Canadian J. Statistics, vol. 17, pp. 107-120, March 1989.
  * [34] S. Liu, “Matrix results on Khatri-Rao and Tracy-Singh products,” Linear Algebra & Its Applications, vol. 289, pp. 267-277, 1999.
  * [35] M. Yuan and Y. Lin, “Model selection and estimation in regression with grouped variables,” J. Royal Statistical Society, Series B (Methodological), vol. 68, pp. 49-67, 2006.
  * [36] J. Friedman, T. Hastie and R. Tibshirani, “A note on the group lasso and a sparse group lasso,” arXiv:1001.0736v1 [math.ST], 5 Jan 2010.
  * [37] N. Simon, J. Friedman, T. Hastie and R. Tibshirani, “A sparse-group lasso,” J. Computational Graphical Statistics, vol. 22, pp. 231-245, 2013.
  * [38] S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1-122, 2010.
  * [39] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 1, pp. 183-202, 2009.
  * [40] C. Leng and C.Y. Tang, “Sparse matrix graphical models,” J. American Statistical Association, vol. 107, pp. 1187-1200, Sep. 2012.
  * [41] T. Tsiligkaridis, A.O. Hero, III, and S. Zhou, “On convergence of Kronecker graphical lasso algorithms,” IEEE Trans. Signal Process., vol. 61, no. 7, pp. 1743-1755, April 2013.
  * [42] X. Lyu, W.W. Sun, Z. Wang, H. Liu, J. Yang and G. Cheng, “Tensor graphical model: Non-convex optimization and statistical inference," IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 42, no. 8, pp. 2024-2037, 1 Aug. 2020.
  * [43] A-L. Barabási and R. Albert, Réka, “Emergence of scaling in random networks,” Science, vol. 286, no. 5439, pp. 509-512, Oct. 1999.
  * [44] S. Lu, J. Kang, W. Gong and D. Towsley, “Complex network comparison using random walks,” in WWW ’14 Companion: Proc. 23rd Intern. Conf. World Wide Web, pp. 727-730, Seoul, Korea, April 2014.
  * [45] S. Zhang, B. Guo, A. Dong, J. He, Z. Xu and S.X. Chen, “Cautionary tales on air-quality improvement in Beijing,” Proc. Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 473, p. 20170457, 2017.
  * [46] W. Chen, F. Wang, G. Xiao, J. Wu and S. Zhang, “Air quality of Beijing and impacts of the new ambient air quality standard,” Atmosphere, vol. 6, pp. 1243-1258, 2015.