\section{Introduction}

In many modern problems, we are interested in studying a network of
entities with multiple attributes rather than a simple univariate
attribute. For example, when an entity represents a person in a social
network, it is widely accepted that the nodal attribute is most
naturally a vector with many personal information including
demographics, interests, and other features, rather than merely a
single attribute, such as a binary vote as assumed in the current
literature of social graph estimation based on Markov random fields
\citep{Banerjee08,kolar08estimating}. In another example, when an
entity represents a gene in a gene regulation network, modern data
acquisition technologies allow researchers to measure the activities
of a single gene in a high-dimensional space, such as an image of the
spatial distribution of the gene expression, or a multi-view snapshot
of the gene activity such as mRNA and protein abundances, rather than
merely a single attribute such as an expression level, which is
assumed in the current literature on gene graph estimation based on
Gaussian graphical models~\citep{peng09partial}. Indeed, it is
somewhat surprising that existing research on graph estimation remains
largely blinded to the analysis of multi-attribute data that are
prevalent and widely studied in the network community. Existing
algorithms and theoretical analysis relies heavily on covariance
selection using graphical lasso, or penalized pseudo-likelihood. They
can not be easily extended to graphs with multi-variate nodal
attributes.

In this paper, we present a study on graph estimation from
multi-attribute data, in an attempt to fill the gap between the
practical needs and existing methodologies from the literature. Under
a Gaussian graphical model, one assumes that a $p$-dimensional random
vector $X \in \RR^p$ follows a multivariate Gaussian distribution with
the mean $\mu$ and covariance matrix $\Sigma$, with each component of
the vector corresponding to a node in the graph.  Based on $n$
independent and identically distributed observations, one can estimate
an undirected graph $G=(V,E)$, where the node set $V$ corresponds to
the $p$ variables, and the edge set $E$ describes the conditional
independence relationships among the variables, that is, variables
$X_a$ and $X_b$ are conditionally independent given all the remaining
variables if $(a,b)\notin E$. Given multi-attribute data, this
approach is clearly invalid, because it naively translates to
estimating one graph per attribute. A subsequent integration of all
such graphs to a summary graph on the entire dataset may lead to
unclear statistical interpretation.

We consider the following new setting for estimating a multi-attribute
graph.  In this setting, we consider a ``stacked'' long random vector
$X=(X^T_1,...,X^T_p)^T$ where $X_1 \in \RR^{k_1}, \ldots, X_p
\in \RR^{k_p}$ are themselves random vectors that jointly follow a
multivariate Gaussian distribution with mean $\mu = (\mu_1, \ldots,
\mu_p)^T$ and covariance matrix $\Sigma^*$, which is partitioned as 
\begin{equation}
  \label{eq:vector_model}
\Sigma^* =   \left(
    \begin{array}{cccc}
     \Sigma_{11}^*  & \cdots & \Sigma_{1p}^* \\
     \vdots &   \ddots & \vdots \\
     \Sigma_{p1}^* & \cdots   & \Sigma_{pp}^* \\
    \end{array}
  \right),
\end{equation}
with $\Sigma^*_{ij} = \Cov(X_i, X_j)$.  Without loss of generality, we
assume $\mu=0$.  Let $G = (V,E)$ be a graph with the vertex set $V =
\{1,\ldots, p\}$ and the set of edges $E \subseteq V \times V$ that
encodes the conditional independence relationships among $(X_a)_{a\in
  V}$. That is, each node $a \in V$ of the graph $G$ corresponds to
the random vector $X_a$ and there is no edge between nodes $a$ and $b$
in the graph if and only if $X_a$ is conditionally independent of
$X_b$ given all the vectors corresponding to the remaining nodes,
$X_{\neg{ab}} = \{X_c\ :\ c \in V\bks\{a,b\}\}$. Such a graph is also
known as a Markov network (of Markov graph), which we shall emphasize
in this paper to compare with an alternative graph over $V$ known as
the association network, which is based on pairwise marginal
independence. Conditional independence can be read from the inverse of
the covariance matrix, as the block corresponding to $X_a$ and $X_b$
will be equal to zero. Let $\Dcal_n = \left\{x_i\right\}_{i=1}^n$ be a
sample of $n$ independent and identically distributed vectors drawn
from $\mathrm{N}(0, \Sigma)$. For a vector $x_i$, we denote $x_{i,a}
\in \RR^{k_a}$ the component corresponding to the node $a \in V$. Our
goal is to estimate the structure of the graph $G$ from the sample
$\Dcal_n$. Note that we allow for different nodes to have different
number of attributes, which is useful in many applications, e.g., when
a node represents a gene pathway in a regulatory network.


Using the standard Gaussian graphical model for univariate nodal
observations, one can estimate the Markov graph for each attribute
individually by estimating the sparsity pattern of the precision
matrix $\Omega = \Sigma^{-1}$.  This is also known as the covariance
selection problem \citep{Dempster72}. For high dimensional problems,
\citet{Meinshausen06} propose a parallel Lasso approach for estimating
Gaussian graphical models by solving a collection of sparse regression
problems. This procedure can be viewed as a pseudo-likelihood based
method. In contrast, \citet{Banerjee08}, \citet{Yuan07}, and
\citet{Friedman08} take a penalized likelihood approach to estimate
the sparse precision matrix $\Omega$. To reduce estimation bias,
\citet{Lam09}, \citet{Jalali12}, and \citet{Shen12} developed the
non-concave penalties to penalize the likelihood function. More
recently, \citet{Yuan10} and \citet{Cai11a} proposed the graphical
Dantzig selector and CLIME, which can be solved by linear programming
and are more amenable to theoretical analysis than the penalized
likelihood approach. Under certain regularity conditions, these
methods have proven to be graph estimation consistent
\citep{Ravikumar11,Yuan10,Cai11a} and scalable software packages, such
as $\mathrm{glasso}$ and $\mathrm{huge}$, were developed to implement
these algorithms \citep{Zhao12b}. However, in the case of
multi-attribute data, it is not clear how to combine estimated graphs
to obtain a single Markov graph reflecting the structure of the
underlying complex system. This is especially the case when nodes in
the graph contain different number of attributes.

In a previous work, \citet{katenka2011multi} proposed a method for
estimating association networks from multi-attribute data using
canonical correlation as a dependence measure between two groups of
attributes. However, association networks are known to confound the
direct interactions with indirect ones as they only represent marginal
associations.  In contrast, we develop a method based on partial
canonical correlation, which give rise to a Markov graph that is
better suited for separating direct interactions from indirect
confounders.  Our work is related to the literature on simultaneous
estimation of multiple Gaussian graphical models under a multi-task
setting
\citep{Guo:09,varoquaux2010brain,honorio2010multi,Chiquet2011inferring,danaher2011joint}.
However, the model given in \eqref{eq:vector_model} is different from
models considered in various multi-task settings and the optimization
algorithms developed in the multi-task literature do not extend to
handle the optimization problem given in our setting.


Unlike the standard procedures for estimating the structure of
Gaussian graphical models (e.g., neighborhood selection
\citep{Meinshausen06} or glasso \citep{Friedman08}), which infer the
partial correlations between pairs of multi-attribute nodes, our
proposed method estimates the {partial canonical correlations} between
pairs of nodes, which leads to a graph estimator over multi-attribute
nodes that bears the same probabilistic independence interpretations
as that of the graph from Gaussian graphical model over univariate
nodes. Under this new framework, the contributions of this paper
include: (i) Computationally, an efficient algorithm is provided to
estimate the multi-attribute Markov graphs; (ii) Theoretically, we
provide sufficient conditions which guarantee consistent graph
recovery; and (iii) Empirically, we apply our procedure to uncover
gene regulatory networks from gene and protein profiles, and to
uncover brain connectivity graph from functional magnetic resonance
imaging data. 

\section{Methodology}

In this section, we propose to estimate the graph by estimating
non-zero partial canonical correlation between the nodes. This leads
to a penalized maximum likelihood objective, for which we develop an
efficient optimization procedure. 

\subsection{Preliminaries}
\label{sec:prelim}

Let $X_a$ and $X_b$ be two multivariate random vectors.  Canonical
correlation is defined between $X_a$ and $X_b$ as
\begin{equation*}
\rho_c(X_a, X_b) = \max_{\ub \in \RR^{k_a}, \vb\in\RR^{k_b}} \Cor(\ub^TX_a, \vb^TX_b).
\end{equation*}
That is, computing canonical correlation between $X_a$ and $X_b$ is
equivalent to maximizing the correlation between two linear
combinations $\ub^TX_a$ and $\vb^TX_b$ with respect to vectors $\ub$
and $\vb$. Canonical correlation can be used to measure association
strength between two nodes with multi-attribute observations. For
example, in \citet{katenka2011multi}, a graph is estimated from
multi-attribute nodal observations by elementwise thresholding the
canonical correlation matrix between nodes, but such a graph estimator
may confound the direct interactions with indirect ones.

In this paper, we exploit the partial canonical correlation to
estimate a graph from multi-attribute nodal observations. A graph
is going to be formed by connecting nodes with non-zero partial
canonical correlation. Let $\hat\Ab = \argmin\ \EE\rbr{\norm{X_a -
  \Ab X_{\neg{ab}}}_2^2}$ and $\hat\Bb = \argmin\ \EE\rbr{\norm{X_b -
  \Bb X_{\neg{ab}}}_2^2}$, then the partial canonical correlation
between $X_a$ and $X_b$ is defined as 
\begin{equation}
\label{eq:partial_correlation}
  \begin{aligned}
\rho_c(X_a, X_b; X_{\neg{ab}}) 
= \max_{\ub\in\RR^{k_a}, \vb\in\RR^{k_b}} \Cor\{\ub^T(X_a
- \hat\Ab X_{\neg{ab}}), \vb^T(X_b - \hat\Bb X_{\neg{ab}})\},    
  \end{aligned}
\end{equation}
 that is, the partial canonical correlation between $X_a$ and $X_b$
is equal to the canonical correlation between the residual vectors of
$X_a$ and $X_b$ after the effect of $X_{\neg{ab}}$ is
removed \citep{rao_partial_1969}.

Let $\Omega^*$ denote the precision matrix under the model in
\eqref{eq:vector_model}. Using standard results for the multivariate
Gaussian distribution (see also Equation (7) in
\citet{rao_partial_1969}), a straightforward calculation shows
that\footnote{Calculation given in Appendix~\ref{sec:appendix_eq_3}}
\begin{equation}
  \label{eq:comp_cov}
  \rho_c(X_a, X_b; X_{\neg{ab}}) \neq 0
  \quad 
  \text{if and only if}
  \quad
   \max_{\ub\in\RR^{k_a}, \vb\in\RR^{k_b}}\ \ub^T\Omega_{ab}^*\vb \neq 0.
\end{equation}
This implies that estimating whether the partial canonical correlation
is zero or not can be done by estimating whether a block of the
precision matrix is zero or not. Furthermore, under the model in
\eqref{eq:vector_model}, vectors $X_a$ and $X_b$ are conditionally
independent given $X_{\neg{ab}}$ if and only if partial canonical
correlation is zero. A graph built on this type of inter-nodal
relationship is known as a {Markov graph}, as it captures both
local and global Markov properties over all arbitrary subsets of nodes
in the graph, even though the graph is built based on pairwise
conditional independence properties. In
$\S$\ref{sec:estimation}, we use the above observations to design
an algorithm that estimates the non-zero partial canonical correlation
between nodes from data $\Dcal_n$ using the penalized maximum
likelihood estimation of the precision matrix.

Based on the relationship given in   \eqref{eq:comp_cov}, we can
motivate an alternative method for estimating the non-zero partial
canonical correlation. Let $\bar{a} = \{b\ :\ b\in V\bks\{a\}\}$
denote the set of all nodes minus the node $a$. Then
\begin{equation*}
  \label{eq:population_block_regression}
  \EE\rbr{X_a\mid X_{\bar a} = x_{\bar a}} =
  \Sigma_{a,\ia}^*\Sigma_{\ia,\ia}^{*,-1}
  x_{\bar a}.
\end{equation*}
Since $\Omega^{*}_{a,\ia} =
-(\Sigma_{aa}^*-\Sigma_{a,\ia}^*\Sigma_{\ia,\ia}^{*,-1}\Sigma_{\ia,a}^*)^{-1}
\Sigma_{a,\ia}^*\Sigma_{\ia,\ia}^{*,-1}$, we observe that a zero block
$\Omega_{ab}$ can be identified from the regression coefficients when
each component of $X_a$ is regressed on $X_{\bar a}$. We do not build
an estimation procedure around this observation, however, we note that
this relationship shows how one would develop a regression
based analogue of the work presented in \citet{katenka2011multi}. 

\subsection{Penalized Log-Likelihood Optimization}
\label{sec:estimation}

Based on the data $\Dcal_n$, we propose to minimize the penalized
negative Gaussian log-likelihood under the model in
\eqref{eq:vector_model},
\begin{equation}
  \label{eq:max_ll_opt}
  \min_{\Omegab \succ \zero}\ \Bigl\{\tr \Sbb \Omegab - \log|\Omegab|
    + \lambda \sum_{a, b} \norm{\Omegab_{ab}}_F\Bigr\}
\end{equation}
where $\Sbb = n^{-1}\sum_{i=1}^n\xb_i\xb_i^T$ is the sample covariance
matrix, $\norm{\Omegab_{ab}}_F)$ denotes the Frobenius norm of
$\Omegab_{ab}$ and $\lambda$ is a user defined parameter that controls
the sparsity of the solution $\hat\Omega$. The Frobenius norm penalty
encourages blocks of the precision matrix to be equal to zero, similar
to the way that the $\ell_2$ penalty is used in the group Lasso
\citep{Yuan06model}.  Here we assume that the same number of samples
is available per attribute. However, the same method can be used in
cases when some samples are obtained on a subset of
attributes. Indeed, we can simply estimate each element of the matrix
$\Sbb$ from available samples, treating non-measured attributes as
missing completely at random \citep[for more details 
see][]{kolar12missing}.
 

The dual problem to \eqref{eq:max_ll_opt} is
\begin{equation}
  \label{eq:max_ll_opt_dual}
    \max_{\Sigmab} \sum_{j \in V} k_j + \log|\Sigmab| \qquad
      \text{subject to } \qquad
      \max_{a,b}\ \norm{\Sbb_{ab} - \Sigmab_{ab}}_F \leq \lambda,
\end{equation}
where $\Sigmab$ is the dual variable to $\Omegab$ and $|\Sigma|$
denotes the determinant of $\Sigma$. Note that the primal problem
gives us an estimate of the precision matrix, while the dual problem
estimates the covariance matrix. The proposed optimization procedure,
described below, will simultaneously estimate the precision matrix and
covariance matrix, without explicitly performing an expensive matrix
inversion.

We propose to optimize the objective function in \eqref{eq:max_ll_opt}
using an inexact block coordinate descent procedure, inspired by
\citet{mazumder11flexible}. The block coordinate descent is an
iterative procedure that operates on a block of rows and columns while
keeping the other rows and columns fixed. We write
\begin{equation*}
  \Omegab = \left(
  \begin{array}{cc}
    \Omegab_{aa} & \Omegab_{a,\ia} \\
    \Omegab_{\ia,a} & \Omegab_{\ia,\ia}\\
  \end{array}
  \right),
\quad
  \Sigmab = \left(
  \begin{array}{cc}
    \Sigmab_{aa} & \Sigmab_{a,\ia} \\
    \Sigmab_{\ia,a} & \Sigmab_{\ia,\ia}\\
  \end{array}
  \right),
\quad
  \Sbb = \left(
  \begin{array}{cc}
    \Sbb_{aa} & \Sbb_{a,\ia} \\
    \Sbb_{\ia,a} & \Sbb_{\ia,\ia}\\
  \end{array}
  \right)
\end{equation*}
and suppose that $(\tilde\Omegab, \tilde\Sigmab)$ are the current
estimates of the precision matrix and covariance matrix. With the
above block partition, we have $\log|\Omegab| =
\log(\Omegab_{\ia,\ia}) + \log(\Omegab_{aa} -
\Omegab_{a,\ia}(\Omegab_{\ia,\ia})^{-1}\Omegab_{\ia,a})$.  In the next
iteration, $\hat \Omegab$ is of the form
\begin{equation*}
  \hat\Omegab = \tilde\Omegab + \left(
  \begin{array}{cc}
    \Deltab_{aa} & \Deltab_{a,\ia} \\
    \Deltab_{\ia,a} & \zero
  \end{array}
\right) =
\left(
  \begin{array}{cc}
    \hat\Omegab_{aa} & \hat\Omegab_{a,\ia} \\
    \hat\Omegab_{\ia,a} & \tilde\Omegab_{\ia,\ia} \\
  \end{array}
\right)
\end{equation*}
and is obtained by minimizing
\begin{equation}
  \label{eq:partial_min}
  \tr\Sbb_{aa}\Omegab_{aa}
\! + \! 2\tr\Sbb_{a,\ia}\Omegab_{\ia,a}
  - \log|\Omegab_{aa} - \Omegab_{a,\ia}(\tilde\Omegab_{\ia,\ia})^{-1}\Omegab_{\ia,a}|
  + \lambda\norm{\Omegab_{aa}}_F + 2\lambda\sum_{b \neq
    a}\norm{\Omegab_{ab}}_F.
\end{equation}
Exact minimization over the variables $\Omegab_{aa}$ and
$\Omegab_{a,\ia}$ at each iteration of the block coordinate descent
procedure can be computationally expensive. Therefore, we propose to
update $\Omegab_{aa}$ and $\Omegab_{a,\ia}$ using one generalized
gradient step update (see \citet{beck09fast}) in each iteration.  Note
that the objective function in \eqref{eq:partial_min} is a sum of a
smooth convex function and a non-smooth convex penalty so that the
gradient descent method cannot be directly applied. Given a step size
$t$, generalized gradient descent optimizes a quadratic approximation
of the objective at the current iterate $\tilde\Omegab$, which results
in the following two updates
\begin{align}
  \label{eq:update_aa}
  \hat\Omegab_{aa} &= \argmin_{\Omegab_{aa}} \Big\{
  \tr(\Sbb_{aa} - \tilde\Sigmab_{aa})\Omegab_{aa} +
  \frac{1}{2t}\norm{\Omegab_{aa} - \tilde\Omegab_{aa}}_F^2 +
  \lambda \norm{\Omegab_{aa}}_F\Bigr\},
  \quad\text{ and} \\
  \label{eq:update_ab}
  \hat\Omegab_{ab} &= \argmin_{\Omegab_{ab}}
  \Bigl\{ \tr(\Sbb_{ab} - \tilde\Sigmab_{ab})\Omegab_{ba} +
  \frac{1}{2t}\norm{\Omegab_{ab} - \tilde\Omegab_{ab}}_F^2 +
  \lambda \norm{\Omegab_{ab}}_F\Bigr\},
  \quad \forall b \in {\bar a}.
\end{align}
If the resulting estimator $\hat\Omegab$ is not positive definite or
the update does not decrease the objective, we halve the step size $t$
and find a new update.  Once the update of the precision matrix
$\hat\Omegab$ is obtained, we update the covariance matrix
$\hat\Sigmab$.  Updates to the precision and covariance matrices can
be found efficiently, without performing expensive matrix inversion,
as we show in Appendix~\ref{appendix:a}
(see~\eqref{eq:update_aa_final}--\eqref{eq:update_sigma}).  Combining
all three steps we get the following algorithm:
\begin{enumerate}
\item
Set the initial estimator $\tilde\Omegab = \diag(\Sbb)$ and
$\tilde\Sigmab = \tilde\Omegab^{-1}$. Set the step size $t=1$.
\item
For each $a \in V$ perform the following:
\begin{tabbing}
\qquad Update $\hat\Omegab$ using 
\eqref{eq:update_aa_final} and \eqref{eq:update_ab_final}.\\
\qquad If $\hat\Omegab$ is not positive definite,
  set $t \leftarrow t/2$ and repeat the update.\\
\qquad Update $\hat\Sigmab$ using \eqref{eq:update_sigma}.
\end{tabbing}
\item
  Repeat Step 2 until the duality gap 
  \[
 \Bigl| \tr (\Sbb \hat \Omegab) - \log|\hat\Omegab|
  + \lambda \sum_{a, b} \norm{\hat\Omegab_{ab}}_F -
  \sum_{j \in V} k_j - \log|\Sigmab|\Bigr| \leq \epsilon,  
  \]
  where $\epsilon$ is a prefixed  precision parameter (for example, $\epsilon=10^{-3}$).
\end{enumerate}
Finally, we form a graph $\hat G = (V, \hat E)$ by connecting nodes
with $\norm{\hat\Omegab_{ab}}_F \neq 0$.

Computational complexity of the procedure is given in
Appendix~\ref{appendix:b}.  Convergence of the above described
procedure to the unique minimum of the objective function in
\eqref{eq:max_ll_opt} does not follow from the standard results on the
block coordinate descent algorithm \citep{tseng01convergence} for two
reasons. First, the minimization problem in \eqref{eq:partial_min} is
not solved exactly at each iteration, since we only update
$\Omegab_{aa}$ and $\Omegab_{a,\ia}$ using one generalized gradient step
update in each iteration. Second, the blocks of variables, over which
the optimization is done at each iteration, are not completely
separable between iterations due to the symmetry of the problem. The
proof of the following convergence result is given in
Appendix~\ref{appendix:c}.

\begin{lemma}
  \label{lem:convergence}
  For every value of $\lambda > 0$, the above described algorithm
  produces a sequence of estimates $\cbr{\tilde\Omegab^{(t)}}_{t \geq
    1}$ of the precision matrix that monotonically decrease the
  objective values given in \eqref{eq:max_ll_opt}. Every element of
  this sequence is positive definite and the sequence converges to the
  unique minimizer $\hat \Omegab$ of \eqref{eq:max_ll_opt}.
\end{lemma}

\subsection{Efficient Identification of Connected Components}

When the target graph $\hat G$ is composed of smaller, disconnected
components, the solution to the problem in \eqref{eq:max_ll_opt} is
block diagonal (possibly after permuting the node indices) and can be
obtained by solving smaller optimization problems. That is, the
minimizer $\hat\Omegab$ can be obtained by solving
\eqref{eq:max_ll_opt} for each connected component independently,
resulting in massive computational gains. We give necessary and
sufficient condition for the solution $\hat\Omegab$ of
\eqref{eq:max_ll_opt} to be block-diagonal, which can be easily
checked by inspecting the empirical covariance matrix $\Sbb$.  

Our first result follows immediately from the Karush-Kuhn-Tucker
conditions for the optimization problem \eqref{eq:max_ll_opt} and
states that if $\hat\Omegab$ is block-diagonal, then it can be
obtained by solving a sequence of smaller optimization problems.
\begin{lemma}
  If the solution to \eqref{eq:max_ll_opt} takes the form 
  $\hat\Omega={\rm diag}(\hat\Omega_1,
  \hat\Omega_2,\ldots,\hat\Omega_l)$, that is, $\hat\Omega$ is a block
  diagonal matrix with the diagonal blocks $\hat\Omega_1, \ldots,
  \hat\Omega_l$, then it can be obtained by solving
  \begin{equation*}
  \min_{\Omegab_{l'} \succ \zero}\ \Bigl\{ \tr \Sbb_{l'} \Omegab_{l'} - \log|\Omegab_{l'}|
    + \lambda \sum_{a, b} \norm{\Omegab_{ab}}_F\Bigr\}
  \end{equation*}
  separately for each $l'=1,\ldots,l$, where $\Sbb_{l'}$ are
  submatrices of $\Sbb$ corresponding to $\Omegab_{l'}$.
\end{lemma}
Next, we describe how to identify diagonal blocks of
$\hat\Omegab$. Let $\Pcal = \{P_1, P_2, \ldots, P_l\}$ be a partition
of the set $V$ and assume that the nodes of the graph are ordered in
a way that if $a \in P_j$, $b \in P_{j'}$, $j < j'$, then $a < b$.
The following lemma states that the blocks of $\hat\Omega$ can be
obtained from the blocks of a thresholded sample covariance matrix.
\begin{lemma}
  \label{lem:nec_suff_block_diag}
  A necessary and sufficient conditions for $\hat\Omegab$ to be block
  diagonal with blocks $P_1, P_2, \ldots, P_l$ is that
  $\norm{\Sbb_{ab}}_F \leq \lambda$ for all $a \in P_j$, $b \in
  P_{j'}$, $j \neq j'$.
\end{lemma}

Blocks $P_1, P_2, \ldots, P_l$ can be identified by forming a $p\times
p$ matrix $\Qb$ with elements $q_{ab} = \ind\{\norm{\Sbb_{ab}}_F >
\lambda\}$ and computing connected components of the graph with
adjacency matrix $\Qb$. The lemma states also that given two penalty
parameters $\lambda_1 < \lambda_2$, the set of unconnected nodes with
penalty parameter $\lambda_1$ is a subset of unconnected nodes with
penalty parameter $\lambda_2$. The simple check above allows us to
estimate graphs on datasets with large number of nodes, if we are
interested in graphs with small number of edges. However, this is
often the case when the graphs are used for exploration and
interpretation of complex systems.
Lemma~\ref{lem:nec_suff_block_diag} is related to existing results
established for speeding-up computation when learning single and
multiple Gaussian graphical models \citep{witten2011new,
  mazumder2011exact, danaher2011joint}. Each condition is different,
since the methods optimize different objective functions.

\section{Consistent Graph Identification}
\label{sec:theory}

In this section, we provide theoretical analysis of the estimator
described in $\S$\ref{sec:estimation}. In particular, we provide
sufficient conditions for consistent graph recovery.  For simplicity
of presentation, we assume that $k_a = k$, for all $a \in V$, that is,
we assume that the same number of attributes is observed for each
node. For each $a=1,\ldots,kp$, we assume that
$(\sigma_{aa}^*)^{-1/2}X_a$ is sub-Gaussian with parameter $\gamma$,
where $\sigma_{aa}^*$ is the $a$th diagonal element of
$\Sigmab^*$. Recall that $Z$ is a sub-Gaussian random variable if
there exists a constant $\sigma \in (0, \infty)$ such that
\begin{equation*}
  \EE\rbr{\exp(tZ)} \leq \exp(\sigma^2t^2),\ \text{for all
  } t \in \RR.
\end{equation*}

Our assumptions involve the Hessian of the function $f(\Ab) =
\tr\Sbb\Ab - \log|\Ab|$ evaluated at the true $\Omegab^*$, $\Hcal =
\Hcal(\Omega^{*}) = (\Omega^{*})^{-1}\otimes(\Omega^{*})^{-1} \in
\RR^{(pk)^2\times(pk)^2}$, with $\otimes$ denoting the Kronecker
product, and the true covariance matrix $\Sigmab^*$. The Hessian and
the covariance matrix can be thought of as block matrices with blocks
of size $k^2 \times k^2$ and $k \times k$, respectively. We will make
use of the operator $\Ccal(\cdot)$ that operates on these block
matrices and outputs a smaller matrix with elements that equal to the
Frobenius norm of the original blocks. For example, $\Ccal(\Sigmab^*)
\in \RR^{p\times p}$ with elements $\Ccal(\Sigmab^*)_{ab} =
\norm{\Sigmab^*_{ab}}_F$. Let $\Tcal = \{(a,b): \norm{\Omegab_{ab}}_F
\neq 0\}$ and $\Ncal = \{(a,b): \norm{\Omegab_{ab}}_F = 0\}$. With
this notation introduced, we assume that the following
{irrepresentable} condition holds. There exists a constant $\alpha \in
[0, 1)$ such that
\begin{equation}
\label{eq:assum:irrepresentable}
\opnorm{
\Ccal\left(
  \Hcal_{\Ncal\Tcal}(\Hcal_{\Tcal\Tcal})^{-1}
\right)}{\infty} \leq 1 - \alpha,
\end{equation}
where $\opnorm{A}{\infty} = \max_i \sum_j|A_{ij}|$.
We will also need the following quantities to specify the results $
\kappa_{\Sigmab^*} = \opnorm{\Ccal(\Sigmab^*)}{\infty}$ and
$\kappa_{\Hcal} =
\opnorm{\Ccal(\Hcal_{\Tcal\Tcal}^{-1})}{\infty}$. These conditions
extend the conditions specified in \citet{Ravikumar11} needed for
estimating graphs from single attribute observations. 

We have the following result that provides sufficient conditions for
the exact recovery of the graph.
\begin{proposition}
\label{prop:sparsistency}
Let $\tau > 2$. We set
the penalty parameter $\lambda$ in \eqref{eq:max_ll_opt} as 
\[
\lambda = 8k\alpha^{-1}
\rbr{128(1+4\gamma^2)^2(\max_a(\sigma_{aa}^*)^2)n^{-1}(2\log(2k) +
  \tau\log(p))}^{1/2}.
\]
If $
  n > C_1s^2k^2(1+8\alpha^{-1})^2(\tau\log p + \log 4 + 2\log k),$
where $s$ is the maximal degree of nodes in $G$, 
$
C_1 =
(48\sqrt{2}(1+4\gamma^2)(\max_a \sigma_{aa}^*)
\max(\kappa_{\Sigmab^*}\kappa_{\Hcal},
\kappa_{\Sigmab^*}^3\kappa_{\Hcal}^2))^2
$
 and
\[
\min_{(a,b)\in\Tcal, a\neq b} \norm{\Omegab_{ab}}_F > 
16\sqrt{2}(1+4\gamma^2)(\max_a \sigma_{aa}^*)
          (1+8\alpha^{-1})\kappa_\Hcal k
          \rbr{\frac{\tau\log p + \log 4 + 2\log k}{n}}^{1/2},
\]
then $ \PP\rbr{\hat G = G} \geq 1 - p^{2-\tau}$.
\end{proposition}

The proof of Proposition~\ref{prop:sparsistency} is given in
Appendix~\ref{appendix:c}. We extend the proof of \citet{Ravikumar11}
to accommodate the Frobenius norm penalty on blocks of the precision
matrix. This proposition specifies the sufficient sample size and a
lower bound on the Frobenius norm of the off-diagonal blocks needed
for recovery of the unknown graph.  Under these conditions and
correctly specified tuning parameter $\lambda$, the solution to the
optimization problem in \eqref{eq:max_ll_opt} correctly recovers the
graph with high probability. In practice, one needs to choose the
tuning parameter in a data dependent way. For example, using the
Bayesian information criterion.  Even though our theoretical analysis
obtains the same rate of convergence as that of \citet{Ravikumar11},
our method has a significantly improved finite-sample performance
(More details will be provided in $\S$\ref{sec:simulation}.). It
remains an open question whether the sample size requirement can be
improved as in the case of group Lasso \citep[see, for
example,][]{lounici10oracle}.  The analysis of \cite{lounici10oracle}
relies heavily on the special structure of the least squares
regression. Hence, their method does not carry over to the more
complicated objective function as in \eqref{eq:max_ll_opt}.


\section{Interpreting Edges}
\label{sec:interpretation}

We propose a post-processing step that will allow us to quantify the
strength of links identified by the method proposed in
$\S$\ref{sec:estimation}, as well as identify important attributes that
contribute to the existence of links.

For any two nodes $a$ and $b$ for which $\Omegab_{ab}\neq0$, we define
$\Ncal(a,b) = \{ c \in V \bks \{a,b\}\ :\ \Omegab_{ac} \neq 0 \text{
  or } \Omegab_{bc} \neq 0 \}$, which is the Markov blanket for the set
of nodes $\{ \Xb_a, \Xb_b \}$. Note that the conditional distribution
of $(\Xb_a^T, \Xb_b^T)^T$ given $\Xb_{\neg{ab}}$ is equal to the
conditional distribution of $(\Xb_a^T, \Xb_b^T)^T$ given
$\Xb_{\Ncal(a,b)}$.  Now,
\begin{equation*}
  \begin{aligned}
\rho_c(\Xb_a, \Xb_b; \Xb_{\neg{ab}}) 
&= \rho_c(\Xb_a, \Xb_b; \Xb_{\Ncal(a,b)})\\
&= \max_{\wb_a\in\RR^{k_a}, \wb_b\in\RR^{k_b}} \Cor(\ub^T(\Xb_a
- \tilde\Ab\Xb_{\Ncal(a,b)}), \vb^T(\Xb_b - \tilde\Bb\Xb_{\Ncal(a,b)})),     
  \end{aligned}
\end{equation*}
where $\tilde\Ab = \argmin\ \EE\rbr{\norm{\Xb_a -
  \Ab\Xb_{\Ncal(a,b)}}_2^2}$ and $\tilde\Bb = \argmin\ \EE\rbr{\norm{\Xb_b
  - \Bb\Xb_{\Ncal(a,b)}}_2^2}$. Let $\bar\Sigmab(a,b) = \Var(\Xb_a,
\Xb_b\mid\Xb_{\Ncal(a,b)})$. Now we can express the partial canonical
correlation as
\[
\rho_c(\Xb_a, \Xb_b; \Xb_{\Ncal(a,b)}) = 
\max_{\wb_a\in\RR^{k_a}, \wb_b\in\RR^{k_a} }\ 
\frac{\wb_a^T\bar\Sigma_{ab}\wb_b}
{\rbr{\wb_a^T\bar\Sigma_{aa}\wb_a}^{1/2}\rbr{\wb_b^T\bar\Sigma_{bb}\wb_b}^{1/2}}
\]
where 
\[
\bar\Sigmab(a,b) = \left(
\begin{array}{cc}
\bar\Sigmab_{aa} & \bar\Sigmab_{ab} \\
\bar\Sigmab_{ba} & \bar\Sigmab_{bb}
\end{array}
\right).
\]
The weight vectors $\wb_a$ and $\wb_b$ can be easily found by solving 
the system of eigenvalue equations
\begin{equation}
  \label{eq:eigenvalue_system}
\left\{
\begin{array}{c}
  \bar\Sigmab_{aa}^{-1}\bar\Sigmab_{ab}\bar\Sigmab_{bb}^{-1}\bar\Sigmab_{ba}\wb_a
  = \phi^2\wb_a \\
  \bar\Sigmab_{bb}^{-1}\bar\Sigmab_{ba}\bar\Sigmab_{aa}^{-1}\bar\Sigmab_{ab}\wb_b
  = \phi^2\wb_b
\end{array}
\right.
\end{equation}
with $\wb_a$ and $\wb_b$ being the vectors that correspond to the
maximum eigenvalue $\phi^2$. Furthermore, we have $\rho_c(\Xb_a,
\Xb_b; \Xb_{\Ncal(a,b)}) = \phi$. Following \citet{katenka2011multi},
the weights $\wb_a$, $\wb_b$ can be used to access the relative
contribution of each attribute to the edge between the nodes $a$ and
$b$. In particular, the weight $(w_{a,i})^2$ characterizes the
relative contribution of the $i$th attribute of node $a$ to
$\rho_c(\Xb_a, \Xb_b; \Xb_{\Ncal(a,b)})$.

Given an estimate $\hat\Ncal(a, b) = \{ c \in V \bks \{a,b\}\ :\
\hat\Omegab_{ac} \neq 0 \text{ or } \hat\Omegab_{bc} \neq 0 \}$ of the
Markov blanket $\Ncal(a, b)$, we form the residual vectors
\[
\rb_{i, a} = \xb_{i,a} - \check{\Ab}\xb_{i, \hat \Ncal(a,b)},
\qquad
\rb_{i, b} = \xb_{i,b} - \check{\Bb}\xb_{i, \hat \Ncal(a,b)},
\]
where $\check{\Ab}$ and $\check{\Bb}$ are the least square estimators
of ${\tilde{\Ab}}$ and ${\tilde{\Bb}}$. Given the residuals, we form
$\check\Sigmab(a,b)$, the empirical version of the matrix $\bar
\Sigmab(a,b)$, by setting
\[
{\check \Sigma}_{aa} = {\Cor}\rbr{\{\rb_{i,a}\}_{i \in [n]}},\quad
{\check \Sigma}_{bb} = {\Cor}\rbr{\{\rb_{i,b}\}_{i \in [n]}},\quad
{\check \Sigma}_{ab} = {\Cor}\rbr{\{\rb_{i,a} \}_{i \in [n]}, \{\rb_{i,a} \}_{i \in [n]}}.
\]
Now, solving the eigenvalue system in \eqref{eq:eigenvalue_system}
will give us estimates of the vectors $\wb_a$, $\wb_b$ and the partial
canonical correlation.

Note that we have described a way to interpret the elements of the
off-diagonal blocks in the estimated precision matrix. The elements of
the diagonal blocks, which correspond to coefficients between
attributes of the same node, can still be interpreted by their
relationship to the partial correlation coefficients.

\section{Simulation Studies}
\label{sec:simulation}

In this section, we perform a set of simulation studies to illustrate
finite sample performance of our method. We demonstrate that the
scalings of $(n,p,s)$ predicted by the theory are sharp. Furthermore,
we compare against three other methods: 1) a method that uses
the ${\rm glasso}$ first to estimate one graph over each of the $k$
individual attributes and then creates an edge in the resulting graph
if an edge appears in at least one of the single attribute graphs, 2)
the method of \cite{Guo:09} and 3) the method of
\cite{danaher2011joint}.  We have also tried applying the ${\rm
  glasso}$ to estimate the precision matrix for the model in
\eqref{eq:vector_model} and then post-processing it, so that an edge
appears in the resulting graph if the corresponding block of the
estimated precision matrix is non-zero. The result of this method is
worse compared to the first baseline, so we do not report it here.

All the methods above require setting one or two tuning parameters
that control the sparsity of the estimated graph. We select these
tuning parameters by minimizing the Bayesian information criterion,
which balances the goodness of fit of the model and its complexity,
over a grid of parameter values. For our multi-attribute method, the
Bayesian information criterion takes the following form
\[
{\rm BIC}(\lambda) = \tr(\Sbb\hat\Omegab) - \log|\hat\Omegab| 
+ \sum_{a < b} \ind\{\hat\Omegab_{ab} \neq \zero\}k_ak_b\log(n).
\]
Other methods for selecting tuning parameters are possible, like
minimization of cross-validation or Akaike information criterion.
However, these methods tend to select models that are too dense.

Theoretical results given in $\S$\ref{sec:theory} characterize the
sample size needed for consistent recovery of the underlying graph.
In particular, Proposition~\ref{prop:sparsistency} suggests that we
need $n = \theta s^2k^2\log(pk)$ samples to estimate the graph
structure consistently, for some $\theta > 0$. Therefore, if we plot
the hamming distance between the true and recovered graph
against $\theta$, we expect the curves to reach zero distance for
different problem sizes at a same point. We verify this on randomly
generated chain and nearest-neighbors graphs.

\begin{figure}[p]
  \centering

  \begin{subfigure}[b]{0.95\textwidth}
    \centering
    \includegraphics[width=\textwidth]{plots/glassoNormal}
    \caption{glasso procedure}
    \label{fig:normal:glasso}
  \end{subfigure}%

 \begin{subfigure}[b]{0.95\textwidth}
    \centering
    \includegraphics[width=\textwidth]{plots/multiTaskNormal}
    \caption{Procedure of \cite{danaher2011joint}}
    \label{fig:normal:mt}
  \end{subfigure}%

 \begin{subfigure}[b]{0.95\textwidth}
    \centering
    \includegraphics[width=\textwidth]{plots/guoNormal}
    \caption{Procedure of \cite{Guo:09}}
    \label{fig:normal:guo}
  \end{subfigure}%

 \begin{subfigure}[b]{0.95\textwidth}
    \centering
    \includegraphics[width=\textwidth]{plots/multiAttNormal}
    \caption{Multi-attribute procedure}
    \label{fig:normal:ma}
  \end{subfigure}%
  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent
    runs. Off-diagonal blocks are full matrices.}
  \label{fig:simulation:normal}
\end{figure}


We generate data as follows. A random graph with $p$ nodes is created
by first partitioning nodes into $p/20$ connected components, each
with $20$ nodes, and then forming a random graph over these $20$
nodes. A chain graph is formed by permuting the nodes and connecting
them in succession, while a nearest-neighbor graph is constructed
following the procedure outlined in \citet{li06gradient}. That is, for
each node, we draw a point uniformly at random on a unit square and
compute the pairwise distances between nodes. Each node is then
connected to $s=4$ closest neighbors.  Since some of nodes will have
more than $4$ adjacent edges, we  randomly remove edges from nodes that
have degree larger than $4$ until the maximum degree of a node in a
network is $4$. Once the graph is created, we construct a
precision matrix, with non-zero blocks corresponding to edges in the
graph. Elements of diagonal blocks are set as $0.5^{|a-b|}$, $0\leq
a,b\leq k$, while off-diagonal blocks have elements with the same
value, $0.2$ for chain graphs and $0.3/k$ for nearest-neighbor
networks. Finally, we add $\rho I$ to the precision matrix, so that
its minimum eigenvalue is equal to $0.5$.  Note that $s=2$ for the
chain graph and $s=4$ for the nearest-neighbor graph. Simulation
results are averaged over 100 replicates.

Figure~\ref{fig:simulation:normal} shows simulation results.  Each row
in the figure reports results for one method, while each column in the
figure represents a different simulation setting. For the first two
columns, we set $k=3$ and vary the total number of nodes in the
graph. The third simulation setting sets the total number of nodes
$p=20$ and changes the number of attributes $k$. In the case of the
chain graph, we observe that for small sample sizes the method of
\cite{danaher2011joint} outperforms all the other methods. We note
that the multi-attribute method is estimating many more parameters,
which require large sample size in order to achieve high
accuracy. However, as the sample size increases, we observe that
multi-attribute method starts to outperform the other methods. In
particular, for the sample size indexed by $\theta = 13$ all the graph
are correctly recovered, while other methods fail to recover the
graph consistently at the same sample size. In the case of
nearest-neighbor graph, none of the methods recover the graph well for
small sample sizes. However, for moderate sample sizes,
multi-attribute method outperforms the other methods.  Furthermore,
as the sample size increases none of the other methods recover the
graph exactly. This suggests that the conditions for consistent
graph recovery may be weaker in the multi-attribute setting.


\begin{figure}[p]
  \centering

  \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoEasy1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoEasy2}
    \end{subfigure}%
    \caption{glasso procedure}
    \label{fig:easy:glasso}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskEasy1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskEasy2}
    \end{subfigure}%
    \caption{Procedure of \cite{danaher2011joint}}
    \label{fig:easy:mt}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoEasy1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoEasy2}
    \end{subfigure}%
    \caption{Procedure of \cite{Guo:09}}
    \label{fig:easy:guo}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttEasy1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttEasy2}
    \end{subfigure}%
    \caption{Multi-attribute procedure}
    \label{fig:easy:ma}
  \end{subfigure
  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent
    runs. Blocks $\Omegab_{ab}$ of the precision matrix $\Omegab$ are
    diagonal matrices. }
  \label{fig:simulation:easy}
\end{figure}

\begin{figure}[p]
  \centering

  \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoHard1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoHard2}
    \end{subfigure}%
    \caption{glasso procedure}
    \label{fig:hard:glasso}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskHard1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskHard2}
    \end{subfigure}%
    \caption{Procedure of \cite{danaher2011joint}}
    \label{fig:hard:mt}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoHard1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoHard2}
    \end{subfigure}%
    \caption{Procedure of \cite{Guo:09}}
    \label{fig:hard:guo}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttHard1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttHard2}
    \end{subfigure}%
    \caption{Multi-attribute procedure}
    \label{fig:hard:ma}
  \end{subfigure}%
  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent runs.
    Off-diagonal blocks $\Omegab_{ab}$ of the precision matrix
    $\Omegab$ have zeros as diagonal elements.}
  \label{fig:simulation:hard}
\end{figure}

\begin{figure}[p]
  \centering

  \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoRS1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoRS2}
    \end{subfigure}%
    \caption{glasso procedure}
    \label{fig:rs::glasso}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskRS1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskRS2}
    \end{subfigure}%
    \caption{Procedure of \cite{danaher2011joint}}
    \label{fig:rs::mt}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoRS1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoRS2}
    \end{subfigure}%
    \caption{Procedure of \cite{Guo:09}}
    \label{fig:rs::guo}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttRS1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttRS2}
    \end{subfigure}%
    \caption{Multi-attribute procedure}
    \label{fig:rs::ma}
  \end{subfigure}%
  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent runs.
    Off-diagonal blocks $\Omegab_{ab}$ of the precision matrix
    $\Omegab$ have elements uniformly sampled from $[-0.3, -0.1]
    \bigcup [0.1,0.3]$.  }
  \label{fig:simulation:RS}
\end{figure}



\subsection{Alternative Structure of Off-diagonal Blocks}

In this section, we investigate performance of different estimation
procedures under different assumptions on the elements of the
off-diagonal blocks of the precision matrix.

First, we investigate a situation where the multi-attribute method does
not perform as well as the methods that estimate multiple graphical
models. One such situation arises when different attributes are
conditionally independent. To simulate this situation, we use the data
generating approach as before, however, we make each block
$\Omegab_{ab}$ of the precision matrix $\Omega$ a diagonal
matrix. Figure \ref{fig:simulation:easy} summarizes results of the
simulation. We see that the methods of \cite{danaher2011joint} and
\cite{Guo:09} perform better, since they are estimating much fewer
parameters than the multi-attribute method. ${\rm glasso}$ does not
utilize any structural information underlying the estimation problem
and requires larger sample size to correctly estimate the graph than
other methods.

A completely different situation arises when the edges between nodes
can be inferred only based on inter-attribute data, that is, when a
graph based on any individual attribute is empty. To generate data
under this situation, we follow the procedure as before, but with the
diagonal elements of the off-diagonal blocks $\Omegab_{ab}$ set to
zero. Figure~\ref{fig:simulation:hard} summarizes results of the
simulation. In this setting, we clearly see the advantage of the
multi-attribute method, compared to other three methods. Furthermore,
we can see that ${\rm glasso}$ does better than multi-graph methods of
\cite{danaher2011joint} and \cite{Guo:09}. The reason is that ${\rm
  glasso}$ can identify edges based on inter-attribute relationships
among nodes, while multi-graph methods rely only on intra-attribute
relationships. This simulation illustrates an extreme scenario where
inter-attribute relationships are important for identifying edges.

So far, off-diagonal blocks of the precision matrix were constructed
to have constant values. Now, we use the same data generating
procedure, but generate off-diagonal blocks of a precision matrix in a
different way. Each element of the off-diagonal block $\Omega_{ab}$ is
generated independently and uniformly from the set $[-0.3, -0.1]
\bigcup [0.1,0.3]$. The results of the simulation are given in
Figure~\ref{fig:simulation:RS}. Again, qualitatively, the results are
similar to those given in Figure~\ref{fig:simulation:normal}, except
that in this setting more samples are needed to recover the graph
correctly.

\subsection{Different Number of Samples per Attribute}

In this section, we show how to deal with a case when different number
of samples is available per attribute. As noted in
$\S$\ref{sec:estimation}, we can treat non-measured attributes as
missing completely at random \citep[see][for more
details]{kolar12missing}. 

Let $R = (r_{il})_{i\in\{1,\ldots,n\},l\in\{1,\ldots,pk\}}\in \RR^{n \times pk}$ be an
indicator matrix, which denotes for each sample point $x_i$ the components
that are observed. Then the sample covariance matrix $S =
(\sigma_{lk}) \in \RR^{pk \times pk}$ is estimated as
$\sigma_{lk} = \rbr{\sum_{i=1}^n r_{i,l}r_{i,k}}^{-1}
{\sum_{i=1}^n r_{i,l}r_{i,k}x_{i,l}x_{i,k}}.$
This estimate is plugged into the objective in \eqref{eq:max_ll_opt}.

\begin{figure}[b]
  \centering

  \includegraphics[width=\textwidth]{plots/biometrikaDiffSample}

  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent
    runs. Additional samples available for the first attribute. }
  \label{fig:simulation:diff_sample}
\end{figure}

We generate a chain graph with $p = 60$ nodes, construct a precision
matrix associated with the graph and $k=3$ attributes, and generate $n
= \theta s^2k^2\log(pk)$ samples, $\theta > 0$. Next, we generate
additional $10\%$, $30\%$ and $50\%$ samples from the same model, but
record only the values for the first attribute. Results of the
simulation are given in Figure~\ref{fig:simulation:diff_sample}.
Qualitatively, the results are similar to those presented in
Figure~\ref{fig:simulation:normal}.


\section{Illustrative Applications to Real Data}

In this section, we illustrate how to apply our method to data arising
in studies of biological regulatory networks and Alzheimer's disease.

\subsection{Analysis of a Gene/Protein Regulatory Network}
\label{sec:real_data:sup}


We provide illustrative, exploratory analysis of data from the
well-known NCI-60 database, which contains different molecular
profiles on a panel of 60 diverse human cancer cell
lines.
Data set consists of protein profiles (normalized reverse-phase lysate
arrays for 92 antibodies) and gene profiles (normalized RNA microarray
intensities from Human Genome U95 Affymetrix chip-set for $> 9000$
genes). We focus our analysis on a subset of 91 genes/proteins for
which both types of profiles are available. These profiles are
available across the same set of $60$ cancer cells. More detailed
description of the data set can be found in \citet{katenka2011multi}.

We inferred three types of networks: a network based on protein
measurements alone, a network based on gene expression profiles and a
single gene/protein network. For protein and gene networks we use the
\texttt{glasso}, while for the gene/protein network, we use our
procedure outlined in $\S$\ref{sec:estimation}. We use the stability
selection~\citep{stability:10} procedure to estimate stable
networks. In particular, we first select the penalty parameter
$\lambda$ using cross-validation, which over-selects the number of
edges in a network. Next, we use the selected $\lambda$ to estimate
100 networks based on random subsamples containing 80\% of the
data-points. Final network is composed of stable edges that appear in
at least 95 of the estimated networks. Table~\ref{tab:net_stats}
provides a few summary statistics for the estimated
networks. Furthermore, protein and gene/protein networks share $96$
edges, while gene and gene/protein networks share $104$ edges. Gene
and protein network share only $17$ edges. Finally, $66$ edges are
unique to gene/protein network. Figure~\ref{fig:node_degree} shows
node degree distributions for the three networks.  We observe that the
estimated networks are much sparser than the association networks in
\citet{katenka2011multi}, as expected due to marginal correlations
between a number of nodes.  The differences in networks require a
closer biological inspection by a domain scientist.


\begin{table}[b]

 \caption{Summary statistics for protein, gene, and gene/protein
   networks ($p = 91$).}

 \begin{tabular}{lrrr}
 ~ & { protein network} & {gene network} & {gene/protein network} \\ 

Number of edges & 122&   214 &   249 \\
Density & 0.03 & 0.05 & 0.06 \\
Largest connected component & 62 & 89 & 82 \\
Avg Node Degree  & 2.68 &  4.70 & 5.47\\
Avg Clustering Coefficient  & 0.0008 &    0.001 &   0.003\\
 \end{tabular}

\label{tab:net_stats}
\end{table}         


\begin{figure}[t]
  \centering
  \includegraphics[width=\columnwidth]{nodeDegree}
  \caption{Node degree distributions for protein, gene and
    gene/protein networks.}
  \label{fig:node_degree}
\end{figure}

\begin{figure}[t]
  \centering
  \includegraphics[width=\columnwidth]{nodeClass}
 
 
 
 
 
 
 
 
 
 
 
 
 
 
\caption{Edge and node classification based on $w_p^2$.}
  \label{fig:node_class}
\end{figure}


We proceed with a further exploratory analysis of the gene/protein
network. We investigate the contribution of two nodal attributes to the 
existence of an edges between the nodes. Following
\citet{katenka2011multi}, we use a simple heuristic based on the weight
vectors to classify the nodes and edges into three classes. For an
edge between the nodes $a$ and $b$, we take one weight vector, say
$\wb_a$, and normalize it to have unit norm. Denote $w_{p}$ the
component corresponding to the protein attribute. Left plot in
Figure~\ref{fig:node_class} shows the values of $w_p^2$ over all
edges. The edges can be classified into three classes based on the
value of $w_p^2$. Given a threshold $T$, the edges for which $w_p^2
\in (0, T)$ are classified as gene-influenced, the edges for which
$w_p^2 \in (1-T, 1)$ are classified as protein influenced, while the
remainder of the edges are classified as mixed type. In the left plot
of Figure~\ref{fig:node_class}, the threshold is set as
$T=0.25$. Similar classification can be performed for nodes after
computing the proportion of incident edges. Let $p_1$, $p_2$ and $p_3$
denote proportions of gene, protein and mixed edges, respectively,
incident with a node. These proportions are represented in a simplex
in the right subplot of Figure~\ref{fig:node_class}. Nodes with mostly
gene edges are located in the lower left corner, while the nodes with
mostly protein edges are located in the lower right corner. Mixed
nodes are located in the center and towards the top corner of the
simplex. Further biological enrichment analysis is possible (see
\citet{katenka2011multi}), however, we do not pursue this here.


\subsection{Uncovering Functional Brain Network}
\label{sec:real_data}

We apply our method to the Positron Emission Tomography dataset, which
contains 259 subjects, of whom 72 are healthy, 132 have mild cognitive
Impairment and 55 are diagnosed as Alzheimer's \& Dementia.  Note that
mild cognitive impairment is a transition stage from normal aging to
Alzheimer's \& Dementia. The data can be obtained from {
 \it http://adni.loni.ucla.edu/}.  The preprocessing is done in the same
way as in \cite{huang2009learning}.

\begin{figure}[!b]
  \centering

    \begin{subfigure}[b]{0.33\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/normGraph1}
      \caption{Healthy subjects}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.33\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/mciGraph1}
      \caption{Mild Cognitive Impairment}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.33\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/adGraph1}
      \caption{Alzheimer's \& Dementia}
    \end{subfigure}%
  
  \caption{Brain connectivity networks}
  \label{fig:brain}
\end{figure}


Each Positron Emission Tomography image contains $91 \times 109 \times
91 = 902,629$ voxels. The effective brain region contains $180,502$
voxels, which are partitioned into $95$ regions, ignoring the regions
with fewer than $500$ voxels. The largest region contains $5,014$
voxels and the smallest region contains $665$ voxels. Our
preprocessing stage extracts $948$ representative voxels from these
regions using the $K$-median clustering algorithm. The parameter $K$ is
chosen differently for each region, proportionally to the initial
number of voxels in that region.  More specifically, for each category
of subjects we have an $n \times (d_1 + \ldots + d_{95})$ matrix, where
$n$ is the number of subjects and $d_1 + \ldots + d_{95} = 902,629$ is
the number of voxels. Next we set $K_i = \lceil d_i / \sum_j d_j
\rceil$, the number of representative voxels in region $i$,
$i=1,\ldots,95$. The representative voxels are identified by running
the $K$-median clustering algorithm on a sub-matrix of size $n \times
d_i$ with $K = K_i$.

We inferred three networks, one for each subtype of subjects using the
procedure outlined in $\S$\ref{sec:estimation}.  Note that for
different nodes we have different number of attributes, which
correspond to medians found by the clustering algorithm.  We use the
stability selection~\citep{stability:10} approach to estimate stable
networks.  The stability selection procedure is combined with our
estimation procedure as follows.  We first select the penalty
parameter $\lambda$ in \eqref{eq:max_ll_opt} using cross-validation,
which overselects the number of edges in a network.  Next, we create
$100$ subsampled data sets, each of which contain $80\%$ of the data
points, and estimate one network for each dataset using the selected
$\lambda$.  The final network is composed of stable edges that appear
in at least $95$ of the estimated networks.

We visualize the estimated networks in Figure~\ref{fig:brain}.
Table~\ref{tab:brain_net_stats} provides a few summary statistics for
the estimated networks. Appendix~\ref{appendix:d} contains names of
different regions, as well as the adjacency matrices for networks.
From the summary statistics, we can observe that in normal subjects
there are many more connections between different regions of the
brain.  Loss of connectivity in Alzheimer's \& Dementia has been
widely reported in the literature \citep{greicius2004default,
  hedden2009disruption, andrews2007disruption, wu2011altered}.

Learning functional brain connectivity is potentially valuable for
early identification of signs of Alzheimer's disease.
\citet{huang2009learning} approach this problem using exploratory data
analysis. The framework of Gaussian graphical models is
used to explore functional brain connectivity. Here we point out that
our approach can be used for the same exploratory task, without the
need to reduce the information in the whole brain to one number. For
example, from our estimates, we observe the loss of connectivity in
the cerebellum region of patients with Alzheimer's disease, which has
been reported previously in \citet{sjobeck2001alzheimer}. As another
example, we note increased connectivity between the frontal lobe and
other regions in the patients, which was linked to compensation for
the lost connections in other regions \citep{stern2006cognitive,
  gould2006brain}. 

\begin{table}[t]

 \caption{Summary statistics for protein, gene, and gene/protein
   networks ($p = 91$)}

 \begin{tabular}{lrrr}
 ~ & { Healthy }
 & {Mild Cognitive} 
& {Alzheimer's \& } \\ 

& subjects & Impairment & Dementia \\

Number of edges & 116 & 84  &   59 \\
Density & 0.030 & 0.020 & 0.014 \\
Largest connected component & 48 & 27 & 25 \\
Avg Node Degree  & 2.40 & 1.73 &  1.2 \\
Avg Clustering Coefficient  & 0.001  & 0.0023 &    0.0007 \\
 \end{tabular}

\label{tab:brain_net_stats}
\end{table}         


\subsection*{Acknowledgments}

{\small
We thank Eric D. Kolaczyk and Natallia V. Katenka for sharing
preprocessed data used in their study with us. Eric P. Xing is
partially supported through the grants NIH R01GM087694 and AFOSR
FA9550010247.  The research of Han Liu is supported by NSF grant
IIS-1116730.
 }


